text:
using Test
using NQCDynamics
using NQCDynamics.Calculators
using LinearAlgebra: tr, Diagonal

@testset "Adiabatic" begin
    model = NQCModels.Free(3)
    @test Calculators.AdiabaticCalculator{Float64}(model, 5) isa Calculators.AdiabaticCalculator
    @test Calculators.RingPolymerAdiabaticCalculator{Float64}(model, 5, 10) isa Calculators.RingPolymerAdiabaticCalculator

    calc = Calculators.AdiabaticCalculator{Float64}(model, 5)
    Calculators.evaluate_potential!(calc, rand(3, 5))
    Calculators.evaluate_derivative!(calc, rand(3, 5))

    calc = Calculators.RingPolymerAdiabaticCalculator{Float64}(model, 5, 10)
    Calculators.evaluate_potential!(calc, rand(3, 5, 10))
    Calculators.evaluate_derivative!(calc, rand(3, 5, 10))
end

@testset "Diabatic" begin
    model = NQCModels.DoubleWell()
    @test Calculators.DiabaticCalculator{Float64}(model, 5) isa Calculators.DiabaticCalculator
    @test Calculators.RingPolymerDiabaticCalculator{Float64}(model, 5, 10) isa Calculators.RingPolymerDiabaticCalculator

    calc = Calculators.DiabaticCalculator{Float64}(model, 1)
    Calculators.evaluate_potential!(calc, rand(1, 1))
    Calculators.evaluate_derivative!(calc, rand(1, 1))
    Calculators.eigen!(calc)
    Calculators.transform_derivative!(calc)
    Calculators.evaluate_nonadiabatic_coupling!(calc)
    Calculators.update_electronics!(calc, rand(1, 1))

    calc = Calculators.RingPolymerDiabaticCalculator{Float64}(model, 1, 10)
    Calculators.evaluate_potential!(calc, rand(1, 1, 10))
    Calculators.evaluate_centroid_potential!(calc, rand(1, 1))

    Calculators.evaluate_V̄!(calc)
    @test calc.V̄[1] ≈ tr(calc.potential[1]) / nstates(model)
    @test (@allocated Calculators.evaluate_V̄!(calc)) == 0

    Calculators.evaluate_traceless_potential!(calc)
    @test calc.traceless_potential[1] ≈ calc.potential[1] .- Diagonal(fill(calc.V̄[1], nstates(model)))
    @test tr(calc.traceless_potential[1]) ≈ 0 atol=eps(Float64)
    @test (@allocated Calculators.evaluate_traceless_potential!(calc)) == 0

    Calculators.evaluate_derivative!(calc, rand(1, 1, 10))
    Calculators.evaluate_centroid_derivative!(calc, rand(1, 1))

    Calculators.evaluate_D̄!(calc)
    @test calc.D̄[1] ≈ tr(calc.derivative[1]) / nstates(model)
    @test (@allocated Calculators.evaluate_D̄!(calc)) == 0

    Calculators.evaluate_traceless_derivative!(calc)
    @test calc.traceless_derivative[1] ≈ calc.derivative[1] .- Diagonal(fill(calc.D̄[1], nstates(model)))
    @test tr(calc.traceless_derivative[1]) ≈ 0 atol=eps(Float64)
    @test (@allocated Calculators.evaluate_traceless_derivative!(calc)) == 0

    Calculators.evaluate_traceless_adiabatic_derivative!(calc)
    @test calc.traceless_adiabatic_derivative[1] ≈ calc.adiabatic_derivative[1] .- Diagonal(fill(calc.D̄[1], nstates(model)))
    @test (@allocated Calculators.evaluate_traceless_adiabatic_derivative!(calc)) == 0

    Calculators.eigen!(calc)
    Calculators.transform_derivative!(calc)
    Calculators.evaluate_nonadiabatic_coupling!(calc)
    Calculators.update_electronics!(calc, rand(1, 1, 10))
end

@testset "LargeDiabatic" begin
    model = NQCModels.MiaoSubotnik()
    @test Calculators.LargeDiabaticCalculator{Float64}(model, 5) isa Calculators.LargeDiabaticCalculator

    calc = Calculators.LargeDiabaticCalculator{Float64}(model, 1)
    Calculators.evaluate_potential!(calc, rand(1, 1))
    Calculators.evaluate_derivative!(calc, rand(1, 1))
    Calculators.eigen!(calc)
    Calculators.transform_derivative!(calc)
    Calculators.evaluate_nonadiabatic_coupling!(calc)
    Calculators.update_electronics!(calc, rand(1, 1))
end

@testset "General constructors" begin
    model = NQCModels.DoubleWell()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.DiabaticCalculator
    @test Calculators.Calculator(model, 1, 2, Float64) isa Calculators.RingPolymerDiabaticCalculator
    model = NQCModels.Free()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.AdiabaticCalculator
    @test Calculators.Calculator(model, 1, 2, Float64) isa Calculators.RingPolymerAdiabaticCalculator
    model = NQCModels.MiaoSubotnik()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.LargeDiabaticCalculator
end

# for field ∈ (:potential, :derivative)
#     @testset "get_$field" begin
#         get_field = Symbol(:get, :_, field)
#         evaluate_field = Symbol(:evaluate, :_, field, :!)
#         eval(quote
#             model = NQCModels.DoubleWell()
#             calc = Calculators.DiabaticCalculator{Float64}(model, 1)
#             r = hcat(1.0)
#             @test !(Calculators.position(calc) ≈ r)
#             Calculators.$get_field(calc, r)
#             @test Calculators.position(calc) ≈ r
#             val1 = Calculators.$get_field(calc, r)
#             val2 = Calculators.$evaluate_field(calc, r)
#             @test val1 === val2 === calc.$field
#         end)
#     end
# end


text:
using DifferentialEquations
using DataFrames
using CSV
using LaTeXStrings
using Statistics
using GLM
using StatsPlots
using Optim
using ForwardDiff
using ColorSchemes
using Gadfly
using Cairo
using Fontconfig
g = Gadfly

#
pyplot(grid=false)
#
fnt = Plots.font("Arial")
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)

crec = parse(Colorant, "#F5EBCE")
crel = parse(Colorant, "#DEE6F0")

function three_state_model(dc, c, k, t)
    dc[1] = k['a'] * c[2] - k['s'] * c[1]
    dc[2] = k['s'] * c[1] - k['a'] * c[2] - k['i'] * c[2] + k['p'] * c[3]
    dc[3] = k['i'] * c[2] - k['p'] * c[3]
end

function isoutofdomain(c, k, t)
    if (c[1]<=1) && (c[2]<=1) && (c[3]<=1) && (c[1]>=0) && (c[2]>=0) && (c[3]>=0)
        return false
    else
        return true
    end
end

function run_recruitment(k_time_pairs, get_irrev=true, c0=[1.0; 0.0; 0.0])
    t = fill(Float64, 0)
    u = fill(Float64[], 0)

    for pair in k_time_pairs
        k = pair[1]
        times = pair[2]
        prob = ODEProblem(three_state_model, c0, times, k)
        sol = solve(prob, isoutofdomain=isoutofdomain, saveat=0.001)
        t = [t; sol.t]
        u = [u; sol.u]
        c0 = sol.u[end]
    end
    as = [x[1] for x in u]
    rs = [x[2] for x in u]
    is = [x[3] for x in u]
    if get_irrev
        data = [as rs is]
    else
        data = [as rs]
    end

    return [t, data]
end

################################################################################
### S2C-D
################################################################################

function get_pulsed_silencing_curve(s=200)
    pulse_pairs = [[k_recruit, (0.0, 1.0)],
                   [k_release, (1.0, 2.0)],
                   [k_recruit, (2.0, 3.0)],
                   [k_release, (3.0, 4.0)],
                   [k_recruit, (4.0, 5.0)],
                   [k_release, (5.0, 6.0)],
                   [k_recruit, (6.0, 7.0)],
                   [k_release, (7.0, 8.0)],
                   [k_recruit, (8.0, 9.0)],
                   [k_release, (9.0, 10.0)],
                   [k_release, (10.0, 15.0)]
                   ]
    gpp = [[0.0, 1.0], [2.0, 3.0], [4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]
    bpp = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0], [9.0, 15.0]]

    psol = run_recruitment(pulse_pairs)

    pp = plot()
    for pair in gpp
        ymin = 0.005
        ymax = 1.05
        pp = plot!(pp,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crec,
                   linealpha = 0,
                   linewidth = 0,
                   label = "",
                   seriestype = :shape)
        #  pp = vspan!(pp, pair, label="", linecolor=crec, fillcolor=crec, linealpha=0, fillalpha=0)
    end
    for pair in bpp
        ymin = 0.005
        ymax = 1.05
        pp = plot!(pp,
                   [pair[1], pair[1], pair[2], pair[2]],
                   [ymin, ymax, ymax, ymin],
                   color = crel,
                   linealpha = 0,
                   linewidth = 0,
                   label = "",
                   seriestype = :shape)
        #  pp = vspan!(pp, pair, label="", linecolor=crel, fillcolor=crel, linealpha=0)
    end
    #  for pair in gpp
    #      pp = vspan!(pp, pair, label="", color=crec)
    #  end
    #  for pair in bpp
    #      pp = vspan!(pp, pair, label="", color=crel)
    #  end

    pp = plot!(pp,
               psol[1], psol[2],
               label=["A" "R" "I"],
               color = [:blue :orange :red],
               xlabel="Time (days)",
               ylabel="Fraction of Cells",
               yticks = 0:0.5:1.0,
               xticks = 0:5:15,
               #  xlims = (-Inf, 15),
               #  ylims=(-Inf, 1.0),
               linewidth=4,
               legend = :outerright,
               order = 0,
               size=(500, s))
    return pp
end

function get_pulsed_vs_continuous(s=200)
    pulse_pairs = [[k_recruit, (0.0, 1.0)],
                   [k_release, (1.0, 2.0)],
                   [k_recruit, (2.0, 3.0)],
                   [k_release, (3.0, 4.0)],
                   [k_recruit, (4.0, 5.0)],
                   [k_release, (5.0, 6.0)],
                   [k_recruit, (6.0, 7.0)],
                   [k_release, (7.0, 8.0)],
                   [k_recruit, (8.0, 9.0)],
                   [k_release, (9.0, 10.0)],
                   [k_release, (10.0, 15.0)]
                   ]
    gpp = [[0.0, 1.0], [2.0, 3.0], [4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]

    cont_pairs = [[k_release, (0.0, 5.0)],
                  [k_recruit, (5.0, 10.0)],
                  [k_release, (10.0, 15.0)]]
    gpc = [[5.0, 10.0]]

    psol = run_recruitment(pulse_pairs)
    csol = run_recruitment(cont_pairs)

    pc = plot(
              psol[1], psol[2][:, 3],
              label="Pulsed",
              color=:purple,
              xlabel="Time (days)",
              ylabel="I(t)",
              ylims=(0.0, 1.0),
              yticks = 0:0.5:1.0,
              linewidth=3,
              legend = :topright,
              size=(500, 300))
    pc = plot!(pc,
               csol[1], csol[2][:, 3],
               label="Continuous",
               color=:deepskyblue,
               linewidth=3)

    for pair in gpp
        pc = vspan!(pc, pair, label="", color=:violet, alpha=0.2)
    end
    for pair in gpc
        pc = vspan!(pc, pair, label="", color=:cyan, alpha=0.2)
    end
    return pc
end

ksi = 1
ksa = 1

ks = 1

k_recruit = Dict('a' => ks/ksa,
                 's' => ks,
                 'i' => ks/ksi,
                 'p' => 0)

k_release = Dict('a' => ks/ksa,
                 's' => 0,
                 'i' => 0,
                 'p' => 0)

pp = get_pulsed_silencing_curve(500)
pc = get_pulsed_vs_continuous(500)

savefig(pp, "fig_s2c.pdf")
savefig(pc, "fig_s2d.pdf")

################################################################################
### S2E
################################################################################

function get_max_delta(ksa, ksi)
    mdelta = 0
    k_recruit = Dict('a' => 0.3,
                     's' => 0.3 * ksa,
                     'i' => 0.3 * ksa / ksi,
                     'p' => 0)
    k_release = Dict('a' => 0.3,
                     's' => 0,
                     'i' => 0,
                     'p' => 0)
    for t in collect(0.01:0.01:4.99)
        kt_pairs = [[k_recruit, (0.0, t)],
                    [k_release, (t, 5.0)]]
        sol_on = run_recruitment(kt_pairs)
        sol_off = run_recruitment(kt_pairs, true, [0.0; 1.0; 0.0])
        vals_on = sol_on[2][:,3]
        vals_off = sol_off[2][:,3]
        deltas = vals_off - vals_on
        md = maximum(deltas)
        mdelta = maximum([md, mdelta])
    end
    return mdelta
end

function analytical_max_delta(ka, ks, ki)
    function irrev_delta(t, ka, ks, ki)
        term1 = -0.5 * (ka + ki + ks + sqrt(-4 * ka * ki + (ka + ki + ks)^2)) * t
        term2 = (-1 + exp(sqrt(-4 * ka * ki + (ka + ki + ks)^2) * t)) * ki
        term3 = sqrt(-4 * ka * ki + (ka + ki + ks)^2)
        return (exp(term1) * term2)/term3
    end

    function tmax_delta(ka, ks, ki)
        term1 = -ka - ki - ks - sqrt(-4 * ki * ks + (ka + ki + ks)^2)
        term2 = -ka - ki - ks + sqrt(-4 * ki * ks + (ka + ki + ks)^2)
        term3 = sqrt(-4 * ki * ks + (ka + ki + ks)^2)
        return log(term1/term2)/term3
    end

    return irrev_delta(tmax_delta(ka, ks, ki), ka, ks, ki)
end

ksa_vals = zeros(0)
ksi_vals = zeros(0)
md_vals = zeros(0)

for ksa in 0.01:0.1:10.01
    for ksi in 0.01:0.1:10.01
        md = analytical_max_delta(1.0, 1.0 * ksa, 1.0 * ksa/ksi)
        push!(ksa_vals, ksa)
        push!(ksi_vals, ksi)
        push!(md_vals, md)
    end
end

data = DataFrame(Dict("Ratio of Reversible Silencing to Activation" => ksa_vals,
                      "Ratio of Reversible to Irreversible Silencing" => ksi_vals,
                      "Gap" => md_vals))

println("plotting")
p = g.plot(data,
         x="Ratio of Reversible to Irreversible Silencing",
         y="Ratio of Reversible Silencing to Activation",
         color="Gap", g.Geom.rectbin,
         g.Guide.yticks(ticks=collect(0:10)), g.Guide.xticks(ticks=collect(0:10)),
         g.Guide.xlabel("ks/ka ratio"), g.Guide.ylabel("ks/ki ratio"),
         g.Coord.cartesian(xmin=0, xmax=10, ymin=0, ymax=10, aspect_ratio=1),
         g.Scale.ContinuousColorScale(p -> get(ColorSchemes.viridis, p)),
         g.Theme(major_label_font="ArialMT", minor_label_font="ArialMT",
               key_title_font="ArialMT", key_label_font="ArialMT", key_label_color="black",
               grid_color="black",
               minor_label_color="black", major_label_color="black", key_title_color="black"))

img = PDF("fig_s2e.pdf")
draw(img, p)

ext: jl


text:
# Note: Run this script without using `mpirun`

using MPI

manager = MPIManager(np=4)
addprocs(manager)

println("Added procs $(procs())")

println("Running 01-hello as part of a Julia cluster")
@mpi_do manager (include("01-hello-impl.jl"); do_hello())

# Interspersed julia parallel call
nheads = @parallel (+) for i=1:10^8
  Int(rand(Bool))
end
println("@parallel nheads $nheads")

println("Running 02-broadcast as part of a Julia cluster")
@mpi_do manager (include("02-broadcast-impl.jl"); do_broadcast())

M = [rand(10,10) for i=1:10]
pmap(svd, M)
println("pmap successful")

println("Running 03-reduce as part of a Julia cluster")
@mpi_do manager (include("03-reduce-impl.jl"); do_reduce())

pids = [remotecall_fetch(myid, p) for p in workers()]
println("julia pids $pids")

println("Running 04-sendrecv as part of a Julia cluster")
@mpi_do manager (include("04-sendrecv-impl.jl"); do_sendrecv())

println("Exiting")
exit()


text:
"""
  CES

CES aggregator of the form

`Y = A * [sum of alpha (X ^ rho)] ^ (1/rho)`
`Q = sum of (alpha X) ^ rho`

Sum of `alpha` = 1. Which is why the `alpha` are not taken to the power `rho`.

Accommodates matrix inputs (T x N)

Best to store AV in the object
- little cost in terms of efficiency
- consistency ensured
"""
mutable struct CES{F1} <: AbstractProductionFunction{F1}
  substElast :: F1
  # relative weights; Nx1
  alphaV :: Vector{F1}
  # Productivities
  AV :: Vector{F1}
end

Base.show(io :: IO,  fS :: CES{F1}) where F1 =
  print(io, "CES production function with $(n_inputs(fS)) inputs and elasticity $(subst_elast(fS))");

"""
	$(SIGNATURES)

Substitution elasticity.
"""
subst_elast(fS :: CES{F1}) where F1 = fS.substElast;

"""
	$(SIGNATURES)

Factor neutral productivities.
"""
productivities(fS :: CES{F1}) where F1 = fS.AV;

"""
	$(SIGNATURES)

Number of inputs.
"""
n_inputs(fS :: CES{F1}) where F1 = length(fS.alphaV);

"""
	$(SIGNATURES)

Does the production function have constant returns to scale?
"""
constant_returns(fS :: CES{F1}) where F1 = true;

# function CES(sElast :: F1, aV :: Vector{F1}, prodV :: Vector{F1}) where F1 <: AbstractFloat

#     fS = new();
#     set_params(fS, substElast = sElast, alphaV = aV, AV = prodV);
#     return fS
#   end
# end

"""
	$(SIGNATURES)

Validate a CES production function.
"""
function validate_prod_fct(this :: CES{F1}) where F1 <: AbstractFloat
    # T, n = size(fS.xM);
    isValid = check_subst_elast(this.substElast)
    # check more +++
end

function check_subst_elast(substElast :: F1) where
  F1 <: AbstractFloat
    return (substElast > 0.01  &&  substElast < 50)
end

## Set parameters in CES
# function set_params(fS :: CES;  substElast = [], alphaV = [], AV = [])
#   # Assign, then check consistency
#   if !isempty(substElast)
#     fS.substElast = substElast;
#   end
#   if !isempty(alphaV)
#     fS.alphaV = alphaV;
#   end
#   if !isempty(AV)
#     fS.AV = AV;
#   end
#   # if !isempty(xM)
#   #   fS.xM = xM;
#   # end

#   # Derived
#   fS.rho = 1.0 - 1.0 / fS.substElast;

#   check(fS);
# end

# # ----  return rho from subst elast
curvature(fS :: CES{F1}) where F1 <: AbstractFloat =
  curvature_from_subst_elast(fS.substElast);

curvature_from_subst_elast(substElast :: F1) where F1 <: Real =
  one(F1) / substElast - one(F1);

## -------  Returns sum(alpha * x ^ rho) ^ (1/rho)
# With rho = 0 (Cobb-Douglas) this returns 1.0.
function ces_q(alphaV :: AbstractVector{F1}, xM :: AbstractMatrix{F1},
  rho :: F1) where F1 <: AbstractFloat

  T = size(xM, 1);
  qV = zeros(F1, T);
  for t = 1 : T
      # Careful here: alphaV * xM could become a matrix
      @views qV[t] = ces_q(alphaV, xM[t,:], rho);
  end
  return qV
end

# One set of inputs.
function ces_q(alphaV :: AbstractVector{F1}, xV :: AbstractVector{F1}, rho :: F1) where F1
  return sum(alphaV .* (xV .^ rho));
end

"""
	$(SIGNATURES)

Output given inputs. If substitution elasticity is close to 1, switch to Cobb-Douglas.
"""
output(fS :: CES{F1}, xM :: AbstractMatrix{F1}) where F1 <: AbstractFloat =
  productivities(fS) .* ces_output(xM, fS.alphaV, subst_elast(fS));
#   qV = ces_q(fS.alphaV, xM, curvature(fS));
#   # This does not overflow when curvature = 1.0
#   yV = productivities(fS) .* (qV .^ (one(F1) / curvature(fS)));
#   if switch_cobb_douglas(subst_elast(fS), yV)
#     # Handle numerical overflow by switching to Cobb-Douglas
#     yV = productivities(fS) .* output_cobb_douglas(xM, fS.alphaV);
#   end
#   return yV
# end

"""
	$(SIGNATURES)

Output from CES production function with weights `alphaV` and substitution elasticity `substElast`. Returns a scalar if xM is a Vector, but not if it is a Matrix with 1 row.

# Arguments

- xM: Can be a `Matrix` where each row is a set of inputs or a `Vector` for a single set of inputs.
"""
function ces_output(xM :: AbstractArray{F1}, alphaV :: AbstractVector{F1}, substElast :: F1) where F1 <: Real

  curv = curvature_from_subst_elast(substElast);
  # This does not overflow when curvature = 1.0
  # If elasticity is close to 1, some yV are Inf. This is checked in the next step.
  yV = ces_q(alphaV, xM, curv) .^ (one(F1) / curv);
  if switch_cobb_douglas(substElast, yV)
    # Handle numerical overflow by switching to Cobb-Douglas
    yV = output_cobb_douglas(xM, alphaV);
  end
  return yV
end

# Decide whether Cobb-Douglas needs to be used
function switch_cobb_douglas(substElast, xV)
  if abs(substElast - 1.0) < 0.01
    return true
  elseif (abs(substElast - 1.0) < 0.05)  &&  any(isinf.(xV))
    return true
  else
    return false
  end
end

"""
	$(SIGNATURES)

Marginal products (TxN).
"""
mproducts(fS :: CES{F1}, xM :: AbstractMatrix{F1}) where F1 <: Real =
  productivities(fS) .* ces_mproducts(xM, fS.alphaV, subst_elast(fS));

#   T = length(fS.AV)
#   N = length(fS.alphaV);
#   rho = curvature(fS);
#   # qV is T x 1
#   qV = ces_q(fS.alphaV, xM, rho);
#   aqV = fS.AV .* (qV .^ (one(F1) / rho - one(F1)));
#   # mpM = similar(xM);
#   # for j = 1 : N
#   #   mpM[:,j] = aqV .* (fS.alphaV[j] .^ rho) .* (xM[:,j] .^ (rho - 1));
#   # end
#   if switch_cobb_douglas(subst_elast(fS), aqV)
#     mpM = fS.AV .* mproducts_cobb_douglas(xM, fS.alphaV);
#   else
#     mpM = aqV .* (xM .^ (rho-one(F1))) .* fS.alphaV';
#   end
#   # @assert (all(abs.(mp2M - mpM) .< 1e-5))
#   return mpM
# end

"""
	$(SIGNATURES)

CES marginal products (TxN).
"""
function ces_mproducts(xM :: AbstractMatrix{F1}, alphaV :: AbstractVector{F1}, substElast :: F1) where F1 <: Real

  T, N = size(xM);
  rho = curvature_from_subst_elast(substElast);
  # qV is T x 1
  qV = ces_q(alphaV, xM, rho);
  aqV = (qV .^ (one(F1) / rho - one(F1)));
  if switch_cobb_douglas(substElast, aqV)
    mpM = mproducts_cobb_douglas(xM, alphaV);
  else
    mpM = aqV .* (xM .^ (rho-one(F1))) .* alphaV';
  end
  return mpM
end

function make_test_ces(T, n, substElast)
  alphaV = collect(range(1.0, 2.0, length = n));
  alphaV = alphaV ./ sum(alphaV);
  if T == 1
    AV = [0.5];
  else
    AV = collect(LinRange(0.5, 9.0, T));
  end

  fS = CES(substElast, alphaV, AV);
  @assert validate_prod_fct(fS)
  return fS
end

# ------------


text:
using Fleck

@safetestset never_fair = "Never fair weather" begin
    using Fleck: Never

    n = Never()
    p = params(n)
    @test p == ()
    @test isa(n, Never)
end


text:
using Documenter, CompScienceMeshes

makedocs(
    clean=false,
    sitename = "CompScienceMeshes.jl Manual")
# deploydocs(
#     julia = "nightly",
#     deps = Deps.pip("mkdocs", "python-markdown-math"),
#     repo = "github.com/krcools/CompScienceMeshes.jl.git",)
deploydocs(
    repo = "github.com/krcools/CompScienceMeshes.jl.git",
)


text:
#####
##### Generate the file julia-latexsubs.el. Invoke from the shell as
#####
##### julia make-julia-latexsubs.jl
#####

@assert VERSION ≥ v"1"          # use a recent Julia version

import REPL

"""
Write Emacs lisp code that populates the hash table named `varname` to `dest`, using
key-value pairs from `src`.
"""
function write_latexsubs_hashtable(src, dest::IO, varname::AbstractString)
    for (k, v) in sort!(collect(src), by = last)
        ks = escape_string(k)
        vs = escape_string(v)
        if occursin(r"^\\U[0-9A-Fa-f]+$", vs)
            # codepoints outside the BMP can be problematic in older Emacsen
            cp = vs[3:end]
            println(dest, "(let ((c (decode-char 'ucs #x$cp)))\n",
                    "  (if c (puthash \"$ks\" (char-to-string c) $(varname))))")
        else
            println(dest, "(puthash \"$ks\" \"$vs\" $(varname))")
        end
    end
end

open("julia-latexsubs.el", "w") do io
    write_latexsubs_hashtable(REPL.REPLCompletions.latex_symbols, io, "julia-latexsubs")
end

@info "generated latex substitutions"

max_stars_repo_name: dvdhansen/julia-emacs

max_issues_repo_name: emmt/julia-emacs

max_forks_repo_name: emmt/julia-emacs