text:
###############################################################################
# Watched Expressions / Snapall
###############################################################################
export @watch, @unwatch, @unwatchall, @snapall

set_autowatch(on::Bool) = global autowatch = on
get_autowatch() = autowatch

watched_exprstrs() = keys(watched_exprs)

watch_exprstr(exprstr) = begin
    watched_exprs[exprstr] = true
end

unwatch_exprstr(exprstr) = begin
    delete!(watched_exprs, exprstr)
end

macro watch(exprs...)
    for expr in exprs
        watch_exprstr(string(expr)) #called at macro expansion time, not run time
    end
    :()
end

macro unwatch(exprs...)
    for expr in exprs
        unwatch_exprstr(string(expr))
    end
    :()
end

macro unwatchall()
    empty!(watched_exprs)
    :()
end

macro snapall(kwexprs...)
    kwargs, kwextras, extra_exprs, arginfo = kwargparse(trace_kwargspec, kwexprs)
    watched_exprs = map(parse, watched_exprstrs())
    exprs = vcat(watched_exprs, extra_exprs)
    condition = kwargs[:iff]
    quote
        if $condition
            @RickTracy.snapexprs $(kwargs[:location]) $(kwargs[:everyN]) $exprs
        end
    end |> esc
end

max_stars_repo_name: JobJob/RickTracy


text:
using PackageCompiler

tutorials_file = joinpath("..", "12_n_site.jl")
PackageCompiler.create_sysimage(
  ["ITensors", "ITensorUnicodePlots", "ITensorGLMakie", "Zygote"];
  sysimage_path="ITensor.so",
  precompile_execution_file=tutorials_file,
)

# Then make an alias like:
#
# alias julia-itensor='julia -J/home/mfishman/.julia/dev/ITensorTutorials/tutorials/package_compile/ITensor.so'
#
# and run `julia-itensor`.

ext: jl


text:
using Distributions
using Klara
using PGUManifoldMC

include("rv_model.jl")

using RvModelKeplerian

DATADIR = "../../data"
SUBDATADIR = "alsmmala"

dataset = readdlm(joinpath(DATADIR, "example2.csv"), ',', header=false); # read observational data
obs_times = dataset[:,1]
obs_rv = dataset[:,2]
sigma_obs = dataset[:,3]
set_times(obs_times);     # set data to use for model evaluation
set_obs( obs_rv);
set_sigma_obs(sigma_obs);

include("utils_ex.jl")
param_true = make_param_true_ex2()
param_perturb_scale = make_param_perturb_scale(param_true)
param_init = 0
param_init = param_true.+0.01*param_perturb_scale.*randn(length(param_true))
println("param_init= ",param_init)

sampler = ALSMMALA(
  # 0.05,
  0.03,
  identitymala=false,
  #update=(sstate, pstate, i, tot) -> rand_exp_decay_update!(sstate, pstate, i, tot, 10., 0.15),
  # update=(sstate, pstate, i, tot) -> mod_update!(sstate, pstate, i, tot, 50),
  update=(sstate, pstate, i, tot) -> rand_exp_decay_update!(sstate, pstate, i, tot, 7.),
  transform=H -> softabs(H, 1000.),
  initupdatetensor=(true, false)
)

nmcmc = 50000
nburnin = 10000
mcrange = BasicMCRange(nsteps=nmcmc, burnin=nburnin)

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

v0 = Dict(:p=>param_init)

p = BasicContMuvParameter(
  :p,
  logtarget=plogtarget,
  diffopts=DiffOptions(mode=:forward, order=2)
)

model = likelihood_model(p, false)
#model = likelihood_model([p], isindexed=false)

#target_accept_rates = [0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.93, 0.95, 0.97, 0.99]
#target_accept_rates = [0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4] #, 0.45 0.5] #, 0.75, 0.8, 0.85, 0.9, 0.93, 0.95, 0.97, 0.99]
target_accept_rates = [0.7]
alsmmala_times = Array(Float64, length(target_accept_rates))
alsmmala_stepsizes = Array(Float64, length(target_accept_rates))
alsmmala_acceptance = Array(Float64, length(target_accept_rates))
alsmmala_esses = Array(Float64, length(target_accept_rates))
alsmmala_iacts = Array(Float64, length(target_accept_rates))
alsmmala_chains = Array(Any, length(target_accept_rates))

for i in 1:length(target_accept_rates)
  target_accept_rate = target_accept_rates[i]
  println("# i= ", i, ": ",target_accept_rate)
  mctuner = PSMMALAMCTuner(
    VanillaMCTuner(verbose=true),
    VanillaMCTuner(verbose=true),
    VanillaMCTuner(verbose=true)
    # AcceptanceRateMCTuner(target_accept_rate, verbose=true)
  )

  job = BasicMCJob(
    model,
    sampler,
    mcrange,
    v0,
    tuner=mctuner,
    outopts=outopts
  )

  tic()
  run(job)
  runtime = toc()

  chain = output(job)
  ratio = acceptance(chain)

  println("# mean[",i,"] = ",mean(chain))
  println("# std[",i,"] =",vec(std(chain.value,2)))
  println("# ess[",i,"] = ",ess(chain))

  alsmmala_times[i] = runtime
  alsmmala_stepsizes[i] = job.sstate.tune.totaltune.step
  alsmmala_acceptance[i] = ratio
  alsmmala_esses[i] = minimum(ess(chain))
  alsmmala_iacts[i] = maximum(iact(chain))
  alsmmala_chains[i] = chain
end


text:
using MLKernels

SUITE["MLKernels"] = BenchmarkGroup()

mlkernelnames = ["SquaredExponentialKernel"]
kernels = Dict{String,MLKernels.Kernel}()
for k in mlkernelnames
    SUITE["MLKernels"][k] = BenchmarkGroup()
    kernels[k] = eval(Meta.parse("MLKernels." * k * "(alpha)"))
end

for k in mlkernelnames
    SUITE["MLKernels"][k]["k(X,Y)"] = @benchmarkable MLKernels.kernelmatrix(
        $(kernels[k]), $X, $Y
    )
    # SUITE["MLKernels"][k][kt]["k!(X,Y)"] = @benchmarkable MLKernels.kernelmatrix!(KXY,$(kernels[k][kt]),$X,$Y) setup=(KXY=copy($KXY))
    SUITE["MLKernels"][k]["k(X)"] = @benchmarkable MLKernels.kernelmatrix($(kernels[k]), $X)
    # SUITE["MLKernels"][k][kt]["k!(X)"] = @benchmarkable MLKernels.kernelmatrix!(KX,$(kernels[k][kt]),$X) setup=(KX=copy($KX))
    # SUITE["MLKernels"][k][kt]["kdiag(X)"] = @benchmarkable MLKernels.kernelmatrix_diag($(kernels[k][kt]),$X)
    # SUITE["MLKernels"][k][kt]["kdiag!(X)"] = @benchmarkable MLKernels.kernelmatrix_diag!(kX,$(kernels[k][kt]),$X) setup=(kX=copy($kX))
end
# results = run(SUITE)


text:
function aria!(ec::EChart; show::Bool = true, description::Union{String, Nothing} = nothing, general::Union{Dict, Nothing} = nothing, series::Union{Dict, Nothing} = nothing, data::Union{Dict, Nothing} = nothing)

    ec.aria = Aria(show = show, description = description, general = general, series = series, data = data)
    return ec

end


text:
include("../LF_common/m_LF1d.jl")
include("../LF_common/m_LF3d.jl")
include("../LF_common/m_Gvectors.jl")

include("../utils/ortho_gram_schmidt.jl")
include("../utils/orthonormalize.jl")

include("init_pot_harm_3d.jl")
include("init_pot_Hps.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")

include("op_nabla2.jl")
include("op_H.jl")

include("calc_Energies.jl")
include("calc_evals.jl")
include("calc_rho.jl")
include("calc_grad.jl")

include("KS_solve_Emin_cg.jl")
include("KS_solve_Emin_pcg.jl")
include("KS_solve_scf.jl")

include("../LF_common/sparse_LF3d.jl")
include("../LF_common/prec_mkl_ilu0.jl")
include("../LF_common/apply_prec_ilu0.jl")

include("../LF_common/solve_poisson_FFT.jl")

include("LDA_VWN.jl")
include("diag_lobpcg.jl")

function test_main( ; pot_type = "harmonic", method = "Emin_cg" )
    # LF parameters
    NN = [25, 25, 25]
    AA = [0.0, 0.0, 0.0]
    BB = [6.0, 6.0, 6.0]

    Npoints = prod(NN)

    # Initialize LF
    LF = init_LF3d_p( NN, AA, BB )
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

    # Initialize G-vectors
    L = BB - AA
    Gv = GvectorsT( NN, diagm(L) )

    if pot_type == "harmonic"
        #
        println("\nHarmonic potential\n")
        #
        # Parameter for potential
        center = AA + 0.5*(BB-AA)
        # Potential
        ω = 2.0
        V_ionic = init_pot_harm_3d( LF, ω, center )
        Nstates = 4
        Focc = 2.0*ones(Nstates)
    elseif pot_type == "Hps"
        #
        println("\nHps HGH G\n")
        #
        # center of the potential
        center = 0.5*(BB-AA)
        # Potential
        V_ionic = init_pot_Hps_HGH_G( Gv, center )
        println("sum(V_ionic): ", sum(V_ionic))
        #
        Nstates = 1
        Focc = 1.0*ones(Nstates)
    elseif pot_type == "Hcoul"
        #
        println("\nHcoul in G\n")
        #
        #
        Vg = zeros(Complex128,Npoints)
        for ig=2:Npoints
            Vg[ig] = -4*pi/Gv.G2[ig]/Gv.Ω
        end
        V_ionic = real( G_to_R(NN,Vg) ) * Npoints
        #
        println("sum(V_ionic): ", sum(V_ionic))
        #
        Nstates = 1
        Focc = 1.0*ones(Nstates)
    end

    if method == "Emin_cg_sparse"
        #
        ∇2 = get_Laplacian3d_kron(LF)
        precH = prec_mkl_ilu0( -0.5*∇2 + spdiagm(V_ionic) )
        #precH = speye(Npoints)  # for testing unpreconditioned code
        Energies, evecs, Potentials = KS_solve_Emin_pcg( LF, Gv, ∇2, precH,
                                        V_ionic, Focc, Nstates, verbose=true )
        evals = calc_evals( LF, ∇2, Potentials, evecs )
        #
    elseif method == "SCF"
        ∇2 = get_Laplacian3d_kron(LF)
        precH = prec_mkl_ilu0( -0.5*∇2 + spdiagm(V_ionic) )
        #precH = speye(Npoints)
        Energies, evals, evecs, Potentials = KS_solve_scf( LF, Gv, ∇2, precH,
                                               V_ionic, Focc, Nstates, verbose=true )
        #
    else
    #
        Energies, evecs, Potentials = KS_solve_Emin_cg( LF, Gv, V_ionic, Focc, Nstates, verbose=true )
        evals = calc_evals( LF, Potentials, evecs )
    end

    println("")
    print_Energies(Energies)
    @printf("\nEigenvalues:\n")
    for i = 1:Nstates
        @printf("%8d %f\n", i, evals[i])
    end

    rho = calc_rho( Focc, evecs )
    @printf( "\nintegRho = %lf\n", sum(rho)*ΔV )

end

#@code_native test_main()
#@time test_main(pot_type="harmonic", method="Emin_cg")
#@time test_main(pot_type="harmonic", method="Emin_cg_sparse")
#@time test_main(pot_type="harmonic", method="SCF")
#@time test_main(pot_type="Hps", method="SCF")
@time test_main(pot_type="Hcoul", method="SCF")


text:
module Precondition

using CUDA
using CUDA.CUSPARSE
using KernelAbstractions
using LightGraphs
using LinearAlgebra
using Metis
using SparseArrays
using TimerOutputs

cuzeros = CUDA.zeros

abstract type AbstractPreconditioner end

struct NoPreconditioner <: AbstractPreconditioner end

mutable struct Preconditioner <: AbstractPreconditioner
    npart::Int64
    nJs::Int64
    partitions::Vector{Vector{Int64}}
    cupartitions
    Js
    cuJs
    map
    cumap
    part
    cupart
    P
    function Preconditioner(J, npart, device=CPU())
        m, n = size(J)
        if npart < 2
            error("Number of partitions `npart` should be at" *
                  "least 2 for partitioning in Metis")
        end
        adj = build_adjmatrix(J)
        g = Graph(adj)
        part = Metis.partition(g, npart)
        partitions = Vector{Vector{Int64}}()
        for i in 1:npart
            push!(partitions, [])
        end
        for (i,v) in enumerate(part)
            push!(partitions[v], i)
        end
        Js = Vector{Matrix{Float64}}(undef, npart)
        nJs = maximum(length.(partitions))
        id = Matrix{Float64}(I, nJs, nJs)
        println("Block Jacobi block size: $nJs")
        for i in 1:npart
            Js[i] = Matrix{Float64}(I, nJs, nJs)
        end
        nmap = 0
        for b in partitions
            nmap += length(b)
        end
        map = Vector{Int64}(undef, nmap)
        part = Vector{Int64}(undef, nmap)
        for b in 1:npart
            for (i,el) in enumerate(partitions[b])
                map[el] = i
                part[el] = b
            end
        end
        row = Vector{Float64}()
        col = Vector{Float64}()
        nzval = Vector{Float64}()

        for b in 1:npart
            for x in partitions[b]
                for y in partitions[b]
                    push!(row, x)
                    push!(col, y)
                    push!(nzval, 1.0)
                end
            end
        end
        P = sparse(row, col, nzval)
        if isa(device, CUDADevice)
            cupartitions = Vector{CuVector{Int64}}(undef, npart)
            for i in 1:npart
                cupartitions[i] = CuVector{Int64}(partitions[i])
            end
            cuJs = Vector{CuMatrix{Float64}}(undef, length(partitions))
            for i in 1:length(partitions)
                cuJs[i] = CuMatrix{Float64}(I, nJs, nJs)
            end
            cumap = cu(map)
            cupart = cu(part)
            P = CuSparseMatrixCSR(P)
        else
            cuJs = nothing
            cupartitions = nothing
            cumap = nothing
            cupart = nothing
        end
        return new(npart, nJs, partitions, cupartitions, Js, cuJs, map, cumap, part, cupart, P)
    end
end

function build_adjmatrix(A)
    rows = Int64[]
    cols = Int64[]
    vals = Float64[]
    @show size(A)
    rowsA = rowvals(A)
    m, n = size(A)
    for i = 1:n
        for j in nzrange(A, i)
            push!(rows, rowsA[j])
            push!(cols, i)
            push!(vals, 1.0)

            push!(rows, i)
            push!(cols, rowsA[j])
            push!(vals, 1.0)
        end
    end
    println("Creating matrix")
    return sparse(rows,cols,vals,size(A,1),size(A,2))
end

function fillblock_gpu!(cuJs, partition, map, rowPtr, colVal, nzVal, part, b)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(partition)
        for j in rowPtr[partition[i]]:rowPtr[partition[i]+1]-1
            if b == part[colVal[j]]
                @inbounds cuJs[map[partition[i]], map[colVal[j]]] = nzVal[j]
            end
        end
    end
    return nothing
end

function fillP_gpu!(cuJs, partition, map, rowPtr, colVal, nzVal, part, b)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(partition)
        for j in rowPtr[partition[i]]:rowPtr[partition[i]+1]-1
            if b == part[colVal[j]]
                @inbounds nzVal[j] += cuJs[map[partition[i]], map[colVal[j]]]
            end
        end
    end
    return nothing
end

function update(J::CuSparseMatrixCSR, p, to)
    m = size(J, 1)
    n = size(J, 2)
    nblocks = length(p.partitions)
    @timeit to "Fill Block Jacobi" begin
        CUDA.@sync begin
            for b in 1:nblocks
                @cuda threads=16 blocks=16 fillblock_gpu!(p.cuJs[b], p.cupartitions[b], p.cumap, J.rowPtr, J.colVal, J.nzVal, p.cupart, b)
            end
        end
    end
    @timeit to "Invert blocks" begin
        CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(p.cuJs, true)
        CUDA.@sync pivot, info, p.cuJs = CUDA.CUBLAS.getri_batched(p.cuJs, pivot)
    end
    p.P.nzVal .= 0.0
    @timeit to "Move blocks to P" begin
        CUDA.@sync begin
            for b in 1:nblocks
                @cuda threads=16 blocks=16 fillP_gpu!(p.cuJs[b], p.cupartitions[b], p.cumap, p.P.rowPtr, p.P.colVal, p.P.nzVal, p.cupart, b)
            end
        end
    end
    return p.P
end

function update(J::SparseMatrixCSC, p, to)
    nblocks = length(p.partitions)
    @timeit to "Fill Block Jacobi" begin
        for b in 1:nblocks
            for i in p.partitions[b]
                for j in J.colptr[i]:J.colptr[i+1]-1
                    if b == p.part[J.rowval[j]]
                        p.Js[b][p.map[J.rowval[j]], p.map[i]] = J.nzval[j]
                    end
                end
            end
        end
    end
    @timeit to "Invert blocks" begin
        for b in 1:nblocks
            p.Js[b] = inv(p.Js[b])
        end
    end
    p.P.nzval .= 0.0
    @timeit to "Move blocks to P" begin
        for b in 1:nblocks
            for i in p.partitions[b]
                for j in p.P.colptr[i]:p.P.colptr[i+1]-1
                    if b == p.part[p.P.rowval[j]]
                        p.P.nzval[j] += p.Js[b][p.map[p.P.rowval[j]], p.map[i]]
                    end
                end
            end
        end
    end
    return p.P
end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

abstract type AbstractRayGenerator{T<:Real} end

"""
    GeometricRayGenerator{T,O<:RayOriginGenerator{T}} <: AbstractRayGenerator{T}

Generates geometric [`Ray`](@ref)s according to the specific implementation of the subclass.
"""
abstract type GeometricRayGenerator{T} <: AbstractRayGenerator{T} end

"""
    OpticalRayGenerator{T} <: AbstractRayGenerator{T}

Generates [`OpticalRay`](@ref)s according to the specific implementation of the subclass.
"""
abstract type OpticalRayGenerator{T} <: AbstractRayGenerator{T} end
export AbstractRayGenerator, GeometricRayGenerator, OpticalRayGenerator


text:
cd(@__DIR__); include("setups/path512.jl")
pyplot(dpi = 200)

## construct full HGLET and Lapped-HGLET dictionaries
@time HGLET_dic = HGLET_dictionary(GP, G_Sig; gltype = :L)
@time LPHGLET_dic = LPHGLET_dictionary(GP, G_Sig; gltype = :L, ϵ = 0.3)

j = 3; k = 2; l = 6;
WH = HGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
WlH = LPHGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'

plt = plot(WH[:, l], c = :black, grid = false, frame = :box, lw = 0.8,
     legendfontsize = 11, legend = false, size = (500, 400), ylim = [-0.15, 0.15])
     xticks!([1; 64:64:N], vcat(string(1), [string(k) for k in 64:64:N]))
savefig(plt, "../figs/Path512_HGLET_j$(j-1)k$(k-1)l$(l-1).png")

plt = plot(WlH[:, l], c = :black, grid = false, frame = :box, lw = 0.8,
     legendfontsize = 11, legend = false, size = (500, 400), ylim = [-0.15, 0.15])
     xticks!([1; 64:64:N], vcat(string(1), [string(k) for k in 64:64:N]))
savefig(plt, "../figs/Path512_LPHGLET_j$(j-1)k$(k-1)l$(l-1).png")


text:
using Optim, Plots
# method for contour
rosenbrock(x1,x2) = (1.0 - x1)^2 + 100.0 * (x2 - x1^2)^2
# method for optim
rosenbrock(θ) = (1.0 - θ[1])^2 + 100.0 * (θ[2] - θ[1]^2)^2
θ = zeros(2)
b1 = range(0.0,stop=1.5,length=100)
b2 = range(0.0,stop=1.5,length=100)
contour(b1,b2,(b1,b2)->rosenbrock(b1,b2),fill=true, c=:viridis)
θs = zeros(100,2)
for i = 1:100
    #global θ = optimize(rosenbrock, θ, GradientDescent(), Optim.Options(iterations=1)).minimizer
    global θ = optimize(rosenbrock, θ, Newton(), Optim.Options(iterations=1)).minimizer
    θs[i,:] = θ'
    display(scatter!(θs[1:i,1], θs[1:i,2], legend=false))
    sleep(0.1)
end
θs

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module layout

import Parameters: @with_kw

import Term: int, get_lr_widths, textlen
import ..renderables: RenderablesUnion, Renderable, AbstractRenderable, RenderableText
import ..style: apply_style
import ..measure: Measure
import ..segment: Segment
using ..box
import ..box: get_lrow, get_rrow
import ..consoles: console_width, console_height

export Padding, vstack, hstack, pad
export Spacer, vLine, hLine

# ---------------------------------------------------------------------------- #
#                                    PADDING                                   #
# ---------------------------------------------------------------------------- #
"""
    Padding

Stores information about ammount of padding.
"""
@with_kw struct Padding
    left::Int
    right::Int
    top::Int
    bottom::Int
end

"""
    pad(text::AbstractString, target_width::Int, method::Symbol)::String

Pad a string to width: `target_width` by adding empty spaces strings " ".
Where the spaces are added depends on the justification `method` ∈ (:left, :center, :right).
"""
function pad(text, target_width::Int, method::Symbol)::String
    # get total padding size
    lw = textlen(text)
    lw >= target_width && return text

    npads = target_width - lw
    if method == :left
        return text * " "^npads
    elseif method == :right
        return " "^npads * text
    else
        nl, nr = get_lr_widths(npads)
        return " "^nl * text * " "^nr
    end
end

"""
    pad(text::AbstractString, left::Int=0, right::Int=0)::String

Pad a string by a fixed ammount.
"""
function pad(text, left::Int=0, right::Int=0)::String
    return " "^left * text * " "^right
end

function pad(segments::Vector{Segment}, left::Int=0, right::Int=0)::Vector{Segment}
    return map(
        (s)->Segment(" "^left * s.text * " "^right),
        segments
    )
end

# ---------------------------------------------------------------------------- #
#                                   STACKING                                   #
# ---------------------------------------------------------------------------- #

vstack(s1::String, s2::String) = s1 * "\n" * s2

"""
    vstack(r1::RenderablesUnion, r2::RenderablesUnion)

Vertically stack two renderables to give a new renderable.
"""
function vstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = r1 isa AbstractRenderable ? r1 : Renderable(r1)
    r2 = r2 isa AbstractRenderable ? r2 : Renderable(r2)

    # get dimensions of final renderable
    w1 = length(r1.segments) > 1 ? max([s.measure.w for s in r1.segments]...) : r1.measure.w
    w2 = length(r2.segments) > 1 ? max([s.measure.w for s in r2.segments]...) : r2.measure.w
    if w1 > w2
        s1 = r1.segments

        s2 = map(
            (s)->Segment(s.text * " ".^(w1-s.measure.w)),
            r2.segments
        )
    elseif w1 < w2
        s1 = map(
            (s)->Segment(s.text * " ".^(w2-s.measure.w)),
            r1.segments
        )
        s2 = r2.segments
    else
        s1, s2 = r1.segments, r2.segments
    end

    # create segments stack
    segments::Vector{Segment} = vcat(s1, s2)
    measure = Measure(max(w1, w2), length(segments))
    return Renderable(segments, measure)
end

"""
    vstack(renderables...)

Vertically stack a variable number of renderables
"""
function vstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = vstack(renderable, ren)
    end
    return renderable
end

"""
    hstack(r1::RenderablesUnion, r2::RenderablesUnion)

Horizontally stack two renderables to give a new renderable.
"""
function hstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = r1 isa AbstractRenderable ? r1 : Renderable(r1)
    r2 = r2 isa AbstractRenderable ? r2 : Renderable(r2)

    # get dimensions of final renderable
    h1 = r1.measure.h
    h2 = r2.measure.h
    Δh = abs(h2 - h1)

    # make sure both renderables have the same number of segments
    if h1 > h2
        s1 = r1.segments
        s2 = vcat(r2.segments, [Segment(" "^(r2.measure.w)) for i in 1:(Δh)])
    elseif h1 < h2
        s1 = vcat(r1.segments, [Segment(" "^(r1.measure.w)) for i in 1:(Δh)])
        s2 = r2.segments
    else
        s1, s2, = r1.segments, r2.segments
    end

    # combine segments
    segments = [Segment(s1.text * s2.text) for (s1, s2) in zip(s1, s2)]

    return Renderable(segments, Measure(r1.measure.w+r2.measure.w, max(r1.measure.h, r2.measure.h)))
end

"""
    hstack(renderables...)

Horizonatlly stack a variable number of renderables
"""
function hstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = hstack(renderable, ren)
    end
    return renderable
end

# --------------------------------- operators -------------------------------- #

"""
    Operators for more coincise layout of renderables
"""

Base.:/(r1::RenderablesUnion, r2::RenderablesUnion) = vstack(r1, r2)

Base.:*(r1::AbstractRenderable, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractString, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractRenderable, r2::AbstractString) = hstack(r1, r2)

# ---------------------------------------------------------------------------- #
#                                LINES & SPACER                                #
# ---------------------------------------------------------------------------- #
abstract type AbstractLayoutElement <: AbstractRenderable end

# ---------------------------------- spacer ---------------------------------- #
"""
        Spacer

A box of empty text with given width and height.
"""
mutable struct Spacer <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

function Spacer(width::Int, height::Int; char::Char = ' ')
    line = char^width
    seg = Segment(line)
    segments = repeat([seg], height)
    return Spacer(segments, Measure(seg.measure.w, height))
end

Spacer(width::Number, height::Number; char::Char = ' ') = Spacer(int(width), int(height); char=char)

# ----------------------------------- vline ---------------------------------- #
"""
    vLine

A multi-line renderable with each line made of a | to create a vertical line
"""
mutable struct vLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

"""
    vLine(height::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a `vLine` given a height and, optionally, style information.
"""
function vLine(
    height::Int; style::String = "default", box::Symbol = :ROUNDED
)
    line = apply_style("["*style*"]" * eval(box).head.left * "[/"*style*"]")
    segments = repeat([Segment(line)], height)
    return vLine(segments, Measure(1, height))
end

"""
    vLine(ren::AbstractRenderable; kwargs...)

Construct a vLine with the same height as a renderable
"""
vLine(ren::AbstractRenderable; kwargs...) = vLine(ren.measure.h; kwargs...)

"""
    vLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Create a `vLine` as tall as the `stdout` console
"""
function vLine(; style::String = "default", box::Symbol = :ROUNDED)
    return vLine(console_height(stdout); style = style, box = box)
end

"""
    hLine

A 1-line renderable made of repeated character from a Box.
"""
mutable struct hLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

"""
    hLine(width::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a styled `hLine` of given width.
"""
function hLine(
    width::Int; style::String = "default", box::Symbol = :ROUNDED
)
    char = eval(box).row.mid
    segments = [Segment(char^width, style)]
    return hLine(segments, Measure(width, 1))
end

"""
    hLine(width::Number, text::String; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Creates an hLine object with texte centered horizontally.
"""
function hLine(
    width::Number,
    text::String;
    style::String = "default",
    box::Symbol = :ROUNDED,
)
    box = eval(box)
    text = apply_style(text)
    tl, tr = get_lr_widths(textlen(text))
    lw, rw = get_lr_widths(width)

    open, close, space =  "[" * style * "]",  "[/" * style * "]", " "
    line = open * get_lrow(box, lw-tl, :top; with_left=false) *
                space * text * space * open *  get_rrow(box, rw-tr, :top; with_right=false) * close

    return hLine([Segment(line, style)], Measure(width, 1))
end

"""
    hLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout`
"""
function hLine(; style::String="default", box::Symbol = :ROUNDED)
    return hLine(console_width(stdout); style = style, box = box)
end

"""
    hLine(text::AbstractString; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout` with centered text.
"""
function hLine(
    text::AbstractString; style::String="default", box::Symbol = :ROUNDED
)
    return hLine(console_width(stdout), text; style = style, box = box)
end

"""
    hLine(ren::AbstractRenderable; kwargs)

Construct an hLine with same width as a renderable
"""
hLine(ren::AbstractRenderable; kwargs...) = hLine(ren.measure.w; kwargs...)
end


text:
include("modules/rendering.jl")
include("modules/sound.jl")
include("modules/twoSourceSound.jl")

using Plots
using Printf
using LaTeXStrings
using .Rendering
using .Sound
using .TwoSourceSound

function genFreqPlane()
	ppu = 10
	imgOpts = PlotOpts(100, 100, 2, ppu)
	sndOpts = SoundPlaneOpts(10, 10, ppu)
	d = 3
	freq = 125.0

	origin1 = (5.0 - d / 2, 0.0)
	origin2 = (5.0 + d / 2, 0.0)

	basePlane = genSoundPlane(sndOpts)
	distPlane1 = genDistPlane(sndOpts, origin1)
	distPlane2 = genDistPlane(sndOpts, origin2)
	ampPlane = genAmplitudePlane(distPlane1, freq) .+ genAmplitudePlane(distPlane2, freq)
	# falloffPlane = genFalloffPlane(distPlane1) .+ genFalloffPlane(distPlane2)
	# dbPlane = genDbPlane(ampPlane, basePlane)

	print("Exporting image...")

	GenPlaneImage(ampPlane, imgOpts)
end

function genDevPlane()
	ppu = 10
	imgOpts = PlotOpts(100, 100, 2, ppu)
	d = 3

	twoSrcSndOpts = TwoSourceOpts(63.0, 125.0, d, (5.0, 0.0))

	function sndGen(x::Float64, y::Float64)
		dev = TwoSourceSound.calcAvgAmplitudeDev(twoSrcSndOpts, (x, y))
		# return 10 * log10(abs(dev))
		return dev
	end

	GenFuncImage(sndGen, imgOpts)
end

function genDGraph(freqRange::Tuple{Float64, Float64}, xRange::Tuple{Float64, Float64}, yRange::Tuple{Float64, Float64}, maxDist::Float64 = 10.0, precision::Float64=0.1)
	x = 0:0.1:maxDist
	plotFunc(d) = TwoSourceSound.calcAvgPlaneAmplitudeDev_CUDA(d, freqRange, xRange, yRange)

	y = Array{Float64}(undef, length(x))
	donecnt = 0.0
	timeStart = round(time())
#Threads.@threads
	for i in 1:length(x)
		y[i] = plotFunc(x[i])
		donecnt += 1
		println("$(@sprintf("%.2f", donecnt / length(x) * 100))%")
	end

	println("Time taken: $(round(time()) - timeStart)s")

	plotly()

	freqStr = "[$(freqRange[1]), $(freqRange[2])]"
	xRangeStr = "[$(xRange[1]), $(xRange[2])]"
	yRangeStr = "[$(yRange[1]), $(yRange[2])]"

	title = "f ∈ $(freqStr); x ∈ $(xRangeStr); y ∈ $(yRangeStr)"
	plot(x, y, xlabel="d, m", ylabel="A", title=title)
	savefig("plot-F$(Int(freqRange[1])):$(Int(freqRange[2]))Hz-X$(Int(xRange[1])):$(Int(xRange[2]))m-Y$(Int(yRange[1])):$(Int(yRange[2]))m-D:$(Int(maxDist))m")
end

print(TwoSourceSound.calcAvgPlaneAmplitudeDev(5.0, (63.0, 127.0), (0.0, 1.0), (0.0, 1.0)))
# genDGraph((63.0, 65.0), (-3.0, 3.0), (0.0, 3.0), 5.0)
# genDevPlane()
# genFreqPlane()


text:
"""
    Fermi.CoupledCluster.RCCSDpT

Fermi struct that holds information about RCCSD(T) wavefunctions

_struct tree:_

**RCCSD(T)** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
struct RCCSDpT{T} <: AbstractCCWavefunction
    CCSD::RCCSD{T}
    correction::T
end

function select_precision(A::String)
    implemented = Dict{Any,Any}("single" => Float32,
                                "double" => Float64)
    try
        return implemented[A]
    catch KeyError
        throw(Fermi.InvalidFermiOptions("Invalid precision: $(A)"))
    end
end

#implementations
include("ijk.jl")


text:
## C E L L

EmptyCell(p::AbstractString) = EmptyCell(GlobalPosition(p))
EmptyCell(p::GlobalPosition) = ConstCell(p,missing)

# construction of a cell given user input as string
Cell(p::AbstractString, formula) = Cell(GlobalPosition(p),formula)
function Cell(p::GlobalPosition, formula::AbstractString)
    isempty(formula) && EmptyCell(p)

    if first(formula) == '='
        content = Meta.parse(chop(formula,head=1,tail=0))
        if isa(content, Expr)
            refs = extractRefsFromExpression(content)

            # sanitise spreadsheet-style references that are not compatible with Julia
            sanitised_refs = replace.(refs,[":" => "to"])
            sanitised_formula = formula
            for r in filter(contains(":"),refs)
                sanitised_r = replace(r, ":" => "to")
                sanitised_formula = replace(sanitised_formula, r => sanitised_r)
            end
            sanitised_content = Meta.parse(chop(sanitised_formula,head=1,tail=0))

            if isempty(refs)
                func = @eval () -> $sanitised_content
            else
                func = @eval $(Meta.parse(join(sanitised_refs,", "))) -> $sanitised_content
            end

            ancestors = GlobalisedRef.(refs,[p.sheet])
            return DynamicCell(p, formula, func, ancestors, missing)
        elseif isa(content,Symbol)
            func = identity
            ancestors = [GlobalisedRef(string(content),p.sheet)]
            return DynamicCell(p, formula, func, ancestors, missing)
        end
    else
        content = Meta.parse(formula)
    end

    value = isa(content,Number) ? content : formula # number or string
    ConstCell(p,value)
end

Base.show(io::IO, c::ConstCell) = Base.print(io, "Cell = " * string(c.value))
Base.show(io::IO, c::DynamicCell) = Base.print(io, "Cell = " * string(c.value) * " " * c.formula)
# Base.show(io::IO, c::CloneCell) = Base.print(io, "Cell = " * string(c.value) * " " * c.cloneOf.formula)

formula(c::DynamicCell) = c.formula
formula(c::ConstCell) = string(c.value)
formula(_::Missing) = missing

value(c::Cell) = c.value
value(_::Missing) = missing


text:
using Documenter, LatticeBoltzmann

makedocs(
    modules = [LatticeBoltzmann],
    format = :html,
    checkdocs = :exports,
    sitename = "LatticeBoltzmann.jl",
    pages = Any["index.md"],
)

deploydocs(repo = "github.com/MarkRedeman/LatticeBoltzmann.jl.git")


text:
using Random
include("structure.jl")
include("initialize.jl")

"""
shuffle_atoms(Backpack::Backpack, perms::Int)

This function creates 'perms' permutations of a given supercell.
"""
function shuffle_atoms(Backpack::Backpack, perms::Int)
    X = Backpack
    seq = [X.element_dict[key] for key in X.supercell.species]
    unique_perms = Array{Vector{Tuple{Int64,Int64}}}(undef, perms)
    unique_perms[1] = seq
    i = 1
    while i < perms
        current_sequence = sublattice_shuffle(seq, length(X.supercell.composition))
        v = 0
        for j = 1:i
            if current_sequence == unique_perms[j]
                v += 1
                break
            end
        end
        if v == 0
            i += 1
            unique_perms[i] = current_sequence
        end
    end
    return unique_perms
end
function sublattice_shuffle(seq, sublattices::Int)
    new_seq = Array{Any}(undef, length(seq))
    for i = 1:sublattices
        indices = findall(x -> x[1] == i, seq)
        new_indices = shuffle(indices)
        for j = 1:length(indices)
            new_seq[indices[j]] = seq[new_indices[j]]
        end
    end
    return new_seq
end

"""
all_permutations(x::T, prefix=T())

Generates all permutations of a given vector.
"""
function all_permutations(x::T, prefix=T()) where T
    if length(x) == 1
        return [[prefix; x]]
    else
        t = T[]
        for i in eachindex(x)
            if i > firstindex(x) && x[i] == x[i-1]
                continue
            end
            append!(t, all_permutations([x[begin:i-1];x[i+1:end]], [prefix; x[i]]))
        end
        return t
    end
end

max_stars_repo_name: ulipoto/WCSRO

max_issues_repo_name: ulipoto/WCSRO

max_forks_repo_name: ulipoto/WCSRO


text:
using LazySets
AFFINE_SYSTEM = ConstrainedLinearControlContinuousSystem

@testset "Bouncing ball" begin
    H = readsxmodel("../examples/bball/bball_flattened.xml", ST=AFFINE_SYSTEM)

    @test nmodes(H) == 1
    @test ntransitions(H) == 1

    # test variables names in the given order
    @test collect(keys(H.ext[:variables])) == [:x, :v]

    S = H.modes[1]
    @test S.A == [0.0 1.0; 0.0 0.0]
    @test S.B == reshape([0.0, -1.0], (2, 1))
    @test S.X[1] isa HalfSpace && S.X[1].a == [-1.0, 0.0] && S.X[1].b == 0.0 # -x <= 0 <=> x >= 0
    @test S.U[1] isa Singleton && S.U[1].element == [1.0]

    R = H.resetmaps[1]
    @test R.A == [1.0 0.0; 0.0 -0.75]
    @test isempty(R.B) && size(R.B) == (2, 0)
    @test length(R.X) == 2
    @test R.X[1] isa Hyperplane && R.X[1].a == [1.0, 0.0] && R.X[1].b == 0.0 # x = 0
    @test R.X[2] isa HalfSpace && R.X[2].a == [0.0, 1.0] && R.X[1].b == 0.0 # v < 0
end


text:
using PhyloNetworks
using DataFrames
using StatsModels

# "phylolm: binary tree, withinspecies var"
net = readTopology("(((t4:0.2537636499,t5:0.2537636499):0.2103870459,t3:0.4641506959):0.5358493041,(t1:0.4807642475,t2:0.4807642475):0.5192357525);")
df = DataFrame(
      y=[11.399108,3.216645,13.648011,4.851454,-4.922803],
      y_sd=[0.2463003,0.3236629,0.2458547,0.4866844,0.3434582],
      y_n=fill(5,5),
      x1=[0.7894387,0.3285286,2.1495131,-1.4935665,-2.3199935],
      x2=[2.8184268,0.4740687,2.5801004,1.9967379,-0.4121874],
      species=net.names
)

n = 5
starnet = readTopology(PhyloNetworks.startree_newick(n))

# for 0 and 1 - can compare results to ordinary BM model
# lambda = 0, compare to star tree
m5 = phylolm(@formula(y ~ x1), df, starnet; reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)
#=
StatsModels.TableRegressionModel{PhyloNetworkLinearModel, Matrix{Float64}}

Formula: y ~ 1 + x1

Model: Brownian motion

Parameter Estimates, using REML:
phylogenetic variance rate: 15.5623
within-species variance: 0.116138

Coefficients:
──────────────────────────────────────────────────────────────────────
               Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────
(Intercept)  6.03325     1.76961  3.41    0.0422   0.401575   11.6649
x1           3.61458     1.09896  3.29    0.0461   0.117202    7.11195
──────────────────────────────────────────────────────────────────────
Log Likelihood: -21.3324812448
AIC: 50.6649624896
=#

m4 = phylolm(@formula(y ~ x1), df, net; model="lambda", reml=true,
      tipnames=:species, withinspecies_var=true, fixedValue=0.0, y_mean_std=true)
#=
StatsModels.TableRegressionModel{PhyloNetworkLinearModel, Matrix{Float64}}

Formula: y ~ 1 + x1

Model: Pagel's lambda

Parameter Estimates, using REML:
phylogenetic variance rate: 15.5623
Lambda: 0
within-species variance: 0.116138

Coefficients:
──────────────────────────────────────────────────────────────────────
               Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────
(Intercept)  6.03325     1.76961  3.41    0.0422   0.401575   11.6649
x1           3.61458     1.09896  3.29    0.0461   0.117202    7.11195
──────────────────────────────────────────────────────────────────────
Log Likelihood: -21.3324812448
AIC: 52.6649624896
=#

# ordinary BM
m3 = phylolm(@formula(y ~ x2), df, net; reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)
#=
StatsModels.TableRegressionModel{PhyloNetworkLinearModel, Matrix{Float64}}

Formula: y ~ 1 + x2

Model: Brownian motion

Parameter Estimates, using REML:
phylogenetic variance rate: 9.30581
within-species variance: 0.116126

Coefficients:
─────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept)  -0.803512    2.11612   -0.38    0.7294   -7.53794    5.93092
x2            4.1111      0.702594   5.85    0.0100    1.87513    6.34707
─────────────────────────────────────────────────────────────────────────
Log Likelihood: -19.6834453226
AIC: 47.3668906453
=#

m2 = phylolm(@formula(y ~ x2), df, net; model="lambda", reml=true,
      tipnames=:species, withinspecies_var=true, fixedValue=1.0, y_mean_std=true)
#=
StatsModels.TableRegressionModel{PhyloNetworkLinearModel, Matrix{Float64}}

Formula: y ~ 1 + x2

Model: Pagel's lambda

Parameter Estimates, using REML:
phylogenetic variance rate: 9.30581
Lambda: 1
within-species variance: 0.116126

Coefficients:
─────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept)  -0.803512    2.11612   -0.38    0.7294   -7.53794    5.93092
x2            4.1111      0.702594   5.85    0.0100    1.87513    6.34707
─────────────────────────────────────────────────────────────────────────
Log Likelihood: -19.6834453226
AIC: 49.3668906453
=#

# lambda = 0.5
m1 = phylolm(@formula(y ~ x1), df, net; model="lambda", reml=true,
      tipnames=:species, withinspecies_var=true, fixedValue=0.5, y_mean_std=true)
#=
StatsModels.TableRegressionModel{PhyloNetworkLinearModel, Matrix{Float64}}

Formula: y ~ 1 + x1

Model: Pagel's lambda

Parameter Estimates, using REML:
phylogenetic variance rate: 21.8812
Lambda: 0.5
within-species variance: 0.11614

Coefficients:
──────────────────────────────────────────────────────────────────────
               Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────
(Intercept)  6.2069      2.5429   2.44    0.0924   -1.88573   14.2995
x1           3.93469     1.42465  2.76    0.0700   -0.59919    8.46858
──────────────────────────────────────────────────────────────────────
Log Likelihood: -21.4050221575
AIC: 52.810044315
=#

# optimize lambda case (where fixedValue is not provided)
m6 = phylolm(@formula(y ~ x1 + x2), df, net; model="lambda", reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)

m7 = phylolm(@formula(y ~ x1 + x2), df, net; reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)


text:
type ExternalResolveBurgerSS <: ExternalAdectionResolver
  n_off :: Int64
  x :: BurgerData
end

type ExternalResolveBurgerLS <: ExternalAdectionResolver
  n_off :: Int64
  x :: BurgerData
end

type ExternalResolveBurgerSL <: ExternalAdectionResolver
  n_off :: Int64
  x :: BurgerData
end

type ExternalResolveBurgerLL <: ExternalAdectionResolver
  n_off :: Int64
  x :: BurgerData
end

export ExternalResolveBurgerLL, ExternalResolveBurgerLS, ExternalResolveBurgerSL, ExternalResolveBurgerSS

include("advectionResolves.jl")
include("advectionExternals.jl")
export advection


text: # File For Back Propagation Algorithm


text:
println("loading packages... ")
# Parallele per tutto il bbin
using ArgParse
using OrthographicProjection

println("packages OK")

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table! s begin
	"source"
		help = "A text file with Potree directory list"
		arg_type = String
		required = true
	"--projectname", "-p"
		help = "Project name"
		arg_type = String
		required = true
	"--output", "-o"
		help = "Output folder"
		arg_type = String
		required = true
	"--bbin"
        help = "Bounding box as 'x_min y_min z_min x_max y_max z_max' or Potree JSON volume model"
		required = true
	"--step"
		help = "Distance between sections"
		arg_type = Float64
		default = 0.
	"--p1"
		help = "Start point"
		arg_type = String
		required = true
	"--p2"
		help = "End point"
		arg_type = String
		required = true
	"--axis"
		help = "A vector in plane"
		arg_type = String
		default = "0 0 1"
	"--thickness"
		help = "Section thickness"
		arg_type = Float64
		default = 0.1
	end

	return parse_args(s)
end

function main()
	args = parse_commandline()

	OrthographicProjection.flushprintln("== params ==")
	for (arg,val) in args
		OrthographicProjection.flushprintln("$arg  =>  $val")
	end

	txtpotreedirs = args["source"]
	project_name = args["projectname"]
	output_folder = args["output"]
	bbin = args["bbin"]
	step = args["step"]
	p1_ = args["p1"]
	p2_ = args["p2"]
	axis_ = args["axis"]
	thickness = args["thickness"]

	# p = tryparse.(Float64,split(plane, " "))
	# @assert length(p) == 4 "$plane: Please described the plane in Hessian normal form"
	# plane = OrthographicProjection.Plane(p[1],p[2],p[3],p[4])

	b = tryparse.(Float64,split(bbin, " "))
	if length(b) == 6
		bbin = OrthographicProjection.AABB(b[4],b[1],b[5],b[2],b[6],b[3])
	end

	p1 = tryparse.(Float64,split(p1_, " "))
	@assert length(p1) == 3 "a 3D point needed"
	p2 = tryparse.(Float64,split(p2_, " "))
	@assert length(p2) == 3 "a 3D point needed"
	axis_y = tryparse.(Float64,split(axis_, " "))
	@assert length(axis_y) == 3 "a 3D axis needed"

	try
		proj_folder, plane, model = OrthographicProjection.preprocess(project_name, output_folder, bbin, p1, p2, axis_y, thickness)
		OrthographicProjection.get_parallel_sections(txtpotreedirs, project_name, proj_folder, bbin, step, plane, model)
	catch y

	end

end

@time main()


text:
using Printf
import Flux
import GridGames

OVERRIDE_MSG = "This should be overridden in concrete implemention of AbstractEnv"

get_marker(::AbstractEnv, ::Integer, ::Integer)   = @assert false OVERRIDE_MSG
get_coord(::AbstractEnv, action::Integer)         = @assert false OVERRIDE_MSG
getstates(::Type{<:AbstractEnv})                  = @assert false OVERRIDE_MSG
getstatesplus(::Type{<:AbstractEnv})              = @assert false OVERRIDE_MSG

makepolicytests(::Type{<:AbstractEnv}) = nothing
legal_actions(g::AbstractEnv) = legal_actions(get_state(g))
svg(s::AbstractEnv) = BoardGameSVG(s)
get_board_state(x::AbstractEnv) = get_state(x).board_state
get_state(x::AbstractEnv) = x.states[end]
get_state(x::AbstractEnv, idx::Integer) = x.states[idx]
num_rows(T::Type{<:AbstractEnv}) = num_rows(getstatetype(T))
num_rows(g::AbstractEnv) = num_rows(typeof(g))
num_cols(T::Type{<:AbstractEnv}) = num_cols(getstatetype(T))
num_cols(g::AbstractEnv) = num_cols(typeof(g))
num_cells(T::Type{<:AbstractEnv}) = num_cells(getstatetype(T))
num_cells(g::AbstractEnv) = num_cells(typeof(g))
image_scale(x::AbstractEnv) = image_scale(getstatetype(x))
image_title_offset(x::AbstractEnv) = image_title_offset(getstatetype(x))
player1_markers(x::AbstractEnv) = player1_markers(getstatetype(x))
player2_markers(x::AbstractEnv) = player2_markers(getstatetype(x))
getstatetype(x::AbstractEnv) = GridGames.getstatetype(typeof(x))
to_play(x::AbstractEnv) = to_play(x, length(x.history) + 1)
to_play(x::AbstractEnv, state_idx::Integer) = isodd(state_idx) ? PLAYER_1 : PLAYER_2
get_pieces(g::AbstractEnv) = get_pieces(getstatetype(g))
get_pieces(T::Type{<:AbstractEnv}) = get_pieces(getstatetype(T))
action_space_size(g::AbstractEnv) = action_space_size(getstatetype(g))
action_space_size(T::Type{<:AbstractEnv}) = action_space_size(GridGames.getstatetype(T))
Base.:-(g1::AbstractEnv, g2::AbstractEnv) = get_state(g1) - get_state(g2)
piece(g::AbstractEnv, name::Symbol) = piece(getstatetype(s), name)
piece(T::Type{<:AbstractEnv}, name::Symbol) = piece(getstatetype(T), name)
init_state(g::AbstractEnv, pieces::Vector{Piece}) = init_state(getstatetype(g), pieces)
init_state(T::Type{<:AbstractEnv}, pieces::Vector{Piece}) = init_state(getstatetype(T), pieces)

function MuZero.policypredsvg(g::AbstractEnv, test_result::PolicyTestResult)
    test = test_result.test
    policy = test_result.prediction
    g = test.game
    sorted_actions = sortperm(policy)
    top3_actions = reverse(sorted_actions)[1:3]
    top3 = map(x->begin
        pr_idx, action = x
        pred_pr = @sprintf "%.2f" policy[action]
        g_copy = copy(g)
        is_max = in(action, test.actions)
        step!(g_copy, action)
        highlight_coord = diff_coord(g_copy, g)
        highlight_color = is_max ? "#AADBC0" : "#ECABA7"
        gridsvg(
            get_state(g_copy);
            title="prediction #$pr_idx",
            title_offset=(-210,-195),
            tr_header="pr=$(pred_pr)",
            tr_offset=(-310,-295),
            highlight=highlight_coord,
            highlight_color=highlight_color,
        )
    end, enumerate(top3_actions))

    GridSVG(top3)
end

function diff_coord(g1::AbstractEnv, g2::AbstractEnv)
    diff = (g1 - g2)
    isempty(diff) ? nothing : diff[1]
end

Base.copy(g::AbstractEnv) = copy_at(g, length(g.states))

function copy_at(g::AbstractEnv, state_idx::Integer)
    typeof(g)(
        copy(gethistory(g)[1:(state_idx-1)]),
        copy(g.rewards[1:(state_idx-1)]),
        copy(g.states[1:state_idx]),
    )
end

# @todo change this to get_image(::AbstractEnv, ::Action)
get_action_image(g::AbstractEnv, a::Integer) = get_action_image(GridGames.getstatetype(g), a)
get_action_image(T::Type{<:AbstractEnv}, a::Integer) = get_action_image(GridGames.getstatetype(T), a)
get_action_image(T::Type{<:AbstractEnv}, a::AbstractEnvAction) = get_action_image(T, a.idx)

function get_image(g::AbstractEnv, state_idx=length(g.states))
    @assert 0 < state_idx <= length(g.states) "state_idx is $state_idx, but the number of states is $(length(g.states))"
    get_image(get_state(g, state_idx))
end

get_reward(g::AbstractEnv, action::Integer) = get_reward(get_state(g), action)

function init_game(T::Type{<:AbstractEnv}, actions::Vector{<:Integer})
    g = T()
    foreach(a->step!(g,a), actions)
    g
end

function coord_to_byte_idx(T::Type{<:AbstractEnv}, coord::Tuple{<:Integer, <:Integer})
    coord_to_byte_idx(getstatetype(T), coord)
end

function get_piece(g::AbstractEnv, name::Symbol)
    get_pieces(g)[get_piece_idx(g, name)]
end

get_piece_idx(g::AbstractEnv, name::Symbol) = get_piece_idx(getstatetype(g), name)

# @todo remove all dependencies on "apply!" (and use step! instead)
apply!(g::AbstractEnv, a::Integer) = step!(g, a)
isterminal(g::AbstractEnv) = length(getlegalactions(g)) == 0
islegal(g::AbstractEnv, action::Integer) = in(action, legal_actions(g))
game_is_won(g::AbstractEnv) = game_is_won(get_state(g))
get_winner(g::AbstractEnv) = get_state(g) |> get_winner
#next_to_play(g::AbstractEnv) = to_play(g) == PLAYER_1 ? PLAYER_2 : PLAYER_1
discount(::AbstractEnv) = Float32(1)

# @todo eventually remove these
is_terminal(g::AbstractEnv) = isterminal(g)
current_image(g::AbstractEnv) = get_image(g)
make_image(g::AbstractEnv, i::Integer) = get_image(g, i)
to_play_at(g::AbstractEnv, i::Integer) = to_play(g, i)
function make_image_label(g::AbstractEnv, image_idx::Integer)
    if 0 < image_idx <= length(g.states)
        actions = gethistory(g)[1:(image_idx-1)]
        return join(actions, ",")
        players = player_history(g)[1:(image_idx-1)]
        return join(map(pair->join(pair), zip(players, actions)), ",")
    end
    @assert false "image_idx is $image_idx, but the number of states is $(length(g.states))"
end
function player_history(g::AbstractEnv)
    players = ["A","B"]
    map(idx->players[(idx-1)%2+1], 1:length(g.history))
end
function reward_matrix(g::AbstractEnv)
    rewards_with_initial_zero = [0., g.rewards...]
    m = zeros(Float32, length(rewards_with_initial_zero), 2)
    for (idx,reward) in enumerate(rewards_with_initial_zero)
        m[idx, 1] = iseven(idx) ? reward : -1*reward
        m[idx, 2] = iseven(idx) ? -1*reward : reward
    end
    m
end
function value_matrix(g::AbstractEnv, from::Integer=1; td_steps)
    to = min(length(gethistory(g)) + 1, from + td_steps)
    value_matrix_support(g, from, to)
end
function value_matrix_support(g::AbstractEnv, from::Integer=1, to=(length(g.history)+1))
    reward_from = from + 1
    reward_to = to
    rewards = reward_matrix(g)[reward_from:reward_to, :]
    function cumvalue(next_value, row)
        i, reward = row
        reward * discount(g)^i + next_value
    end
    mat = accumulate(
        cumvalue,
        enumerate(eachrow(reverse(rewards; dims=1)));
        init=[0., 0.]) |> reverse
    [mat..., [0.,0.]]
end

# function step!(g::AbstractEnv, actions::Vector{<:Integer})
#     map(a->step!(g, a), actions)
# end

# function step!(g::AbstractEnv, action::Integer)
#     if !isterminal(g) && islegal(g, action)
#         push!(g.rewards, get_reward(g, action))
#         push!(g.states, next_state(g, action))
#         push!(g.history, action)
#         return g.rewards[end]
#     end
#     nothing
# end

next_state(g::AbstractEnv, action::Integer) = next_state(get_state(g), action)

function Base.getindex(g::AbstractEnv, row_idx::Integer, col_idx::Integer)
    get_state(g)[row_idx, col_idx]
end

function Base.getindex(g::AbstractEnv, ::Colon, col_idx::Integer)
    get_state(g)[:, col_idx]
end

function Base.getindex(g::AbstractEnv, row_idx::Integer, ::Colon)
    get_state(g)[row_idx, :]
end

function Base.getindex(g::AbstractEnv, row_idx::Integer, col_range::UnitRange{<:Integer})
    get_state(g)[row_idx, col_range]
end

function Base.getindex(g::AbstractEnv, row_range::UnitRange{<:Integer}, col_idx::Integer)
    get_state(g)[row_range, col_idx]
end

function Base.show(io::IO, g::AbstractEnv)
    game_type_name = string(typeof(g))
    state_type_name = string(getstatetype(g))
    print(io, "$game_type_name{$state_type_name{\"$(encode(get_state(g)))\"}}")
end

# function Base.show(io::IO, m::MIME"text/plain", g::AbstractEnv)
#     show(io, "text/plain", get_state(g))
#     println(io)
#     if game_is_won(g)
#         println(io, "Winner: Player $(get_winner(g))")
#     else
#         println(io, "To Play: $(to_play(g))")
#     end
#     println(io, "History: $(join(g.history, ", "))")
# end

byte_idx_to_coord(g::AbstractEnv, byte_idx::Integer) = byte_idx_to_coord(getstatetype(g), byte_idx)
byte_idx_to_coord(T::Type{<:AbstractEnv}, byte_idx::Integer) = byte_idx_to_coord(getstatetype(T), byte_idx)

# function gridsvg(g::AbstractEnv; main_title="image", next_action=nothing, header=nothing, tr_header=nothing)
#     function cell_fn(_,_,row,col)
#         piece_idx = s[row,col]
#         if piece_idx == 0
#             return nothing
#         end
#         m = all_markers(s)[piece_idx]
#         (value=m.icon_symbol, color=m.icon_color, size=m.icon_size)
#     end
#     GridSVG(
#         Grid(num_cols(s), num_rows(s), cell_fn),
#         main_title=main_title,
#         main_title_offset=image_title_offset(s),
#         tr_header=tr_header,
#         header=header,
#         aspect_ratio=(num_cols(s), num_rows(s)),
#         scale=image_scale(s))
# end

function Base.show(io::IO, m::MIME"image/svg+xml", g::AbstractBoardGameSVG)
    show(io, m, gridsvg(g.x))
end

actions(g::AbstractEnv) = actions(getstatetype(g))
actions(T::Type{<:AbstractEnv}) = actions(getstatetype(T))
uniformpolicy(T::Type{<:AbstractEnv}) = (s)->rand(actions(T))

function get_boardgame_config(game_type::Type{<:AbstractEnv};
    network_type, action_space_size, max_moves, root_dirichlet_alpha,
    lr_init, support_size, observation_shape,
    stacked_observations, hidden_state_size, checkpoint_interval,
    num_actors, weight_decay, save_interval, num_unroll_steps)

    function visit_softmax_temperature(num_moves, num_steps)
        if num_moves < 30
            return 1.0
        end
        0.0  # Play according to the max.
    end

    tictactoe_visit_softmax_temperature(num_moves, num_steps) = 1.0

    return MuZeroConfig(
        game_type,
        network_type=network_type,
        action_space_size=action_space_size,
        max_moves=max_moves,
        discount=1.0,
        root_dirichlet_alpha=root_dirichlet_alpha,
        num_simulations=100,
        training_steps=2000000,
        batch_size=64,
        td_steps=max_moves, # Always use Monte Carlo return.
        num_actors=num_actors,
        lr_init=lr_init,
        lr_decay_steps=400e3,
        visit_softmax_temperature_fn=tictactoe_visit_softmax_temperature,
        known_bounds=KnownBounds(-1, 1),
        support_size=support_size,
        observation_shape=observation_shape,
        stacked_observations=stacked_observations,
        hidden_state_size=hidden_state_size,
        checkpoint_interval=checkpoint_interval,
        weight_decay=weight_decay,
        save_interval=save_interval,
        num_unroll_steps=num_unroll_steps,
    )
end

function getconfig(g::Type{<:AbstractEnv})
    get_boardgame_config(
        g,
        network_type=Incremental,
        max_moves=9,
        root_dirichlet_alpha=0.3,
        action_space_size=9, # Number of all possible actions in the game
        support_size = 10, # Value and reward are scaled (with almost sqrt) and encoded on a vector with a range of -support_size to support_size. Choose it so that support_size <= sqrt(max(abs(discounted reward)))
        observation_shape = (3, 3, 3),
        stacked_observations = 0, # Number of previous observations and previous actions to add to the current observation
        hidden_state_size = 32, # The 1-dimensional size of the encoded game state
        checkpoint_interval = 100, #  Number of training steps before using the model for self-playing
        num_actors = 3,
        weight_decay = 1e-3,
        lr_init=0.002f0,
        save_interval = 10000,
        num_unroll_steps=1,
    )
end

# function Base.show(io::IO, m::MIME"image/svg+xml", g::AbstractEnvSVG)
#     game = g.x
#     svgs = map(action_idx->gridsvg(game, action_idx), 2:length(game.history))
#     fig = GridSVG(svgs, :, 2;
#         main_title="Game",
#         main_title_offset=(70,40),
#         main_title_size=40)
#     show(io, m, fig)
# end

function gridsvg(g::AbstractEnv)
    title_offset = (-210, -200)
    title_size = 30

    svgs = map(x->begin
        (action_idx, state) = x
        curr = state
        prev = g.states[action_idx]
        highlight_coord = diff_coord(curr, prev)
        gridsvg(
            state;
            highlight=highlight_coord,
            title="action $action_idx",
            title_offset=title_offset,
            title_size=title_size)
    end, enumerate(g.states[2:end]))

    if isempty(svgs)
        push!(svgs,
            gridsvg(
                g.states[1];
                title="initial state",
                title_offset=title_offset,
                title_size=title_size
            ))
    end

    fig = GridSVG(svgs, :, 2;
        main_title="Game",
        main_title_offset=(70,40),
        main_title_size=40)

    # state = get_state(g, action_idx)
    # player = to_play(g, action_idx)
    # next_action = g.history[action_idx]
    # gridsvg(state;
    #     title="action $action_idx",
    #     title_offset=(-210, -200),
    #     title_size=30)
end


text:
export primaldualsolver

mutable struct PrimalDualSolver{T,S} <: AbstractLinearSolver
    S::Matrix{T}
    reg::Vector{Regularization}
    regName::Vector{String}
    gradientOp::S
    u::Vector{T}
    c::Vector{T}
    cO::Vector{T}
    y1::Vector{T}
    y2::Vector{T}
    σ::T
    τ::T
    ϵ::T
    PrimalDualGap::T
    enforceReal::Bool
    enforcePositive::Bool
    iterations::Int64
    shape::NTuple{2,Int64}
end

function PrimalDualSolver(S::Matrix{T}; b=nothing, λ=1e-4, reg = nothing
			 , regName = "L1"
                         , gradientOp = nothing
                         , enforceReal::Bool=false
                         , enforcePositive::Bool=false
                         , iterations::Int64=10
                         , σ=1.0
                         , τ=1.0
                         , ϵ=1e-10
                         , PrimalDualGap=1.0
                         , shape::NTuple{2,Int64}=(0,0)
                         , kargs...) where T

  if reg == nothing
    if typeof(λ)<:AbstractFloat && typeof(regName)==String
        reg = [Regularization(regName, T(λ); kargs...)]
    elseif typeof(λ)<:AbstractVector && typeof(regName)<:AbstractVector{String}
        reg = Regularization(regName, T(λ); kargs...)
    else
        error("could not initialite regularization")
    end
  end

  M,N = size(S)

  if shape == (0,0)
    shape = (Int(sqrt(N)),Int(sqrt(N)))
  else
    shape = shape
  end

  if typeof(regName)==String
      regName = [regName]
  end

  if regName[1] == "L1"
    gradientOp = opEye(T,N) #UniformScaling(one(T))
  elseif regName[1] == "TV"
    gradientOp = gradientOperator(T,shape)
  end

  if b != nothing
    u = b
  else
    u = zeros(T,M)
  end

  c = zeros(T,N)
  cO = zeros(T,N)
  y1 = zeros(T,M)
  y2 = zeros(T,size(gradientOp*c,1))

  return PrimalDualSolver(S,reg,regName,gradientOp,u,c,cO,y1,y2,T(σ),T(τ),T(ϵ),T(PrimalDualGap),enforceReal,enforcePositive,iterations,shape)
end

function init!(solver::PrimalDualSolver; S::Matrix{T}=solver.S, u::Vector{T}=T[], c::Vector{T}=T[],
    PrimalDualGap::T=solver.PrimalDualGap) where {T,R}

  solver.u[:] .= u
  solver.PrimalDualGap = (1/2)*(norm(solver.u,2))^2

  # start vector
  if isempty(c)
    solver.c[:] .= zero(T)
  else
    solver.c[:] .= c
  end
  solver.cO[:] .= zero(T)
  solver.y1[:] .= zero(T)
  solver.y2[:] .= zero(T)

end

function solve(solver::PrimalDualSolver, u::Vector{T}; S::Matrix{T}=solver.S, startVector::Vector{T}=eltype(S)[]
              , solverInfo=nothing, PrimalDualGap::T=solver.PrimalDualGap, kargs...) where {T}

  # initialize solver parameters
  init!(solver; S=S, u=u, c=startVector, PrimalDualGap=solver.PrimalDualGap)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.c,solver.cO,solver.y1,solver.y2)

  # perform the iterations
  for item in solver
    solverInfo != nothing && storeInfo(solverInfo,solver.c,solver.cO,solver.y1,solver.y2)
  end

  return solver.c
end

function iterate(solver::PrimalDualSolver, iteration::Int=0)
  if done(solver,iteration) return nothing end

  # updating dual variables
  for i=1:length(solver.reg)
      solver.y1 .= (solver.y1 + solver.σ*(solver.S*solver.c - solver.u))./(1+solver.σ)
      if solver.regName[1] == "L1"
          solver.y2 .= ProxL1Conj(solver.y2 + solver.σ*solver.gradientOp*solver.c, solver.reg[1].λ, solver.shape)
      elseif solver.regName[1] == "TV"
          solver.y2 .= ProxTVConj(solver.y2 + solver.σ*solver.gradientOp*solver.c, solver.reg[1].λ, solver.shape)
      end
  end

  # updating primal variable
  for i=1:length(solver.reg)
      solver.c += - solver.τ*(adjoint(solver.S)*solver.y1 + adjoint(solver.gradientOp)*solver.y2)
  end

  applyConstraints(solver.c, nothing, solver.enforceReal, solver.enforcePositive)

  # updating convergence measure
  for i=1:length(solver.reg)
    solver.PrimalDualGap = abs((1/2)*norm(solver.S*solver.c-solver.u)^2 + solver.reg[1].λ*norm(solver.c,1) + (1/2)*norm(solver.y1,2)^2 + dot(solver.y1,solver.u))
  end

  return solver.y1, iteration+1

end

function converged(solver::PrimalDualSolver)
  if (solver.PrimalDualGap) >= solver.ϵ
    return false
  end

  return true
end

@inline done(solver::PrimalDualSolver,iteration::Int) = converged(solver) || iteration>=solver.iterations

# Proximal map of the convex conjugate of the l1 norm
function ProxL1Conj(x::Vector{T},α::T,shape::NTuple{2,Int64}) where T
   m,n = shape
   p1 = reshape(x,m,n)

   # threshold p1
   for j=1:n, i=1:m
       p1[i,j] = sign(p1[i,j])*min(norm(p1[i,j]),α)
   end
   for i=1:m
       p1[i,n] = sign(p1[i,n])*min(abs(p1[i,n]),α)
   end

   return vcat(vec(p1))
end

# Gradient
function BB(u::Vector{T}, shape::NTuple{2,Int64}) where T
 m,n = shape
 u = reshape(u,m,n)

 p1 = zeros(T,m,n);  p2 = zeros(T,m,n);  p3 = zeros(T,m,n);  p4 = zeros(T,m,n)
 p5 = zeros(T,m,n);  p6 = zeros(T,m,n);  p7 = zeros(T,m,n);  p8 = zeros(T,m,n)

    p1[1:m-1,1:n] = u[1:m-1,1:n]-u[2:m,1:n]
    p2[1:m,1:n-1] = u[1:m,1:n-1]-u[1:m,2:n]
    p3[1:m-1,1:n-1] = u[1:m-1,1:n-1]-u[2:m,2:n]
    p4[1:m-1,2:n] = u[1:m-1,2:n]-u[2:m,1:n-1]
    p5[1:m-2,1:n-1] = u[1:m-2,1:n-1]-u[3:m,2:n]
    p6[1:m-2,2:n] = u[1:m-2,2:n]-u[3:m,1:n-1]
    p7[1:m-1,1:n-2] = u[1:m-1,1:n-2]-u[2:m,3:n]
    p8[1:m-1,3:n] = u[1:m-1,3:n]-u[2:m,1:n-2]

 return vcat(vec(p1),vec(p2),vec(p3),vec(p4),vec(p5),vec(p6),vec(p7),vec(p8))
end

# Divergence
function BBS(u::Array{T},shape::NTuple{2,Int64}) where T
   m,n = shape
   p1 = reshape(u[1:m*n],m,n);            p2 = reshape(u[m*n+1:2*n*m],m,n)
   p3 = reshape(u[2*m*n+1:3*m*n],m,n);    p4 = reshape(u[3*m*n+1:4*m*n],m,n)
   p5 = reshape(u[4*m*n+1:5*m*n],m,n);    p6 = reshape(u[5*m*n+1:6*m*n],m,n)
   p7 = reshape(u[6*m*n+1:7*m*n],m,n);    p8 = reshape(u[7*m*n+1:end],m,n)

   x = zeros(T, m, n)

   # points in the corners
   x[1,1] = p1[1,1]+p2[1,1]+p3[1,1]+p4[1,1]+p5[1,1]+p6[1,1]+p7[1,1]+p8[1,1]
   x[m,1] = p2[m,1]-p1[m-1,1]+p3[m,1]+p4[m,1]-p4[m-1,2]+p5[m,1]+p6[m,1]-p6[m-2,2]+p7[m,1]+p8[m,1]-p8[m-1,3]
   x[1,n] = p1[1,n]-p2[1,n-1]+p3[1,n]+p4[1,n]+p5[1,n]+p6[1,n]+p7[1,n]+p8[1,n]
   x[m,n] = -p1[m-1,n]-p2[m,n-1]-p3[m-1,n-1]-p5[m-2,n-1]-p7[m-1,n-2]

   # remaining points with i=1 (first row)
   x[1,2:n-1] = p1[1,2:n-1]+p2[1,2:n-1]-p2[1,1:n-2]+p3[1,2:n-1]+p4[1,2:n-1]+p5[1,2:n-1]+p6[1,2:n-1]+p7[1,2:n-1]+p8[1,2:n-1]
   # remaining points with i=2 (second row)
   x[2,3:n-2] = p1[2,3:n-2]-p1[1,3:n-2]+p2[2,3:n-2]-p2[2,2:n-3]+p3[2,3:n-2]-p3[1,2:n-3]+p4[2,3:n-2]-p4[1,4:n-1]+p5[2,3:n-2]+p6[2,3:n-2]+p7[2,3:n-2]-p7[1,1:n-4]+p8[2,3:n-2]-p8[1,5:n]
   #remaining points with i=m (last row)
   x[m,3:n-2] = -p1[m-1,3:n-2]+p2[m,3:n-2]-p2[m,2:n-3]-p3[m-1,2:n-3]+p4[m,3:n-2]-p4[m-1,4:n-1]+p5[m,3:n-2]-p5[m-2,2:n-3]+p6[m,3:n-2]-p6[m-2,4:n-1]+p7[m,3:n-2]-p7[m-1,1:n-4]+p8[m,3:n-2]-p8[m-1,5:n]
   # remaingin points with j=1 (first column)
   x[3:m,1] = p1[3:m,1]-p1[2:m-1,1]+p2[3:m,1]+p3[3:m,1]+p4[3:m,1]-p4[2:m-1,2]+p5[3:m,1]+p6[3:m,1]-p6[1:m-2,2]+p7[3:m,1]+p8[3:m,1]-p8[2:m-1,3]
   # remaingin points with j=2 (second column)
   x[3:m,2] = p1[3:m,2]-p1[2:m-1,2]+p2[3:m,2]-p2[3:m,1]+p3[3:m,2]-p3[2:m-1,1]+p4[3:m,2]-p4[2:m-1,3]+p5[3:m,2]-p5[1:m-2,1]+p6[3:m,2]-p6[1:m-2,3]+p7[3:m,2]+p8[3:m,2]-p8[2:m-1,4]
   #remaingin points with j=n (last column)
   x[3:m,n] = p1[3:m,n]-p1[2:m-1,n]-p2[3:m,n-1]-p3[2:m-1,n-1]+p4[3:m,n]+p5[3:m,n]-p5[1:m-2,n-1]+p6[3:m,n]+p7[3:m,n]-p7[2:m-1,n-2]+p8[3:m,n]
   #remaingin points with j=n-1
   x[3:m,n-1] = p1[3:m,n-1]-p1[2:m-1,n-1]+p2[3:m,n-1]-p2[3:m,n-2]+p3[3:m,n-1]-p3[2:m-1,n-2]+p4[3:m,n-1]-p4[2:m-1,n]+p5[3:m,n-1]-p5[1:m-2,n-2]+p6[3:m,n-1]-p6[1:m-2,n]+p7[3:m,n-1]-p7[2:m-1,n-3]+p8[3:m,n-1]
   #remaingin points with j=n-2
   x[3:m,n-2] = p1[3:m,n-2]-p1[2:m-1,n-2]+p2[3:m,n-2]-p2[3:m,n-3]+p3[3:m,n-2]-p3[2:m-1,n-3]+p4[3:m,n-2]-p4[2:m-1,n-1]+p5[3:m,n-2]-p5[1:m-2,n-3]+p6[3:m,n-2]-p6[1:m-2,n-1]+p7[3:m,n-2]-p7[2:m-1,n-4]+p8[3:m,n-2]-p8[2:m-1,n]

   #all remaining (inner) matrix elements
   x[3:m-1,3:n-3] = p1[3:m-1,3:n-3]-p1[2:m-2,3:n-3]+p2[3:m-1,3:n-3]-p2[3:m-1,2:n-4]+p3[3:m-1,3:n-3]-p3[2:m-2,2:n-4]+p4[3:m-1,3:n-3]-p4[2:m-2,4:n-2]+p5[3:m-1,3:n-3]-p5[1:m-3,2:n-4]+p6[3:m-1,3:n-3]-p6[1:m-3,4:n-2]+p7[3:m-1,3:n-3]-p7[2:m-2,1:n-5]+p8[3:m-1,3:n-3]-p8[2:m-2,5:n-1]

   return vec(x)
end

function gradientOperator(::Type{T},shape::NTuple{2,Int64}) where T
   M,N = shape
   ncol = M*N
   nrow = 8*M*N
   return LinearOperator{T}(nrow,ncol,false,false,x->BB(x,shape),nothing,x->BBS(x,shape))
end

# Proximal map of the convex conjugate of the debiasing function
function ProxTVConj(x::Vector{T},α::T,shape::NTuple{2,Int64}) where T
   m,n = shape
   p1 = reshape(x[1:m*n],m,n)
   p2 = reshape(x[m*n+1:2*m*n],m,n)
   p3 = reshape(x[2*m*n+1:3*m*n],m,n)
   p4 = reshape(x[3*m*n+1:4*m*n],m,n)
   p5 = reshape(x[4*m*n+1:5*m*n],m,n)
   p6 = reshape(x[5*m*n+1:6*m*n],m,n)
   p7 = reshape(x[6*m*n+1:7*m*n],m,n)
   p8 = reshape(x[7*m*n+1:end],m,n)

   # threshold p1
   for j=1:n, i=1:m
       p1[i,j] = (sqrt(5)-2)*sign(p1[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p1[i,n] = (sqrt(5)-2)*sign(p1[i,n])*min(abs(p1[i,n]),α)
   end

   # threshold p2
   for j=1:n, i=1:m
       p2[i,j] = (sqrt(5)-2)*sign(p2[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for j=1:n
       p2[m,j] = (sqrt(5)-2)*sign(p2[m,j])*min(abs(p2[m,j]),α)
   end

   # threshold p3
   for j=1:n, i=1:m
       p3[i,j] = (sqrt(5)-3*sqrt(2)/2)*sign(p3[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p3[i,n] = (sqrt(5)-3*sqrt(2)/2)*sign(p3[i,n])*min(abs(p3[i,n]),α)
   end
   for j=1:n
       p3[m,j] = (sqrt(5)-3*sqrt(2)/2)*sign(p3[m,j])*min(abs(p3[m,j]),α)
   end

   # threshold p4
   for j=1:n, i=1:m
       p4[i,j] = (sqrt(5)-3*sqrt(2)/2)*sign(p4[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p4[i,n] = (sqrt(5)-3*sqrt(2)/2)*sign(p4[i,n])*min(abs(p4[i,n]),α)
   end
   for j=1:n
       p4[m,j] = (sqrt(5)-3*sqrt(2)/2)*sign(p4[m,j])*min(abs(p4[m,j]),α)
   end

   # threshold p5
   for j=1:n, i=1:m
       p5[i,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p5[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p5[i,n] = ((1+sqrt(2)-sqrt(5))/2)*sign(p5[i,n])*min(abs(p5[i,n]),α)
   end
   for j=1:n
       p5[m,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p5[m,j])*min(abs(p5[m,j]),α)
   end

   # threshold p6
   for j=1:n, i=1:m
       p6[i,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p6[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p6[i,n] = ((1+sqrt(2)-sqrt(5))/2)*sign(p6[i,n])*min(abs(p6[i,n]),α)
   end
   for j=1:n
       p6[m,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p6[m,j])*min(abs(p6[m,j]),α)
   end

   # threshold p7
   for j=1:n, i=1:m
       p7[i,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p7[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p7[i,n] = ((1+sqrt(2)-sqrt(5))/2)*sign(p7[i,n])*min(abs(p7[i,n]),α)
   end
   for j=1:n
       p7[m,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p7[m,j])*min(abs(p7[m,j]),α)
   end

   # threshold p8
   for j=1:n, i=1:m
       p8[i,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p8[i,j])*min(norm([p1[i,j],p2[i,j],p3[i,j],p4[i,j],p5[i,j],p6[i,j],p7[i,j],p8[i,j]]),α)
   end
   for i=1:m
       p8[i,n] = ((1+sqrt(2)-sqrt(5))/2)*sign(p8[i,n])*min(abs(p8[i,n]),α)
   end
   for j=1:n
       p8[m,j] = ((1+sqrt(2)-sqrt(5))/2)*sign(p8[m,j])*min(abs(p8[m,j]),α)
   end

   return vcat(vec(p1),vec(p2),vec(p3),vec(p4),vec(p5),vec(p6),vec(p7),vec(p8))
end


text:
"""
    WaveSurface <: NeutralLandscapeMaker

    WaveSurface(; direction=360rand(), periods=1)
    WaveSurface(direction, [periods=1])

Creates a sinusoidal landscape with a `direction` and a number of `periods`. If
neither are specified, there will be a single period of random direction.
"""
@kwdef struct WaveSurface <: NeutralLandscapeMaker
    direction::Float64 = 360rand()
    periods::Int64 = 1
    function WaveSurface(direction::T, periods::K = 1) where {T <: Real, K <: Integer}
        @assert periods >= 1
        new(mod(direction, 360.0), periods)
    end
end

function _landscape!(mat, alg::WaveSurface)
    rand!(mat, PlanarGradient(alg.direction))
    mat .= sin.(mat .* (2π * alg.periods))
end


text:
# Selecting.

mutable struct SelectNode <: TabularNode
    over::Union{SQLNode, Nothing}
    list::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}

    function SelectNode(; over = nothing, list, label_map = nothing)
        if label_map !== nothing
            return new(over, list, label_map)
        end
        n = new(over, list, OrderedDict{Symbol, Int}())
        for (i, l) in enumerate(n.list)
            name = label(l)
            if name in keys(n.label_map)
                err = DuplicateLabelError(name, path = [l, n])
                throw(err)
            end
            n.label_map[name] = i
        end
        n
    end
end

SelectNode(list...; over = nothing) =
    SelectNode(over = over, list = SQLNode[list...])

"""
    Select(; over; list)
    Select(list...; over)

The `Select` node specifies the output columns.

```sql
SELECT \$list...
FROM \$over
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Select(Get.person_id);

julia> print(render(q))
SELECT "person_1"."person_id"
FROM "person" AS "person_1"
```
"""
Select(args...; kws...) =
    SelectNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Select), pats::Vector{Any}) =
    dissect(scr, SelectNode, pats)

function PrettyPrinting.quoteof(n::SelectNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Select))
    if isempty(n.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.list, ctx))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::SelectNode) =
    label(n.over)

rebase(n::SelectNode, n′) =
    SelectNode(over = rebase(n.over, n′), list = n.list, label_map = n.label_map)


text:
module Deps

import Conda, JSON, TOML, Pkg, Dates, ..PythonCall

### META

const _meta_file = Ref("")

meta_file() = _meta_file[]

function load_meta()
    fn = meta_file()
    isfile(fn) ? open(JSON.parse, fn) : Dict{String,Any}()
end

save_meta(meta) = open(io -> JSON.print(io, meta), meta_file(), "w")

function get_meta(keys...)
    meta = load_meta()
    for key in keys
        if haskey(meta, key)
            meta = meta[key]
        else
            return
        end
    end
    meta
end

function set_meta(args...)
    length(args) < 2 && error("setmeta() takes at least 2 arguments")
    here = meta = load_meta()
    for key in args[1:end-2]
        here = get!(Dict{String,Any}, here, key)
    end
    here[args[end-1]] = args[end]
    save_meta(meta)
end

### CONDA

const _conda_env = Ref("")

conda_available() = !isempty(_conda_env[])

"""
    conda_env()

The path to the Conda environment in which Python dependencies are managed.
"""
function conda_env()
    ans = _conda_env[]
    isempty(ans) && error("Conda is not available")
    return ans
end

conda_script_dir() = Conda.script_dir(conda_env())

const CONDA_ENV_STACK = Vector{Dict{String,String}}()

function conda_activate()
    e = conda_env()
    push!(CONDA_ENV_STACK, copy(ENV))
    try
        # We run `conda shell.$shell activate $e` to get the shell script which activates
        # the environment, then parse its instructions (taking advantage of the fact that
        # conda outputs these instructions in quite a structured fashion) and make the
        # necessary modifications.
        #
        # Another approach is to use `conda shell.$shell+json activate $e` to get the
        # instructions as JSON, but as of July 2020 the latest version of conda on 32-bit
        # linux does not support +json.
        #
        # Another approach is to actually execute the activate script in a subshell and then
        # extract the resulting environment. This would capture environment variables set
        # in scripts, which are currently ignored, and so replicate the conda environment
        # precisely. Just need to work out how to get at the resulting environment.
        if Sys.iswindows()
            shell = "powershell"
            exportvarregex = r"^\$Env:([^ =]+) *= *\"(.*)\"$"
            setvarregex = exportvarregex
            unsetvarregex = r"^(Remove-Item +\$Env:/|Remove-Variable +)([^ =]+)$"
            runscriptregex = r"^\. +\"(.*)\"$"
        else
            shell = "posix"
            exportvarregex = r"^\\?export ([^ =]+)='(.*)'$"
            setvarregex = r"^([^ =]+)='(.*)'"
            unsetvarregex = r"^\\?unset +([^ ]+)$"
            runscriptregex = r"^\\?\. +\"(.*)\"$"
        end
        for line in eachline(Conda._set_conda_env(`$(Conda.conda) shell.$shell activate $e`))
            if (m = match(exportvarregex, line)) !== nothing
                @debug "conda activate export var" k=m.captures[1] v=m.captures[2]
                ENV[m.captures[1]] = m.captures[2]
            elseif (m = match(setvarregex, line)) !== nothing
                @debug "ignoring conda activate set var" k=m.captures[1] v=m.captures[2]
            elseif (m = match(unsetvarregex, line)) !== nothing
                @debug "conda activate unset var" k=m.captures[1]
                delete!(ENV, m.captures[1])
            elseif (m = match(runscriptregex, line)) !== nothing
                @debug "ignoring conda activate script" file=m.captures[1]
            elseif !isempty(strip(line))
                @warn "ignoring conda activate line" line
            end
        end
        @assert realpath(ENV["CONDA_PREFIX"]) == realpath(e)
        return
    catch
        conda_deactivate()
        rethrow()
    end
end

function conda_deactivate()
    env = pop!(CONDA_ENV_STACK)
    for k in collect(keys(ENV))
        delete!(ENV, k)
    end
    merge!(ENV, env)
    return
end

conda_run(cmd::Cmd) = Conda.runconda(cmd, conda_env())

conda_run_root(cmd::Cmd) = Conda.runconda(cmd)

### PYTHON

python_dir() = Conda.python_dir(conda_env())

python_exe() = joinpath(python_dir(), Sys.iswindows() ? "python.exe" : "python")

### PIP

pip_exe() = joinpath(conda_script_dir(), Sys.iswindows() ? "pip.exe" : "pip")

pip_enable() = Conda.pip_interop(true, conda_env())

function pip_run(cmd::Cmd)
    env = conda_env()
    @info("Running $(`pip $cmd`) in $(env) environment")
    run(Conda._set_conda_env(`$(pip_exe()) $cmd`, env))
    nothing
end

### RESOLVE

function can_skip_resolve()
    # resolve if the conda environment doesn't exist yet
    isdir(conda_env()) || return false
    # resolve if we haven't resolved before
    deps = get_meta("jldeps")
    deps === nothing && return false
    # resolve whenever the PythonCall version changes
    version = get(deps, "version", nothing)
    version === nothing && return false
    version != string(PythonCall.VERSION) && return false
    # resolve whenever any of the environments in the load_path changes
    timestamp = get(deps, "timestamp", nothing)
    timestamp === nothing && return false
    timestamp = max(timestamp, stat(meta_file()).mtime)
    load_path = get(deps, "load_path", nothing)
    load_path === nothing && return false
    load_path != Base.load_path() && return false
    for env in Base.load_path()
        proj = Base.env_project_file(env)
        dir = nothing
        if proj isa String
            dir = dirname(proj)
        elseif proj === true
            dir = env
        end
        dir === nothing && continue
        isdir(dir) || continue
        stat(dir).mtime < timestamp || return false
        fn = joinpath(dir, "PythonCallDeps.toml")
        if isfile(fn)
            stat(fn).mtime < timestamp || return false
        end
    end
    return true
end

"""
    user_deps_file()

The path to the `PythonCallDeps.toml` file in the active environment.
"""
function user_deps_file()
    for env in Base.load_path()
        proj = Base.env_project_file(env)
        if proj isa String
            return joinpath(dirname(proj), "PythonCallDeps.toml")
        elseif proj === true
            return joinpath(env, "PythonCallDeps.toml")
        end
    end
    error()
end

function deps_files()
    dirs = String[]
    for env in Base.load_path()
        proj = Base.env_project_file(env)
        if proj isa String
            push!(dirs, dirname(proj))
        elseif proj === true
            push!(dirs, env)
        end
    end
    for p in values(Pkg.dependencies())
        push!(dirs, p.source)
    end
    ans = String[]
    for dir in dirs
        fn = joinpath(dir, "PythonCallDeps.toml")
        if isfile(fn)
            push!(ans, fn)
        end
    end
    return ans
end

"""
    resolve(; create=true, force=false)

Resolve all Python dependencies.

If `create=true` then a new Conda environment is created and activated.
Otherwise, the existing one is updated.

By default, if no dependencies have actually changed, then resolving them is skipped.
Specify `force=true` to skip this check and force resolving dependencies.
"""
function resolve(; create=true, force=false)
    if !force && can_skip_resolve()
        # skip installing any dependencies
        create && conda_activate()

    else
        # find python dependencies
        all_deps_files = deps_files()
        conda_channels = String[]
        conda_packages = String[]
        pip_packages = String[]
        pip_indexes = String[]
        scripts = String[]
        for fn in all_deps_files
            @info "Found PythonCall dependencies at '$(fn)'"
            deps = TOML.parsefile(fn)
            if haskey(deps, "conda")
                if haskey(deps["conda"], "channels")
                    union!(conda_channels, deps["conda"]["channels"])
                end
                if haskey(deps["conda"], "packages")
                    union!(conda_packages, deps["conda"]["packages"])
                end
            end
            if haskey(deps, "pip")
                if haskey(deps["pip"], "indexes")
                    union!(pip_indexes, deps["pip"]["indexes"])
                end
                if haskey(deps["pip"], "packages")
                    union!(pip_packages, deps["pip"]["packages"])
                end
            end
            if haskey(deps, "script")
                if haskey(deps["script"], "expr")
                    push!(scripts, deps["script"]["expr"])
                end
                if haskey(deps["script"], "file")
                    push!(scripts, read(deps["script"]["file"], String))
                end
            end
        end

        # canonicalise dependencies
        sort!(unique!(conda_channels))
        sort!(unique!(conda_packages))
        sort!(unique!(pip_packages))
        sort!(unique!(pip_indexes))
        sort!(unique!(scripts))

        # determine if any dependencies have changed
        env = conda_env()
        skip = !force && isdir(env)
        if skip
            depinfo = get_meta("jldeps")
            skip &= (
                conda_channels == get(depinfo, "conda_channels", nothing) &&
                conda_packages == get(depinfo, "conda_packages", nothing) &&
                pip_indexes == get(depinfo, "pip_indexes", nothing) &&
                pip_packages == get(depinfo, "pip_packages", nothing) &&
                scripts == get(depinfo, "scripts", nothing)
            )
        end

        # create and activate the conda environment with the desired packages
        # if update=true, just install the packages
        if skip
            # nothing has changed
            create && conda_activate()

        else
            conda_args = String[]
            for channel in conda_channels
                push!(conda_args, "--channel", channel)
            end
            append!(conda_args, conda_packages)
            if create || !isdir(env)
                ispath(env) && conda_run_root(`env remove --yes --prefix $env`)
                conda_run_root(`create --yes --no-default-packages --no-channel-priority --prefix $env $conda_args`)
                conda_activate()
            else
                conda_run(`install --yes --no-channel-priority --satisfied-skip-solve $conda_args`)
            end

            # install pip packages
            if !isempty(pip_packages)
                pip_enable()
                pip_args = String[]
                for index in pip_indexes
                    push!(pip_args, "--extra-index-url", index)
                end
                append!(pip_args, pip_packages)
                pip_run(`install $pip_args`)
            end

            # run scripts
            for script in scripts
                @info "Executing `$script`"
                eval(Meta.parse(script))
            end
        end

        # record what we did
        depinfo = Dict(
            "timestamp" => time(),
            "load_path" => Base.load_path(),
            "version" => string(PythonCall.VERSION),
            "files" => all_deps_files,
            "conda_packages" => conda_packages,
            "conda_channels" => conda_channels,
            "pip_packages" => pip_packages,
            "pip_indexes" => pip_indexes,
            "scripts" => scripts,
        )
        set_meta("jldeps", depinfo)
    end

    return
end

### INTERACTIVE

function spec_split(spec)
    spec = strip(spec)
    i = findfirst(c->isspace(c) || c in ('!','<','>','='), spec)
    if i === nothing
        return (spec, "")
    else
        return (spec[1:prevind(spec,i)], spec[i:end])
    end
end

"""
    status()

Display the status of dependencies of the current Julia project.
"""
function status()
    file = user_deps_file()
    printstyled("Status ", bold=true)
    println(file)
    if !isfile(file)
        println("(does not exist yet)")
        return
    end
    deps = TOML.parsefile(file)
    if haskey(deps, "conda")
        if haskey(deps["conda"], "channels") && !isempty(deps["conda"]["channels"])
            printstyled("Conda channels:", bold=true)
            println()
            for channel in deps["conda"]["channels"]
                println("  ", channel)
            end
        end
        if haskey(deps["conda"], "packages") && !isempty(deps["conda"]["packages"])
            printstyled("Conda packages:", bold=true)
            println()
            for spec in deps["conda"]["packages"]
                pkg, ver = spec_split(spec)
                print("  ", pkg, " ")
                printstyled(ver, color=:light_black)
                println()
            end
        end
    end
    if haskey(deps, "pip")
        if haskey(deps["pip"], "indexes") && !isempty(deps["pip"]["indexes"])
            printstyled("Pip indexes:", bold=true)
            println()
            for channel in deps["pip"]["indexes"]
                println("  ", channel)
            end
        end
        if haskey(deps["pip"], "packages") && !isempty(deps["pip"]["packages"])
            printstyled("Pip packages:", bold=true)
            println()
            for spec in deps["pip"]["packages"]
                pkg, ver = spec_split(spec)
                print("  ", pkg, " ")
                printstyled(ver, color=:light_black)
                println()
            end
        end
    end
    if haskey(deps, "script")
        if haskey(deps["script"], "expr")
            printstyled("Script expr:", bold=true)
            println()
            println("  ", deps["script"]["expr"])
        end
        if haskey(deps["script"], "file")
            printstyled("Script file:", bold=true)
            println()
            println("  ", deps["script"]["file"])
        end
    end
end

function getconvert!(::Type{T}, x, k) where {T}
    v = get!(T, x, k)
    if v isa T
        return v::T
    else
        v = convert(T, v)
        x[k] = v
        return x[k]::T
    end
end

"""
    add(...)

Add Python dependencies to the current Julia project.

Keyword arguments (all optional):
- `conda_channels`: An iterable of conda channels to use.
- `conda_packages`: An iterable of conda packages to install.
- `pip_indexes`: An iterable of pip indexes to use.
- `pip_packages`: An iterable of pip packages to install.
- `script_expr`: An expression to evaluate in the `Deps` module.
- `script_file`: The path to a Julia file to evaluate in the `Deps` module.
- `resolve=true`: When true, immediately resolve the dependencies. Otherwise, the
  dependencies are not resolved until you call [`resolve`](@ref) or load PythonCall in a
  new Julia session.
- `create=false`: When true, creates the environment from scratch when resolving.

The conda and pip packages can include version specifiers, such as `python>=3.6`.
"""
function add(; conda_channels=nothing, conda_packages=nothing, pip_indexes=nothing, pip_packages=nothing, script_expr=nothing, script_file=nothing, resolve=true, create=false)
    file = user_deps_file()
    deps = isfile(file) ? TOML.parsefile(file) : Dict{String,Any}()
    if conda_channels !== nothing
        sort!(union!(getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "conda"), "channels"), conda_channels))
    end
    if conda_packages !== nothing
        sort!(union!(getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "conda"), "packages"), conda_packages))
    end
    if pip_indexes !== nothing
        sort!(union!(getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "pip"), "indexes"), pip_indexes))
    end
    if pip_packages !== nothing
        sort!(union!(getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "pip"), "packages"), pip_packages))
    end
    if script_expr !== nothing
        getconvert!(Dict{String,Any}, deps, "script")["expr"] = script_expr
    end
    if script_file !== nothing
        getconvert!(Dict{String,Any}, deps, "script")["file"] = script_file
    end
    open(io->TOML.print(io, deps), file, "w")
    if resolve
        Deps.resolve(force=true, create=create)
    end
    return
end

"""
    rm(...)

Remove Python dependencies from the current Julia project.

Keyword arguments (all optional):
- `conda_channels`: An iterable of conda channels to remove.
- `conda_packages`: An iterable of conda packages to remove.
- `pip_indexes`: An iterable of pip indexes to remove.
- `pip_packages`: An iterable of pip packages to remove.
- `script_expr=false`: When true, remove the script expression.
- `script_file=false`: When true, remove the script file.
- `resolve=true`: When true, immediately resolve the dependencies. Otherwise, the
  dependencies are not resolved until you call [`resolve`](@ref) or load PythonCall in a
  new Julia session.
- `create=true`: When true, creates the environment from scratch when resolving.
"""
function rm(; conda_channels=nothing, conda_packages=nothing, pip_indexes=nothing, pip_packages=nothing, script_expr=false, script_file=false, resolve=true, create=true)
    file = user_deps_file()
    deps = isfile(file) ? TOML.parsefile(file) : Dict{String,Any}()
    if conda_channels !== nothing
        filter!(x -> x ∉ conda_channels, getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "conda"), "channels"))
    end
    if conda_packages !== nothing
        filter!(x -> spec_split(x)[1] ∉ conda_packages, getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "conda"), "packages"))
    end
    if pip_indexes !== nothing
        filter!(x -> x ∉ pip_indexes, getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "pip"), "indexes"))
    end
    if pip_packages !== nothing
        filter!(x -> spec_split(x)[1] ∉ pip_packages, getconvert!(Vector{String}, getconvert!(Dict{String,Any}, deps, "pip"), "packages"))
    end
    if script_expr
        delete!(getconvert!(Dict{String,Any}, deps, "script"), "expr")
    end
    if script_file !== nothing
        delete!(getconvert!(Dict{String,Any}, deps, "script"), "file")
    end
    open(io->TOML.print(io, deps), file, "w")
    if resolve
        Deps.resolve(force=true, create=create)
    end
    return
end

end


text:
__precompile__(true)
module Nothings

using Compat

export allownothing, disallownothing, isnothing, nothings,
       NothingException, levels, coalesce

if VERSION >= v"0.7.0-DEV.2762"
    using Base: Nothing, nothing
else
    export Nothing, nothing

    """
        Nothing

    A type with no fields whose singleton instance [`nothing`](@ref) is used
    to represent nothing values.
    """
    struct Nothing end

    """
        nothing

    The singleton instance of type [`Nothing`](@ref) representing a nothing value.
    """
    const nothing = Nothing()

    Base.show(io::IO, x::Nothing) = print(io, "nothing")

end

"""
 NothingException(msg)

Exception thrown when a [`nothing`](@ref) value is encountered in a situation
where it is not supported. The error message, in the `msg` field
may provide more specific details.
"""
struct NothingException <: Exception
    msg::AbstractString
end

import Base: ==, !=, <, *, <=, !, +, -, ^, /, &, |, xor

Base.showerror(io::IO, ex::NothingException) =
    print(io, "NothingException: ", ex.msg)

isnothing(::Any) = false
isnothing(::Nothing) = true

Base.promote_rule(::Type{T}, ::Type{Nothing}) where {T} = Union{T, Nothing}
Base.promote_rule(::Type{T}, ::Type{Union{S,Nothing}}) where {T,S} = Union{promote_type(T, S), Nothing}
Base.promote_rule(::Type{T}, ::Type{Any}) where {T} = Any
Base.promote_rule(::Type{Any}, ::Type{Nothing}) = Any
Base.promote_rule(::Type{Nothing}, ::Type{Any}) = Any
Base.promote_rule(::Type{Nothing}, ::Type{Nothing}) = Nothing

Base.convert(::Type{Union{T, Nothing}}, x) where {T} = convert(T, x)

# Comparison operators
==(::Nothing, ::Nothing) = nothing
==(::Nothing, b) = nothing
==(a, ::Nothing) = nothing
# != must be defined explicitly since fallback expects a Bool
!=(::Nothing, ::Nothing) = nothing
!=(::Nothing, b) = nothing
!=(a, ::Nothing) = nothing
Base.isequal(::Nothing, ::Nothing) = true
Base.isequal(::Nothing, b) = false
Base.isequal(a, ::Nothing) = false
<(::Nothing, ::Nothing) = nothing
<(::Nothing, b) = nothing
<(a, ::Nothing) = nothing
Base.isless(::Nothing, ::Nothing) = false
Base.isless(::Nothing, b) = false
Base.isless(a, ::Nothing) = true
if VERSION < v"0.7.0-DEV.300"
    <=(::Nothing, ::Nothing) = nothing
    <=(::Nothing, b) = nothing
    <=(a, ::Nothing) = nothing
end

# Unary operators/functions
for f in (:(!), :(+), :(-), :(Base.identity), :(Base.zero), :(Base.one), :(Base.oneunit),
        :(Base.abs), :(Base.abs2), :(Base.sign),
        :(Base.acos), :(Base.acosh), :(Base.asin), :(Base.asinh), :(Base.atan), :(Base.atanh),
        :(Base.sin), :(Base.sinh), :(Base.cos), :(Base.cosh), :(Base.tan), :(Base.tanh),
        :(Base.exp), :(Base.exp2), :(Base.expm1), :(Base.log), :(Base.log10), :(Base.log1p),
        :(Base.log2), :(Base.exponent), :(Base.sqrt), :(Base.gamma), :(Base.lgamma),
        :(Base.iseven), :(Base.ispow2), :(Base.isfinite), :(Base.isinf), :(Base.isodd),
        :(Base.isinteger), :(Base.isreal), :(Base.isimag), :(Base.isnan), :(Base.isempty),
        :(Base.iszero), :(Base.transpose), :(Base.float))
    @eval $(f)(d::Nothing) = nothing
end

for f in (:(Base.zero), :(Base.one), :(Base.oneunit))
    @eval function $(f)(::Type{Union{T, Nothing}}) where T
        T === Any && throw(MethodError($f, (Any,)))
        $f(T)
    end
end

# Binary operators/functions

for f in (:(+), :(-), :(*), :(/), :(^),
        :(Base.div), :(Base.mod), :(Base.fld), :(Base.rem), :(Base.min), :(Base.max))
    @eval ($f)(::Nothing, ::Nothing) = nothing
end

for f in (:(+), :(*), :(Base.min), :(Base.max))
    @eval begin
        ($f)(a::Number, o::Nothing) = a
        ($f)(o::Nothing, a::Number) = a
    end
end

Base.:(-)(a::Number, o::Nothing) =  a
Base.:(-)(o::Nothing, a::Number) = -a
Base.:(^)(a::Number, o::Nothing) =  a
Base.:(^)(o::Nothing, a::Number) = nothing

for f in (:(/), :(Base.div), :(Base.fld), :(Base.cld))
    @eval begin
        ($f)(a::Number, o::Nothing) = a
        ($f)(o::Nothing, a::Number) = inv(a)
    end
end

for f in (:(Base.mod), :(Base.rem))
    @eval begin
        ($f)(a::Number, o::Nothing) = a
        ($f)(o::Nothing, a::Number) = nothing
    end
end

# to avoid ambiguity warnings
(^)(::Nothing, ::Integer) = nothing

# Bit operators
(&)(::Nothing, ::Nothing) = nothing
(&)(a::Nothing, b::Bool) = b
(&)(b::Bool, a::Nothing) = b
(&)(::Nothing, b::Integer) = b
(&)(b::Integer, ::Nothing) = b
(|)(::Nothing, ::Nothing) = nothing
(|)(a::Nothing, b::Bool) = b
(|)(b::Bool, a::Nothing) = b
(|)(::Nothing, b::Integer) = b
(|)(b::Integer, ::Nothing) = b
xor(::Nothing, ::Nothing) = nothing
xor(a::Nothing, b::Bool) = b
xor(b::Bool, a::Nothing) = b
xor(a::Nothing, b::Integer) = b
xor(b::Integer, a::Nothing) = b

# String functions
*(d::Nothing, x::AbstractString) = x
*(d::AbstractString, x::Nothing) = x

end # module Nothings