text:
struct PlumeEmulator{T <: Real} <: AbstractMeltRate
    α :: T         #calibration coefficient
    λ1 :: T        #liquidus slope
    λ2 :: T        #liquidus intercept
    λ3 :: T        #liquidus pressure coefficient
    E0 :: T        #entrainment coefficient
    Cd :: T        #drag coeffecient
    Γ_TS :: T      #combined Stanton number
    L :: T         #latent heat of fusion for ice
    c :: T         #specific heat capacity of ocean water
    βs :: T        #haline contraction coefficient
    βt :: T        #thermal expansion coefficient
    g :: T         #gravity
    ρi :: T        #ice density
    ρw :: T        #water density
    Ta :: Function #Function specifying ambient temperature
    Sa :: Function #Function specifying ambient salinity
end

"""
    PlumeEmulator(; <kwargs>)

Construct a plume emulator to prescribe the melt rate in WAVI

Keyword arguments
=================
- α: calibration coefficient: constant prefactor used to calibrate melt rate. Defaults to 0.73 according to Favier 2019 (10.5194/gmd-12-2255-2019)
- λ1: liquidus slope (units ∘C)
- λ2: liquidus intercept (units ∘C)
- λ3: liquidus pressure coefficient (units K/m)
- E0: entrainment coefficient (dimensionless)
- Cd: drag coefficient (dimensionless)
- Γ_TS: combined Stanton number (dimensionless)
- L: latent heat of fusion for ice (units: J/kg)
- c: specific heat capacity of water (units J/kg/K)
- βs: haline contraction coefficient (units psu^(-1))
- βt: thermal expansion coefficient (units C^{-1})
- g: gravitational acceleration (units m / s^2)
- ρi: ice density (units kg/m^3)
- ρw: water density (units kg/m^3)
- Ta: ambient temperature function, defaults to the ISOMIP warm0 scenario. Ta must be a function of a single variable (depth) [time dependence not yet implemented in WAVI.jl]
- Sa: ambient salnity function, defaults to the ISOMIP warm0 scenario.  Sa must be a function of a single variable (depth) [time dependence not yet implemented in WAVI.jl]
"""
function PlumeEmulator(;
                        α = 0.73,
                        λ1 = -5.73e-2,
                        λ2 = 8.32e-4,
                        λ3 = 7.61e-4,
                        E0 = 3.6e-2,
                        Cd = 2.5e-3,
                        Γ_TS = 0.0118,
                        L = 3.35e5,
                        c = 3.974e3,
                        βs = 7.86e-4,
                        βt = 3.87e-5,
                        g = 9.81,
                        ρi = 918.0,
                        ρw = 1028.0,
                        Ta = isomip_warm0_temp,
                        Sa = isomip_warm0_salinity)

    #check that Ta and Sa accept a single argument
    try Ta(0.0)
    catch
        ArgumentError("Ambient temperature function Ta must accept a single argument")
    end
    try Sa(0.0)
    catch
        ArgumentError("Ambient salinity function Sa must accept a single argument")
    end

    return PlumeEmulator(α, λ1, λ2, λ3,E0, Cd, Γ_TS, L, c, βs, βt,g, ρi, ρw, Ta, Sa)

end

function update_melt_rate!(plume_emulator::PlumeEmulator, fields, grid)
    @unpack basal_melt, h, b, grounded_fraction = fields.gh #get the ice thickness and grounded fraction
    set_plume_emulator_melt_rate!(basal_melt,
                                h,
                                grounded_fraction,
                                b,
                                grid.dx,
                                grid.dy,
                                plume_emulator)
    return nothing
end

"""
    set_plume_emulator_melt_rate(melt, h, grounded_fraction)

Set the melt rate in according to Lazeroms 2018. All plume parameters passed in the plume model emulator (pme)
"""
function set_plume_emulator_melt_rate!(melt, h, grounded_fraction,bathy, dx, dy,pme)
    nx, ny = size(h)
    zbf = @. -(pme.ρi/pme.ρw)*h.*(1-grounded_fraction) .+ grounded_fraction*bathy #ice draft if floating everywhere
    ∂zb∂x, ∂zb∂y = get_slope(zbf,dx, dy) #returns the partial derivates of base in both directions

    #loop over each grid point
    for i = 1:nx
        for j = 1:ny
            melt[i,j] =  pme_pointwise_melt(grounded_fraction,zbf, ∂zb∂x,∂zb∂y, pme, i, j)
        end
    end

    return melt
end

"""
    pme_pointwise_melt(grounded_fraction,zbf, ∂zb∂x,∂zb∂y, pme, i, j)

Return the pme prediction of melt rate at single grid point (i,j).
____________
Arguments:
- grounded_fraction: matrix of grounded fraction at grid points
- zbf: ice shelf basal elevation
- ∂zb∂x: slope of ice shelf draft at grid pts in x
- ∂zb∂y: slope of ice shelf draft at grid pts in y
- pme: plume model emulator structure. Holds all parameters.
- i,j: specific grid pt
"""
function pme_pointwise_melt(grounded_fraction,zbf, ∂zb∂x,∂zb∂y, pme, i, j)
    if grounded_fraction[i,j] == 0. #i.e. if its floating
        #initialize directions
        dirs = get_plume_directions() #array of directions
        nd = size(dirs)[1] #number of search directions
        norm_dirs = dirs ./ norm.(eachrow(dirs)) #normalized directions

        #initialize slope and gl depth
        zgl = zeros(1, nd) #grounding line depth
        s = zeros(1, nd) #basal slopes
        iscount = trues(1, nd) #do these points count or not?
        iscount_s = trues(1,nd)
        iscount_zgl = trues(1,nd)

        #loop over directions
        for k = 1:size(dirs)[1]
            #compute the slope in this direction
            s[k] = -dot(norm_dirs[k,:], [∂zb∂x[i,j] ∂zb∂y[i,j]]) #note negative sign
            if s[k] <= 0
                iscount[k] = false
                iscount_s[k] = false
            end #[this syntax count be tightened]

            #compute gl depth in this direction
            zgl[k] = get_pme_gl_depth(dirs[k,:], grounded_fraction, zbf, i, j)
            if isnan(zgl[k]) #no gl point found in this direction
                iscount[k] = false
                iscount_zgl[k] = false
            elseif zgl[k] >= zbf[i,j] #if gronding line higher than shelf point
                iscount[k] = false
                iscount_zgl[k] = false
            end
        end

        #compute effective gl depth and angle
        N = sum(iscount)
        zgl_eff = sum(zgl[iscount])/N;
        α = atan(sum(s[iscount])/N);

        #update the melt rate
        if N > 0
            Sa_gl = pme.Sa(zgl_eff)
            Sa_local = pme.Sa(zbf[i,j])
            Ta_local = pme.Ta(zbf[i,j])
            Tf0 = pme.λ1*Sa_gl + pme.λ2 + pme.λ3 * zgl_eff
            M0 = lazeroms_melt_rate_prefactor(pme.βs, pme.βt, pme.g, pme.Cd, pme.Γ_TS, pme.L, pme.c, pme.λ1, pme.λ3, pme.E0, α, Sa_local, Ta_local)
            x = lazeroms_dimensionless_distance(pme.λ3, zbf[i,j], zgl_eff, Ta_local, Tf0)
            melt = M0 * pme.α * (Ta_local - Tf0)^2 * lazeroms_dimensionless_melt(x) * 365.25* 24 * 60^2
        else
            melt = 0.
        end
    else
        melt = 0.
    end
    return melt #, zgl, s, iscount, iscount_s, iscount_zgl
end

"""
    get_slope(h, dx, dy)

Returns the partial derivatives of h with respect to x and y as matrices. Uses centred FD at interior points and one sided FD at edge points (both second order)

"""
function get_slope(h, dx, dy)
    ∂h∂x = zeros(size(h))
    ∂h∂x[2:end-1,:] = (-h[1:end-2,:] + h[3:end,:])/ 2/dx
    ∂h∂x[1,:] = (-3/2 *h[1,:] + 2*h[2,:] - 1/2 * h[3,:] )/dx
    ∂h∂x[end,:] = (1/2 * h[end-2, :] - 2 *h[end-1,:] + 3/2 * h[end,:])/dx

    ∂h∂y = zeros(size(h))
    ∂h∂y[:,2:end-1]= (-h[:,1:end-2] + h[:,3:end])/2 /dy
    ∂h∂y[:,1] = (-3/2 * h[:,1] + 2 * h[:,2] - 1/2 * h[:,3])/dy
    ∂h∂y[:,end] = (1/2 * h[:,end-2] - 2*h[:,end-1] + 3/2 * h[:,end] )/dy;

    return ∂h∂x, ∂h∂y
end

function get_plume_directions()
    directions = [0 1; #N
                  1 2; #NNE
                  1 1; #NE
                  2 1; #ENE
                  1 0; #E
                  2 -1; #ESE
                  1 -1; #SE
                  1 -2; #SSE
                  0 -1; #S
                  -1 -2; #SSW
                  -1 -1; #SW
                  -2 -1; #WSW
                  -1 0; #W
                  -2 1; #WNW
                  -1 1; #NW
                  -1 2] #NNW
    return directions
end

function lazeroms_melt_rate_prefactor(βs,βt, g, Cd, Γ_TS, L, c,λ1, λ3,  E0, α, Sa, Ta)
    cp1 = (L/c)/sqrt(Cd)/Γ_TS * βt/βs / Sa
    cp2 = -λ1 * βt/βs
    ct = cp2/cp1
    M0 = sqrt(βs * Sa * g) / sqrt(λ3 * (L/c)^3) *
         ((1 -  cp1 * sqrt(Cd)*Γ_TS)/(Cd + E0 * sin(α)))^(1/2) *
         (sqrt(Cd) * Γ_TS * E0 * sin(α)/ (sqrt(Cd) * Γ_TS + ct + E0 * sin(α)))^(3/2)
    return M0
end

"""
    lazeroms_dimensionless_distance(λ3, zb, zgl, Ta, Tf0)

The dimensionless distance in the Lazeroms plume emulator formulation (see equation 28b in Lazeroms et al. 2019)
Note that we include only the first term in the expansion (i.e. no contribution from term C_ϵ)
"""
lazeroms_dimensionless_distance(λ3, zb, zgl, Ta, Tf0) = λ3 * (zb - zgl) / (Ta - Tf0)

"""
    lazeroms_dimensionless_melt(x)

Dimensionless melt function, equation 26 in Lazeroms et al. 2019
"""
lazeroms_dimensionless_melt(x) = (3*(1-x)^(4/3) - 1)*(1- (1-x)^(4/3))^(1/2) / 2 / sqrt(2)

"""
    get_pme_gl_depth(direction, grounded_fraction, zbf, i,j)

Get the effective grounding line depth at shelf point i,j in the specified direction.
Note that this function does not yet apply the interpolation described in Lazeroms et al. 2018
"""
function get_pme_gl_depth(direction, grounded_fraction, zbf, i,j)
    nx, ny = size(grounded_fraction)

    #initialize
    zgl = NaN
    isterm_gl = false #have we found a grounded point?
    isterm_domain = false #have we gone outsdie of the domain?
    p = i
    q = j

    #workflow: (1) check we aren't grounded (before step as may be grounded to begin with)
    #          (2) make a step
    #          (3) check we haven't gone out of the domain
    #          (4) loop if still in domain and not grounded
    while ~(isterm_domain) && ~(isterm_gl)
        #check whether search point is grounded, trigger terminated if so
        grounded_fraction[p,q] > 1 - 1e-3 ?  isterm_gl = true : isterm_gl = false
        if isterm_gl #we've found a grounded point
            zgl = zbf[p,q]
        end

        #update search point
        p = p + direction[1]
        q = q + direction[2]

        #have we gone out of the domain
        isin_domain(p,q,nx,ny) ? isterm_domain = false : isterm_domain = true #if the first step is out of the domain, terminate
    end #end while loop
    return zgl
end

isin_domain(p,q,nx,ny) = ((p <= nx) && (p>0) && (q >0) && (q <= ny)) ? true : false


text:
include("phylo.jl")

function singly_cuts(f, t :: Tree)
    let root_clade = nothing
        function g(parent, child, root_clade)
            if issingly(root_clade - child._clades)
                f(reroot_root_clades(parent, child, root_clade))
            end
        end

        function h(t)
            if issingly(t)
                f(t)
            end
        end

        function cuts_(t :: Tree)
            if isleaf(t)

            else
                if (isroot(t))
                    root_clade = t._clades
                    h(t)
                end
                h(t.left)
                h(t.right)
                if !isroot(t)
                    g(t, t.left, root_clade)
                    g(t, t.right, root_clade)
                end
                cuts_(t.left)
                cuts_(t.right)
            end
        end
        cuts_(t)
    end
end

function istrivial(t :: Tree)
    (t._clades.data |> length) < 4
end

 function singly_cuts(t :: Tree)
    res = Tree[]
    singly_cuts(t) do s
        push!(res, s)
    end
    res
 end

ext: jl

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
# unit test he SVRG solver
include("../src/DMD_Util.jl")
include("../src/DMD_Type.jl")
include("../src/DMD_Func.jl")
include("../src/DMD_BFGS.jl")

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 8;      # number of modes
T = Float64;
sigman = 1e-2;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
Xdat = X + sigman*Random.randn!(copy(X))
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------

#
# Initial Guess
#------------------------------------------------------------------------------
# create initial guess using trap. rule with exact dmd
 a0, B0 = dmdexactestimate(m, n, k, Xdat, t);
#
# Create Object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad);
 copyto!(params.a, a0);
# copyto!(params.B, B0);
#Random.randn!(params.ar);
#
# Apply Solver
#------------------------------------------------------------------------------
sigma   = 1e0;
itm = 1000;
tol = 1e-10;
ptf = 1;
opts = DMD_BFGS_Options(sigma, params, itm=itm, tol=tol, ptf=ptf);

obj_his, err_his = solveDMD_withBFGS(params, opts);

println(besterrperm(params.a,at))

println(aBFunc(params)/lossFunc(X))

params2 = DMDParams(k,Xdat,t,lossFunc,lossGrad);
copyto!(params2.a,at)
println(aBFunc(params2)/lossFunc(X))
copyto!(params2.a,a0)
println(aBFunc(params2)/lossFunc(X))


text:
# Quadrature rule

"""
"""
struct CellQuadrature{DS} <: CellDatum
  cell_quad::AbstractArray{<:Quadrature}
  cell_point::AbstractArray{<:AbstractArray{<:Point}}
  cell_weight::AbstractArray{<:AbstractArray{<:Real}}
  trian::Triangulation
  domain_style::DS
end

"""
"""
function CellQuadrature(trian::Triangulation,degree::Integer)
  ctype_to_reffe, cell_to_ctype = compress_cell_data(get_cell_reffe(trian))
  ctype_to_quad = map(r->Quadrature(get_polytope(r),degree),ctype_to_reffe)
  ctype_to_point = map(get_coordinates,ctype_to_quad)
  ctype_to_weight = map(get_weights,ctype_to_quad)
  cell_quad = expand_cell_data(ctype_to_quad,cell_to_ctype)
  cell_point = expand_cell_data(ctype_to_point,cell_to_ctype)
  cell_weight = expand_cell_data(ctype_to_weight,cell_to_ctype)
  CellQuadrature(cell_quad,cell_point,cell_weight,trian,ReferenceDomain())
end

function CellQuadrature(trian::AppendedTriangulation,degree1::Integer,degree2::Integer)
  quad1 = CellQuadrature(trian.a,degree1)
  quad2 = CellQuadrature(trian.b,degree2)
  lazy_append(quad1,quad2,trian)
end

function CellQuadrature(trian::AppendedTriangulation,degree::Integer)
  CellQuadrature(trian,degree,degree)
end

function lazy_append(
  quad1::CellQuadrature,
  quad2::CellQuadrature,
  trian::AppendedTriangulation=lazy_append(quad1.trian,quad2.trian))

  @notimplementedif DomainStyle(quad1) != DomainStyle(quad2)
  cell_quad = lazy_append(quad1.cell_quad,quad2.cell_quad)
  cell_point = lazy_append(quad1.cell_point,quad2.cell_point)
  cell_weight = lazy_append(quad1.cell_weight,quad2.cell_weight)
  CellQuadrature(cell_quad,cell_point,cell_weight,trian,DomainStyle(quad1))
end

get_data(f::CellQuadrature) = f.cell_quad
get_triangulation(f::CellQuadrature) = f.trian
DomainStyle(::Type{CellQuadrature{DS}}) where DS = DS()

function change_domain(a::CellQuadrature,::ReferenceDomain,::PhysicalDomain)
  @notimplemented
end

function change_domain(a::CellQuadrature,::PhysicalDomain,::ReferenceDomain)
  @notimplemented
end

function get_cell_points(a::CellQuadrature)
  CellPoint(a.cell_point,a.trian,a.domain_style)
end

function integrate(f::CellField,quad::CellQuadrature)

  trian_f = get_triangulation(f)
  trian_x = get_triangulation(quad)

  if have_compatible_domains(trian_f,trian_x)
    nothing
  elseif have_compatible_domains(trian_f,get_background_triangulation(trian_x))
    nothing
  elseif have_compatible_domains(trian_x,get_background_triangulation(trian_f))
    @unreachable """\n
    CellField objects defined on a sub-triangulation cannot be integrated
    with a CellQuadrature defined on the underlying background mesh.

    This happens e.g. when trying to integrate a CellField defined on a Neumann boundary
    with a CellQuadrature defined on the underlying background mesh.
    """
  else
    @unreachable """\n
    Your are trying to integrate a CellField using a CellQuadrature defined on incompatible
    triangulations. Verify that either the two objects are defined in the same triangulation
    or that the triangulaiton of the CellField is the background triangulation of the CellQuadrature.
    """
  end

  b = change_domain(f,quad.trian,quad.domain_style)
  x = get_cell_points(quad)
  bx = b(x)
  if quad.domain_style == PhysicalDomain()
    lazy_map(IntegrationMap(),bx,quad.cell_weight)
  else
    cell_map = get_cell_map(quad.trian)
    cell_Jt = lazy_map(∇,cell_map)
    cell_Jtx = lazy_map(evaluate,cell_Jt,quad.cell_point)
    lazy_map(IntegrationMap(),bx,quad.cell_weight,cell_Jtx)
  end
end

function integrate(a,quad::CellQuadrature)
  b = CellField(a,quad.trian,quad.domain_style)
  integrate(b,quad)
end

# Some syntactic sugar

struct Integrand
  object
end

const ∫ = Integrand

(*)(a::Integrand,b::CellQuadrature) = integrate(a.object,b)
(*)(b::CellQuadrature,a::Integrand) = integrate(a.object,b)

# Cell measure
"""
Contributions added to the cells of the background Triangulation.
"""
function get_cell_measure(trian::Triangulation)
  quad = CellQuadrature(trian,0)
  cell_to_dV = integrate(1,quad)
  cell_to_bgcell = get_cell_to_bgcell(trian)
  bgtrian = get_background_triangulation(trian)
  bgcell_to_dV = zeros(num_cells(bgtrian))
  _meas_K_fill!(bgcell_to_dV,cell_to_dV,cell_to_bgcell)
  bgcell_to_dV
end

function _meas_K_fill!(bgcell_to_dV,cell_to_dV,cell_to_bgcell)
  cache = array_cache(cell_to_dV)
  for (cell, bgcell) in enumerate(cell_to_bgcell)
    dV = getindex!(cache,cell_to_dV,cell)
    bgcell_to_dV[bgcell] += dV
  end
end


text:
using SeisPlot, SeisAcoustic, FFTW, LinearAlgebra

# velocity and density model
nz = 101; nx = 301
vel = 3000 * ones(nz, nx);  # m/s
vel[51:end,:] .= 3500;  # m/s
rho = 2000 * ones(nz, nx);  # kg/m^3

# top boundary condition
free_surface = false    #(pml or free_surface)

# grid size
h = 10

# organize these parameters into a structure
fparams = FdParams(rho, vel, free_surface, h; data_format=Float64);

# generate observations of one frequency slice by a point source
isz = 5; isx = 151; dt = 0.001;
src = Source(isz, isx, dt, fparams; amp=1/rho[isz,isx]/vel[isz,isx]^2*10^5);

# total modelling length
tmax = 2.0;
nt   = floor(Int64, tmax/dt) + 1

# working frequency
df = 1.0 / (dt * nt);
f_work = 20.0;
iw = floor(Int64, f_work/df) + 1;
f_work = (iw-1) * df;
omega  = 2.0 * pi * f_work;
wlet   = zeros(nt)
copyto!(wlet, src.it_min, src.p)
fwave  = fft(wlet)

# LU decomposition of Helmhotz operator
H = get_helmholtz_LU(fparams, omega);

# source vector
s = zeros(Complex{Float64}, fparams.Nz*fparams.Nx)
s[src.src2spt] = fwave[iw]

# produce modelling true wave field
u = zeros(Complex{Float64}, fparams.Nz * fparams.Nx);
ldiv!(u, H, s);
SeisPlotTX(real(reshape(u, fparams.Nz, fparams.Nx)), wbox=9, hbox=3)

# specify the monochromatic recordings
irx = collect(1:2:nx);
irz = 2 * ones(Int64, length(irx));
dobs= MonochromaticRecordings(irz, irx, omega, fparams);
sample_spt2rec!(dobs, u);

# initial homogenous velocity model
vel = 3000 * ones(nz, nx)
fparams = FdParams(rho, vel, free_surface, h; data_format=Float64);

# helmholtz operator of homogeneous model
H = get_helmholtz_LU(fparams, omega);

# source side wavefield
u = zeros(Complex{Float64}, fparams.Nz * fparams.Nx);
ldiv!(u, H, s);
SeisPlotTX(real(reshape(u, fparams.Nz, fparams.Nx)), wbox=9, hbox=3)

dsyn = MonochromaticRecordings(irz, irx, omega, fparams);
sample_spt2rec!(dsyn, u);

# compute the residue
dres = get_residue(dsyn, dobs);
g = velocity_gradient(dres, u, H, fparams);
g = reshape(g, fparams.nz, fparams.nx)
SeisPlotTX(real(g), wbox=9, hbox=3)

# compute the gradient numerically
vel0    = 3000 * ones(nz, nx);
delta_m = 1e-7;

# one model parameter's gradient
iz = 37; ix = 137;

# positive velocity model partubation
vel0[iz,ix] = vel0[iz,ix] + delta_m

# model parameter
fparams = FdParams(rho, vel0, free_surface, h; data_format=Float64);
H = get_helmholtz_LU(fparams, omega);

# source side wavefield
u = zeros(Complex{Float64}, fparams.Nz * fparams.Nx);
ldiv!(u, H, s);

d1 = MonochromaticRecordings(irz, irx, omega, fparams);
sample_spt2rec!(d1, u);
SeisPlotTX(real(reshape(u, fparams.Nz, fparams.Nx)), wbox=9, hbox=3)

# compute the residue
dres = get_residue(d1, dobs);
J_plus = 1/2.0 * (norm(dres.p))^2

# negative velocity model partubation
vel0[iz,ix] = vel0[iz,ix] - 2*delta_m

fparams = FdParams(rho, vel0, free_surface, h; data_format=Float64);
H = get_helmholtz_LU(fparams, omega);

# source side wavefield
u = zeros(Complex{Float64}, fparams.Nz * fparams.Nx);
ldiv!(u, H, s);

d2 = MonochromaticRecordings(irz, irx, omega, fparams);
sample_spt2rec!(d2, u);

# compute the residue
dres = get_residue(d2, dobs);
J_minus = 1/2.0 * (norm(dres.p))^2;

# compute gradient numerically
g_num = (J_plus - J_minus) / (2 * delta_m)
g[iz,ix]
(g[iz,ix] - g_num) / g_num


text:
const MANTISSA_MASK = Base.significand_mask(Float64)
const EXP_MASK = Base.exponent_mask(Float64) >> Base.significand_bits(Float64)

memcpy(d, doff, s, soff, n) = ccall(:memcpy, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Int), d + doff - 1, s + soff - 1, n)
memmove(d, doff, s, soff, n) = ccall(:memmove, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Int), d + doff - 1, s + soff - 1, n)

pow5_bitcount(::Type{Float16}) = 30
pow5_bitcount(::Type{Float32}) = 61
pow5_bitcount(::Type{Float64}) = 121

pow5_inv_bitcount(::Type{Float16}) = 30
pow5_inv_bitcount(::Type{Float32}) = 59
pow5_inv_bitcount(::Type{Float64}) = 122

qinvbound(::Type{Float16}) = 4
qinvbound(::Type{Float32}) = 9
qinvbound(::Type{Float64}) = 21

qbound(::Type{Float16}) = 15
qbound(::Type{Float32}) = 31
qbound(::Type{Float64}) = 63

"""
    Ryu.log10pow2(e::Integer)

Computes `floor(log10(2^e))`. This is valid for all `e < 1651`.
"""
log10pow2(e) = (e * 78913) >> 18

"""
    Ryu.log10pow5(e::Integer)

Computes `floor(log10(5^e))`. This is valid for all `e < 2621`.
"""
log10pow5(e) = (e * 732923) >> 20

"""
    Ryu.pow5bits(e)

Computes `e == 0 ? 1 : ceil(log2(5^e))`. This is valid for `e < 3529` (if performend in `Int32` arithmetic).
"""
pow5bits(e) = ((e * 1217359) >> 19) + 1

""""
     Ryu.mulshift(m::U, mula, j) where {U<:Unsigned}

Compute `(m * mul) >> j`, where `j >= 8*sizeof(U)`. The type of the results is the larger of `U` or `UInt32`.
"""
@inline function mulshift(m::U, mul, j) where {U<:Unsigned}
    W = widen(U)
    nbits = 8*sizeof(U)
    return ((((W(m) * (mul % U)) >> nbits) + W(m) * (mul >> nbits)) >> (j - nbits)) % promote_type(U,UInt32)
end

indexforexp(e) = div(e + 15, 16)
pow10bitsforindex(idx) = 16 * idx + 120
lengthforindex(idx) = div(((Int64(16 * idx) * 1292913986) >> 32) + 1 + 16 + 8, 9)

"""
    Ryu.pow5(x, p)

Return `true` if `5^p` is a divisor of `x`.
"""
@inline function pow5(x, p)
    count = 0
    while true
        q = div(x, 5)
        r = x - 5 * q
        r != 0 && return count >= p
        x = q
        count += 1
    end
end

"""
    Ryu.pow2(x, p)

Return `true` if `2^p` is a divisor of `x`. In other words, if the trailing `p` bits of `x` are zero.
"""
pow2(x, p) = (x & ((Int64(1) << p) - 1)) == 0

"""
    Ryu.decimallength(v)

The number of decimal digits of the integer `v`.
"""
@inline function decimallength(v)
    v >= 10000000000000000 && return 17
    v >= 1000000000000000 && return 16
    v >= 100000000000000 && return 15
    v >= 10000000000000 && return 14
    v >= 1000000000000 && return 13
    v >= 100000000000 && return 12
    v >= 10000000000 && return 11
    v >= 1000000000 && return 10
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end
@inline function decimallength(v::UInt32)
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end
@inline function decimallength(v::UInt16)
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end

@inline function mulshiftinvsplit(::Type{Float64}, mv, mp, mm, i, j)
    @inbounds mul = DOUBLE_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftinvsplit(::Type{Float32}, mv, mp, mm, i, j)
    @inbounds mul = FLOAT_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftinvsplit(::Type{Float16}, mv, mp, mm, i, j)
    @inbounds mul = HALF_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float64}, mv, mp, mm, i, j)
    @inbounds mul = DOUBLE_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float32}, mv, mp, mm, i, j)
    @inbounds mul = FLOAT_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float16}, mv, mp, mm, i, j)
    @inbounds mul = HALF_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

"""
    Ryu.umul256(a::UInt128, bHi::UInt64, bLo::UInt64)::Tuple{UInt128, UInt128}

Compute `p = a*b` where `b = bLo + bHi<<64`, returning the result as `pLo, pHi` where `p = pLo + pHi<<128`.
"""
@inline function umul256(a, bHi, bLo)
    aLo = a % UInt64
    aHi = (a >> 64) % UInt64

    b00 = UInt128(aLo) * bLo
    b01 = UInt128(aLo) * bHi
    b10 = UInt128(aHi) * bLo
    b11 = UInt128(aHi) * bHi

    b00Lo = b00 % UInt64
    b00Hi = (b00 >> 64) % UInt64

    mid1 = b10 + b00Hi
    mid1Lo = mid1 % UInt64
    mid1Hi = (mid1 >> 64) % UInt64

    mid2 = b01 + mid1Lo
    mid2Lo = mid2 % UInt64
    mid2Hi = (mid2 >> 64) % UInt64

    pHi = b11 + mid1Hi + mid2Hi
    pLo = (UInt128(mid2Lo) << 64) | b00Lo
    return pLo, pHi
end

"""
    Ryu.umul256_hi(a::UInt128, bHi::UInt64, bLo::UInt64)::UInt128

Compute `pHi = (a*b)>>128` where `b = bLo + bHi<<64`.
"""
@inline umul256_hi(a, bHi, bLo) = umul256(a, bHi, bLo)[2]

"""
    Ryu.mulshiftmod1e9(m, mula, mulb, mulc, j)::UInt32

Compute `(m * mul) >> j % 10^9` where `mul = mula + mulb<<64 + mulc<<128`, and `j >= 128`.
"""
@inline function mulshiftmod1e9(m, mula, mulb, mulc, j)
    b0 = UInt128(m) * mula
    b1 = UInt128(m) * mulb
    b2 = UInt128(m) * mulc
    mid = b1 + ((b0 >> 64) % UInt64)
    s1 = b2 + ((mid >> 64) % UInt64)
    v = s1 >> (j - 128)
    multiplied = umul256_hi(v, 0x89705F4136B4A597, 0x31680A88F8953031)
    shifted = (multiplied >> 29) % UInt32
    return (v % UInt32) - UInt32(1000000000) * shifted
end

@inline function append_n_digits(olength, digits, buf, pos)
    i = 0
    while digits >= 10000
        c = digits % 10000
        digits = div(digits, 10000)
        c0 = (c % 100) << 1
        c1 = div(c, 100) << 1
        unsafe_copyto!(buf, pos + olength - i - 2, DIGIT_TABLE, c0 + 1, 2)
        unsafe_copyto!(buf, pos + olength - i - 4, DIGIT_TABLE, c1 + 1, 2)
        i += 4
    end
    if digits >= 100
        c = (digits % 100) << 1
        digits = div(digits, 100)
        unsafe_copyto!(buf, pos + olength - i - 2, DIGIT_TABLE, c + 1, 2)
        i += 2
    end
    if digits >= 10
        c = digits << 1
        unsafe_copyto!(buf, pos + olength - i - 2, DIGIT_TABLE, c + 1, 2)
        i += 2
    else
        buf[pos] = UInt8('0') + digits
        i += 1
    end
    return pos + i
end

@inline function append_d_digits(olength, digits, buf, pos, decchar)
    i = 0
    while digits >= 10000
        c = digits % 10000
        digits = div(digits, 10000)
        c0 = (c % 100) << 1
        c1 = div(c, 100) << 1
        unsafe_copyto!(buf, pos + olength + 1 - i - 2, DIGIT_TABLE, c0 + 1, 2)
        unsafe_copyto!(buf, pos + olength + 1 - i - 4, DIGIT_TABLE, c1 + 1, 2)
        i += 4
    end
    if digits >= 100
        c = (digits % 100) << 1
        digits = div(digits, 100)
        unsafe_copyto!(buf, pos + olength + 1 - i - 2, DIGIT_TABLE, c + 1, 2)
        i += 2
    end
    if digits >= 10
        c = digits << 1
        buf[pos] = DIGIT_TABLE[c + 1]
        buf[pos + 1] = decchar
        buf[pos + 2] = DIGIT_TABLE[c + 2]
        i += 3
    else
        buf[pos] = UInt8('0') + digits
        buf[pos + 1] = decchar
        i += 2
    end
    return pos + i
end

@inline function append_c_digits(count, digits, buf, pos)
    i = 0
    while i < count - 1
        c = (digits % 100) << 1
        digits = div(digits, 100)
        unsafe_copyto!(buf, pos + count - i - 2, DIGIT_TABLE, c + 1, 2)
        i += 2
    end
    if i < count
        buf[pos + count - i - 1] = UInt8('0') + (digits % 10)
        i += 1
    end
    return pos + i
end

@inline function append_nine_digits(digits, buf, pos)
    if digits == 0
        for _ = 1:9
            buf[pos] = UInt8('0')
            pos += 1
        end
        return pos
    end
    i = 0
    while i < 5
        c = digits % 10000
        digits = div(digits, 10000)
        c0 = (c % 100) << 1
        c1 = div(c, 100) << 1
        unsafe_copyto!(buf, pos + 7 - i, DIGIT_TABLE, c0 + 1, 2)
        unsafe_copyto!(buf, pos + 5 - i, DIGIT_TABLE, c1 + 1, 2)
        i += 4
    end
    buf[pos] = UInt8('0') + digits
    i += 1
    return pos + i
end

const BIG_MASK = (big(1) << 64) - 1

const POW10_SPLIT = collect(Iterators.flatten(map(0:63) do idx
    pow10bits = pow10bitsforindex(idx)
    map(0:lengthforindex(idx)-1) do i
        v = (div(big(1) << pow10bits, big(10)^(9 * i)) + 1) % ((big(10)^9) << 136)
        return (UInt64(v & BIG_MASK), UInt64((v >> 64) & BIG_MASK), UInt64((v >> 128) & BIG_MASK))
    end
end))

function generateinversetables()
    POW10_OFFSET_2 = Vector{UInt16}(undef, 68 + 1)
    MIN_BLOCK_2 = fill(0xff, 68 + 1)
    POW10_SPLIT_2 = Tuple{UInt64, UInt64, UInt64}[]
    lowerCutoff = big(1) << (54 + 8)
    for idx = 0:67
        POW10_OFFSET_2[idx + 1] = length(POW10_SPLIT_2)
        i = 0
        while true
            v = ((big(10)^(9 * (i + 1)) >> (-(120 - 16 * idx))) % (big(10)^9) << (120 + 16))
            if MIN_BLOCK_2[idx + 1] == 0xff && ((v * lowerCutoff) >> 128) == 0
                i += 1
                continue
            end
            if MIN_BLOCK_2[idx + 1] == 0xff
                MIN_BLOCK_2[idx + 1] = i
            end
            v == 0 && break
            push!(POW10_SPLIT_2, ((v & BIG_MASK) % UInt64, ((v >> 64) & BIG_MASK) % UInt64, ((v >> 128) & BIG_MASK) % UInt64))
            i += 1
        end
    end
    POW10_OFFSET_2[end] = length(POW10_SPLIT_2)
    MIN_BLOCK_2[end] = 0x00

    return POW10_OFFSET_2, MIN_BLOCK_2, POW10_SPLIT_2
end

const POW10_OFFSET_2, MIN_BLOCK_2, POW10_SPLIT_2 = generateinversetables()

"""
    Ryu.pow5invsplit(T, i)

Compute `floor(2^k/5^i)+1`, where `k = pow5bits(i) - 1 + pow5_inv_bitcount(T)`. The result
is an unsigned integer twice as wide as `T` (i.e. a `UInt128` if `T == Float64`), with
`pow5_inv_bitcount(T)` significant bits.
"""
function pow5invsplit(::Type{T}, i) where {T<:AbstractFloat}
    W = widen(uinttype(T))
    pow = big(5)^i
    inv = div(big(1) << (ndigits(pow, base=2) - 1 + pow5_inv_bitcount(T)), pow) + 1
    return W(inv)
end

"""
    Ryu.pow5split(T, i)

Compute `floor(5^i/2^k)`, where `k = pow5bits(i) - pow5_bitcount(T)`. The result is an
unsigned integer twice as wide as `T` (i.e. a `UInt128` if `T == Float64`), with
`pow5_bitcount(T)` significant bits.
"""
function pow5split(::Type{T}, i) where {T<:AbstractFloat}
    W = widen(uinttype(T))
    pow = big(5)^i
    return W(pow >> (ndigits(pow, base=2) - pow5_bitcount(T)))
end

const DOUBLE_POW5_INV_SPLIT = map(i->pow5invsplit(Float64, i), 0:291)
const FLOAT_POW5_INV_SPLIT = map(i->pow5invsplit(Float32, i), 0:30)
const HALF_POW5_INV_SPLIT = map(i->pow5invsplit(Float16, i), 0:17)

const DOUBLE_POW5_SPLIT = map(i->pow5split(Float64, i), 0:325)
const FLOAT_POW5_SPLIT = map(i->pow5split(Float32, i), 0:46)
const HALF_POW5_SPLIT = map(i->pow5split(Float16, i), 0:23)

const DIGIT_TABLE = UInt8[
  '0','0','0','1','0','2','0','3','0','4','0','5','0','6','0','7','0','8','0','9',
  '1','0','1','1','1','2','1','3','1','4','1','5','1','6','1','7','1','8','1','9',
  '2','0','2','1','2','2','2','3','2','4','2','5','2','6','2','7','2','8','2','9',
  '3','0','3','1','3','2','3','3','3','4','3','5','3','6','3','7','3','8','3','9',
  '4','0','4','1','4','2','4','3','4','4','4','5','4','6','4','7','4','8','4','9',
  '5','0','5','1','5','2','5','3','5','4','5','5','5','6','5','7','5','8','5','9',
  '6','0','6','1','6','2','6','3','6','4','6','5','6','6','6','7','6','8','6','9',
  '7','0','7','1','7','2','7','3','7','4','7','5','7','6','7','7','7','8','7','9',
  '8','0','8','1','8','2','8','3','8','4','8','5','8','6','8','7','8','8','8','9',
  '9','0','9','1','9','2','9','3','9','4','9','5','9','6','9','7','9','8','9','9'
]

const POW10_OFFSET = UInt16[
  0, 2, 5, 8, 12, 16, 21, 26, 32, 39,
  46, 54, 62, 71, 80, 90, 100, 111, 122, 134,
  146, 159, 173, 187, 202, 217, 233, 249, 266, 283,
  301, 319, 338, 357, 377, 397, 418, 440, 462, 485,
  508, 532, 556, 581, 606, 632, 658, 685, 712, 740,
  769, 798, 828, 858, 889, 920, 952, 984, 1017, 1050,
  1084, 1118, 1153, 1188
]

max_stars_repo_name: taqtiqa-mark/julia

max_issues_repo_name: taqtiqa-mark/julia

max_forks_repo_name: taqtiqa-mark/julia


text:
using Liblet

# Create an Automaton from a given grammar

A = Automaton(Grammar("""
        S -> a A
        S -> a B
        A -> b B
        A -> b C
        B -> c A
        B -> c C
        C -> a
    """))


text:
using LinearAlgebra
using ADCME

n = 101 # number of grid nodes in [0,1]
h = 1/(n-1)
x = LinRange(0,1,n)[2:end-1]

b = Variable(10.0) # create a Variable for `b`
A = diagm(0=>2/h^2*ones(n-2), -1=>-1/h^2*ones(n-3), 1=>-1/h^2*ones(n-3)) # discrete Laplacian matrix
B = b*A + diagm(0=>ones(n-2))  # coefficient matrix
f = @. 4*(2 + x - x^2) # right hand side
u = B\f # solve the equation
ue = u[div(n+1,2)] # extract values at x=0.5

loss = (ue-1.0)^2 # form the loss function

# Optimization
opt = ScipyOptimizerInterface(loss)
sess = Session(); init(sess)
ScipyOptimizerMinimize(sess, opt)

println("Estimated b = ", run(sess, b))


text:
# PBE correlation (without LDA part)
# iflag=1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
function XC_c_pbe( rho, grho )

    ga = 0.0310906908696548950
    be = 0.06672455060314922
    third = 1.0/3.0
    pi34 = 0.6203504908994
    xkf = 1.919158292677513
    xks = 1.128379167095513

    rs = pi34/rho^third
    ec, vc = XC_c_pw( rho )

    kf = xkf/rs
    ks = xks * sqrt(kf)
    t = sqrt(grho) / (2.0 * ks * rho)

    expe = exp(-ec/ga)
    af = be / ga * (1.0 / (expe - 1.0) )

    bf = expe * (vc - ec)

    y = af * t * t

    xy = (1.0 + y) / (1.0 + y + y * y)

    qy = y * y * (2.0 + y) / (1.0 + y + y * y)^2

    s1 = 1.0 + be / ga * t * t * xy
    h0 = ga * log(s1)

    dh0 = be * t * t / s1 * ( -7.0 / 3.0 * xy - qy * (af * bf / be - 7.0 / 3.0) )

    ddh0 = be/(2.0 * ks * ks * rho) * (xy - qy) / s1

    sc = rho * h0
    v1c = h0 + dh0
    v2c = ddh0

    return sc, v1c, v2c

end


text:
function expand_uses( sys )  #=

expands the <use> ... </use> sections by substituting their component definitions and then returns the result.

This function requires that the Julia package EzXML has been loaded prior to the execution.

Note! At this moment, variables in
  - multiphase elements:
    lhs and rhs pins, matrix parameters
  - control blocks (search for elseif elm.name == "ctrl"):
    input and output variables, one-line Julia code
are not replaced according to <use> substitutions.

Author: Taku Noda
Started on Nov. 5, 2020
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

cmpnnts = Component[]  # all components will be stored.
subsyss = Subsystem[]  # all subsystems will be stored.

sys_root = sys.root  # getting the root element.

# reading the component and subsystem definitions:
num_cmpnnts = 0
num_subsyss = 0
for elm in eachelement( sys_root )
  tag = elm.name
  if tag == "component"
    name = elm[ "name" ]
    type = elm[ "type" ]
    calc = ""
    for elm1 in eachelement( elm )
      if elm1.name == "calc"
        content = strip( elm1.content, [ ' ', '\t', '\n' ] )
        if !isempty( content )
          if content[ end ] == ';'
            calc = content
          else
            calc = content*";"
          end
        end
        break
      end
    end
    cmpnnt = Component( name, type, calc, elm )
    push!( cmpnnts, cmpnnt )
    num_cmpnnts = num_cmpnnts + 1
  elseif tag == "subsystem"
    name = elm[ "name" ]
    type = elm[ "type" ]
    calc = ""
    for elm1 in eachelement( elm )
      if elm1.name == "calc"
        content = strip( elm1.content, [ ' ', '\t', '\n' ] )
        if !isempty( content )
          if content[ end ] == ';'
            calc = content
          else
            calc = content*";"
          end
        end
        break
      end
    end
    subsys = Subsystem( name, type, calc, elm )
    push!( subsyss, subsys )
    num_subsyss = num_subsyss + 1
  end
end

uses = Use[]  # all use statements will be stored.

# collecting all use statements:
num_uses = 0
# from components:
for itr_cmpnnts = 1:num_cmpnnts
  cmpnnt_root = ( cmpnnts[ itr_cmpnnts ] ).root
  for elm in eachelement( cmpnnt_root )
    if elm.name == "use"
      path = nodepath( elm )
      if itr_cmpnnts == 1
        isdf_error( "A component must be defined before it is used, and thus the first component cannot have a <use> ... </use> section.", path )
      end
      component = elm[ "component" ]
      cmpnnt_id = 0
      for jtr_cmpnnts = 1:itr_cmpnnts
        if jtr_cmpnnts == itr_cmpnnts
          isdf_error( "The definition of the component \"$component\" cannot be found before this use statement. A component must be defined before it is used. Or, a recursive use of a component is not allowed.", path )
        end
        if component == ( cmpnnts[ jtr_cmpnnts ] ).name
          cmpnnt_id = jtr_cmpnnts
          break
        end
      end
      as = elm[ "as" ]
      subst = ""
      content = strip( elm.content, [ ' ', '\t', '\n' ] )
      if !isempty( content )
        if content[ end ] == ';'
          subst = content
        else
          subst = content*";"
        end
      end
      use = Use( cmpnnt_id, as, elm, subst )
      push!( uses, use )
      num_uses = num_uses + 1
    end
  end
end
# from subsystems:
for itr_subsyss = 1:num_subsyss
  subsys_root = ( subsyss[ itr_subsyss ] ).root
  for elm in eachelement( subsys_root )
    if elm.name == "use"
      path = nodepath( elm )
      component = elm[ "component" ]
      cmpnnt_id = 0
      for itr_cmpnnts = 1:num_cmpnnts
        if component == ( cmpnnts[ itr_cmpnnts ] ).name
          cmpnnt_id = itr_cmpnnts
          break
        end
      end
      if cmpnnt_id == 0
        isdf_error( "The definition of the component \"$component\" cannot be found before this use statement. A component must be defined before it is used.", path )
      end
      as = elm[ "as" ]
      subst = ""
      content = strip( elm.content, [ ' ', '\t', '\n' ] )
      if !isempty( content )
        if content[ end ] == ';'
          subst = content
        else
          subst = content*";"
        end
      end
      use = Use( cmpnnt_id, as, elm, subst )
      push!( uses, use )
      num_uses = num_uses + 1
    end
  end
end

# expanding the use statements:
for itr_uses = 1:num_uses
  use = uses[ itr_uses ]
  use_elm = use.use_elm
  use_as = use.as
  cmpnnt_root = ( cmpnnts[ use.cmpnnt_id ] ).root
  io = IOBuffer();
  println( io, cmpnnt_root )
  cmpnnt1_doc = parsexml( String( take!( io ) ) )
  cmpnnt1_root = cmpnnt1_doc.root
  # appending use substitutions on top of calc in cmpnnt1_root:
  hascalc = false
  for elm in eachelement( cmpnnt1_root )
    if elm.name == "calc"
      elm.content = ( use.subst )*( elm.content )
      hascalc = true
    end
  end
  elms = elements( cmpnnt1_root )
  if !hascalc
    calc = ElementNode( "calc" )
    calc.content = use.subst
    linkprev!( elms[ 1 ], calc )
  end
  #
  curr_elm = use_elm
  num_elms = length( elms )
  for itr_elms = 1:num_elms
    elm = elms[ itr_elms ]
    unlink!( elm )
    if elm.name != "calc"
      elm[ "name" ] = use_as*'/'*( elm[ "name" ] )
      lhs = elm[ "lhs" ]
      if lhs[ 1 ] != '_'
        elm[ "lhs" ] = use_as*'/'*( elm[ "lhs" ] )
      end
      rhs = elm[ "rhs" ]
      if rhs[ 1 ] != '_'
        elm[ "rhs" ] = use_as*'/'*( elm[ "rhs" ] )
      end
    end
    curr_elm = linknext!( curr_elm, elm )
  end
  unlink!( use_elm )
end

# unlinking used components:
for itr_cmpnnts = 1:num_cmpnnts
  cmpnnt_root = ( cmpnnts[ itr_cmpnnts ] ).root
  unlink!( cmpnnt_root )
end

dq = Char( 34 )  # double quotation.
bsdq = Char( 92 )*dq  # backslash + double quotation.

for itr_subsyss = 1:num_subsyss
  subsys_root = ( subsyss[ itr_subsyss ] ).root
  calc = ""
  for elm in eachelement( subsys_root )
    path = nodepath( elm )
    if elm.name == "calc"
      calc = elm.content
    elseif elm.name == "branch"
      if haskey( elm, "val" )
        calc1 = calc*elm[ "val" ]
        calc1 = replace( calc1, r" |\t|\n" => "" )
        calc1_p = try
          Meta.parse( calc1 )
        catch
          isdf_error( "An error occurred when the following calculation code was parsed: $calc1.", path )
        end
        val = try
          eval( calc1_p )
        catch
          isdf_error( "An error occurred when the following calculation was evaluated: $calc1.", path )
        end
        elm[ "val" ] = string( val )
      end
      if haskey( elm, "lhs" )
        lhs = elm[ "lhs" ]
        if lhs[ 1 ] == '_'
          calc1 = calc*lhs
          calc1 = replace( calc1, r" |\t|\n" => "" )
          calc1_p = try
            Meta.parse( calc1 )
          catch
            isdf_error( "An error occurred when the following calculation code was parsed: $calc1.", path )
          end
          lhs = try
            eval( calc1_p )
          catch
            isdf_error( "An error occurred when the following calculation was evaluated: $calc1.", path )
          end
          elm[ "lhs" ] = lhs
        end
      end
      if haskey( elm, "rhs" )
        rhs = elm[ "rhs" ]
        if rhs[ 1 ] == '_'
          calc1 = calc*rhs
          calc1 = replace( calc1, r" |\t|\n" => "" )
          calc1_p = try
            Meta.parse( calc1 )
          catch
            isdf_error( "An error occurred when the following calculation code was parsed: $calc1.", path )
          end
          rhs = try
            eval( calc1_p )
          catch
            isdf_error( "An error occurred when the following calculation was evaluated: $calc1.", path )
          end
          elm[ "rhs" ] = rhs
        end
      end
    # elseif elm.name == "ctrl"
    end
  end
  elms = elements( subsys_root )
  for itr_elms = 1:length( elms )
    elm = elms[ itr_elms ]
    if elm.name == "calc"
      unlink!( elm )
    end
  end
end

return sys

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
using DataFrames
using DiscreteValueIteration
using DrWatson

@testset "simrunning.jl" begin
    pp = simple_pp()
    trace = simple_short_trace(pp)
    traces = [simple_short_trace(pp), empty_trace(pp), simple_short_trace(pp)]
    N = length(traces)

    pp_params = Dict(:test=>1)
    name = "test_problem"
    data = @dict(traces, pp, pp_params, name)

    @test isa(PMDPs.flatrate(pp, traces, MersenneTwister(1)), DataFrame)

    PMDPs.process_data(data, PMDPs.flatrate)
    @test isfile(datadir("results", "test_problem", "flatrate_N=$(N)_test=1.bson"))

    PMDPs.process_data(data, PMDPs.hindsight)
    @test isfile(datadir("results", "test_problem", "hindsight_N=$(N)_test=1.bson"))

    # PMDPs.vi(pp, pp_params, traces, rnd)
    PMDPs.process_data(data, PMDPs.vi)
    @test isfile(datadir("results", "test_problem", "vi_N=$(N)_test=1.bson"))
    @test isfile(datadir("vi_policies", "test_problem", "vi_test=1_.bson" ))
    PMDPs.process_data(data, PMDPs.vi)

    PMDPs.process_data(data, PMDPs.mcts)
    @test isfile(datadir("results", "test_problem", "mcts_N=$(N)_test=1.bson"))

    rm(datadir("results", "test_problem"), recursive=true)
    rm(datadir("vi_policies", "test_problem"), recursive=true)
end


text:
include("UQ.jl")

"""
    SSampling(m::Integer, scale::Real=10,  N::Integer=100, N_samples::Integer=20, N_matrices::Integer=500; randomization="A")

    This function gather samples of (x_exact - x_approximate)^{T} \Sigma^{-1} (x_exact - x_approximate).
    The exact solution is sampled either from x ~ N(0, I) (`randomization`= anything but "A") or x ~ N(0, A^{-1}) (`randomization`="A").
    The matrix A = U D U^{T} drawn as follows. U ~ uniform over O(N), D_{ii} ~ exp distribution sith scale `s`
    `N` fixes dimension, `N_samples` -- number of x generated for each matrix, `N_matrices` -- number of matrices generated.

    (see the article, for more details)
    The output `Statistics` is an array of the size (5, `N_samples`*`N_matrices`).
    `Statistics[1, :]` -- samples, correspond to the mode of inverse gamma distribution (with additional observation) [compare with chi^2]
    `Statistics[2, :]` -- samples, correspond to the mode of inverse gamma distribution (without additional observation) [compare with chi^2]
    `Statistics[3, :]` -- samples, correspond to the hierarchical modelling (with additional observation) [compare with F]
    `Statistics[4, :]` -- samples, correspond to the hierarchical modelling (without additional observation) [compare with F]
    `Statistics[5, :]` -- \Sigma = A^{-1} - V (V^T A V)^-1 V^T, i.e. Bayesian CG [compare with chi^2]
"""
function SSampling(m::Integer, scale::Real=10,  N::Integer=100, N_samples::Integer=20, N_matrices::Integer=500; randomization="A")
    Statistics = zeros(5, N_samples*N_matrices)
    for j=1:N_matrices
        A, s, U = randS(N, scale)
        inv_A = U*Diagonal(s.^(-1))*transpose(U)
        L = U*Diagonal(s.^(-1/2))
        if randomization == "A"
            draw = (n=N, l=L) -> randx(zeros(n), l)
        else
            draw = (n=N) -> randx(zeros(n))
        end

        for i=1:N_samples
            x = draw()
            b = A*x
            V = Krylov(A, b, m)
            y = projection(zeros(N), A, b, V, V)
            sigma = inv_A - V*inv(transpose(V)*A*V)*transpose(V)
            eps = 0
            k = 0
            while k == 0
                try
                    k += 1
                    Statistics[5, i + (j-1)*N_samples] = transpose(x - y)*pinv(sigma + I*eps)*(x - y)
                catch LAPACKException
                    k -= 1
                    eps = 1e-6
                end
            end

            delta = inv(transpose(V)*A*V)*transpose(V)*b
            alpha_2 = m
            beta_2 = norm(delta)^2
            s_2 = beta_2/(alpha_2 + 2)

            z = draw()
            b = A*z
            V = Krylov(A, b, m)
            P = I - V*inv(transpose(V)*A*V)*transpose(V)*A
            sample = P*z
            alpha_1 = (N-m)
            beta_1 = norm(sample)^2
            s_1 = beta_1/(alpha_1 + 2)

            squared_error = norm(x - y)^2
            Statistics[1, i + (j-1)*N_samples] = squared_error/s_1
            Statistics[2, i + (j-1)*N_samples] = squared_error/s_2
            Statistics[3, i + (j-1)*N_samples] = (alpha_1/beta_1)*squared_error/(N-m)
            Statistics[4, i + (j-1)*N_samples] = (alpha_2/beta_2)*squared_error/(N-m)
        end
    end
    return Statistics
end

"""
    ASampling(m::Integer, scale::Real=10,  N::Integer=100, N_samples::Integer=20, N_matrices::Integer=500; randomization="A")

    This function gather samples of (x_exact - x_approximate)^{T} \Sigma^{-1} (x_exact - x_approximate).
    The exact solution is sampled either from x ~ N(0, I) (`randomization`= anything but "A") or b ~ N(0, I) (`randomization`="A").
    The matrix A = U D V^{T} drawn as follows. U ~ uniform over O(N), V ~ uniform over O(N), D_{ii} ~ exp distribution sith scale `s`
    `N` fixes dimension, `N_samples` -- number of x generated for each matrix, `N_matrices` -- number of matrices generated.

    (see the article, for more details)
    The output `Statistics` is an array of the size (5, `N_samples`*`N_matrices`).
    `Statistics[1, :]` -- samples, correspond to the mode of inverse gamma distribution (with additional observation) [compare with chi^2]
    `Statistics[2, :]` -- samples, correspond to the mode of inverse gamma distribution (without additional observation) [compare with chi^2]
    `Statistics[3, :]` -- samples, correspond to the hierarchical modelling (with additional observation) [compare with F]
    `Statistics[4, :]` -- samples, correspond to the hierarchical modelling (without additional observation) [compare with F]
    `Statistics[5, :]` -- \Sigma = A^{-1} - V (V^T A V)^-1 V^T, i.e. Bayesian CG [compare with chi^2]
"""
function ASampling(m::Integer, scale::Real=10,  N::Integer=100, N_samples::Integer=20, N_matrices::Integer=500; randomization="A", symmetric=true)
    Statistics = zeros(5, N_samples*N_matrices)
    for j=1:N_matrices
        if symmetric
            A, s, U = randS(N, scale)
            L = U*Diagonal(s.^(-1))*transpose(U)
            inv_A = L*transpose(L)
        else
            A, s, U, V = randA(N, scale)
            L = V*Diagonal(s.^(-1))*transpose(U)
            inv_A = L*transpose(L)
        end

        if randomization == "A"
            draw = (n=N, l=L) -> randx(zeros(n), l)
        else
            draw = (n=N) -> randx(zeros(n))
        end

        for i=1:N_samples
            x = draw()
            b = A*x
            V = Krylov(A, b, m)
            y = projection(zeros(N), A, b, V, A*V)
            sigma = inv_A - V*inv(transpose(A*V)*A*V)*transpose(V)
            eps = 0
            k = 0
            while k == 0
                try
                    k += 1
                    Statistics[5, i + (j-1)*N_samples] = transpose(x - y)*pinv(sigma + I*eps)*(x - y)
                catch LAPACKException
                    k -= 1
                    eps = 1e-6
                end
            end

            delta = inv(transpose(A*V)*A*V)*transpose(A*V)*b
            alpha_2 = m
            beta_2 = norm(delta)^2
            s_2 = beta_2/(alpha_2 + 2)

            z = draw()
            b = A*z
            V = Krylov(A, b, m)
            P = I - V*inv(transpose(A*V)*A*V)*transpose(A*V)*A
            sample = P*z
            beta_1 = norm(sample)^2
            alpha_1 = (N-m)
            s_1 = beta_1/(alpha_1 + 2)

            squared_error = norm(x - y)^2
            Statistics[1, i + (j-1)*N_samples] = squared_error/s_1
            Statistics[2, i + (j-1)*N_samples] = squared_error/s_2
            Statistics[3, i + (j-1)*N_samples] = (alpha_1/beta_1)*squared_error/(N-m)
            Statistics[4, i + (j-1)*N_samples] = (alpha_2/beta_2)*squared_error/(N-m)
        end
    end
    return Statistics
end

max_stars_repo_name: VLSF/BayesKrylov

max_issues_repo_name: VLSF/BayesKrylov

max_forks_repo_name: VLSF/BayesKrylov


text:
module HexGreedImages

using Images, Reel, ImageDraw

using ...PCG.Geometry
using ...PCG.Geometry: Point
using ...PCG.Types
using ...PCG.Storages
using ...PCG.Storages.LinearStorages
using ...PCG.Topologies
using ...PCG.Topologies.HexGreedTopologies
using ...PCG.Universes
using ...PCG.Universes: Universe
using ...PCG.Operations
using ...PCG.Recorders.GreedImages

export HexSprite

# TODO: optimize algorithm
function recursive_fill!(image, point, color)

    try
        if image[point.y, point.x] == color
            return
        end
    catch e
        return
    end

    image[point.y, point.x] = color

    for delta in [Point(0, 1), Point(0, -1), Point(1, 0), Point(-1, 0)]
        recursive_fill!(image, point + delta, color)
    end

end

function HexSprite(color, cell_size::Size)
    image_size = ceil(Int64, CELL_SIZE * cell_size) + Size(2, 2)

    image = fill(RGBA(0, 0, 0, 0), yx(image_size))

    center = Point(image_size / 2)

    cell = HexGreedIndex(0, 0, 0)

    polygon = Polygon([xy(round(Int64, center + point * cell_size))
                       for point in cell_corners(cell)])

    draw!(image, polygon, color)

    recursive_fill!(image, round(Int64, center), color)

    return Sprite(color, image)
end

function GreedImages.node_position(recorder::GreedImageRecorder, index::HexGreedIndex, canvas_size::Size)
    point = Point(canvas_size / 2) + (cell_center(index) - Point(1.0, 1.0)) * recorder.cell_size
    return round(Int64, point)
end

function GreedImages.canvas_size(topology::HexGreedTopology, recorder::GreedImageRecorder)::Size
    return ceil(Int64, CELL_SIZE * recorder.cell_size * (topology.radius * 2 + 2))
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
using RecyclePermanent
using RNGPool
using JLD2

include("../models/gaussianModel.jl")

setRNGs(12345)

θ = Gaussianθ(0.0, 1.5, 2.5, 0.5)
model = makeGaussianModel(θ)
ys = simulateModel(model, 1000)

@save "data/gaussianData.jld" ys θ


text:
@with_kw struct BlochBoundaryCondition <: AbstractBoundaryCondition
    k⃗::Vec
end

function BlochBoundaryCondition(planewave_source::PlaneWaveSource)
    # I think we need to include the refractive index of the surrounding medium
    BlochBoundaryCondition(k⃗₀(planewave_source))
end


text:
# -------- Opening a window ---------------
# Forward reference for @cfunction
function windowEventWatcher end
const window_event_watcher_cfunc = Ref(Ptr{Nothing}(0))

const window_paused = Threads.Atomic{UInt8}(0) # Whether or not the game should be running (if lost focus)
const frame_timer = WallTimer()

function makeWinRenderer(title = "My Julia SDL Game",
        min_win_dims::Union{Nothing,AbstractDims} = nothing)
    global winWidth, winHeight, winWidth_highDPI, winHeight_highDPI

    win = SDL2.CreateWindow(title,
        Int32(SDL2.WINDOWPOS_CENTERED()), Int32(SDL2.WINDOWPOS_CENTERED()), winWidth[], winHeight[],
        UInt32(SDL2.WINDOW_ALLOW_HIGHDPI|SDL2.WINDOW_OPENGL|SDL2.WINDOW_RESIZABLE|SDL2.WINDOW_SHOWN));
    if min_win_dims !== nothing
        SDL2.SetWindowMinimumSize(win, min_win_dims.w, min_win_dims.h)
    end
    window_event_watcher_cfunc[] = @cfunction(windowEventWatcher, Cint, (Ptr{Nothing}, Ptr{SDL2.Event}))
    SDL2.AddEventWatch(window_event_watcher_cfunc[], win);

    # Find out how big the created window actually was (depends on the system):
    winWidth[], winHeight[], winWidth_highDPI[], winHeight_highDPI[] = getWindowSize(win)
    #cam.w[], cam.h[] = winWidth_highDPI, winHeight_highDPI

    renderer = SDL2.CreateRenderer(win, Int32(-1), UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))
    SDL2.SetRenderDrawBlendMode(renderer, UInt32(SDL2.BLENDMODE_BLEND))
    return win,renderer
end

# This huge function handles all window events. I believe it needs to be a
# callback instead of just the regular pollEvent because the main thread is
# paused while resizing, whereas this callback continues to trigger.
function windowEventWatcher(data_ptr::Ptr{Cvoid}, event_ptr::Ptr{SDL2.Event})::Cint
    global winWidth, winHeight, cam, window_paused, renderer, win
    ev = unsafe_load(event_ptr, 1)
    ee = ev._Event
    t = UInt32(ee[4]) << 24 | UInt32(ee[3]) << 16 | UInt32(ee[2]) << 8 | ee[1]
    t = SDL2.Event(t)
    if (t == SDL2.WindowEvent)
        event = unsafe_load( Ptr{SDL2.WindowEvent}(pointer_from_objref(ev)) )
        winevent = event.event;  # confusing, but that's what the field is called.
        if (winevent == SDL2.WINDOWEVENT_RESIZED || winevent == SDL2.WINDOWEVENT_SIZE_CHANGED)
            curPaused = window_paused[]
            window_paused[] = 1  # Stop game playing so resizing doesn't cause problems.
            winID = event.windowID
            eventWin = SDL2.GetWindowFromID(winID);
            if (eventWin == data_ptr)
                w,h,w_highDPI,h_highDPI = getWindowSize(eventWin)
                winWidth[], winHeight[] = w, h
                winWidth_highDPI[], winHeight_highDPI[] = w_highDPI, h_highDPI
                cam.w[], cam.h[] = winWidth[], winHeight[]
                recenterButtons!()
            end
            # Note: render after every resize event. I tried limiting it with a
            # timer, but it's hard to tune (too infrequent and the screen
            # blinks) & it didn't seem to reduce cpu significantly.
            render(sceneStack[end], renderer, eventWin)
            SDL2.GL_SwapWindow(eventWin);
            window_paused[] = curPaused  # Allow game to resume now that resizing is done.
        elseif (winevent == SDL2.WINDOWEVENT_FOCUS_LOST || winevent == SDL2.WINDOWEVENT_HIDDEN || winevent == SDL2.WINDOWEVENT_MINIMIZED)
            # Stop game playing so resizing doesn't cause problems.
            #if !debug  # For debug builds, allow editing while playing
                window_paused[] = 1
            #end
        elseif (winevent == SDL2.WINDOWEVENT_FOCUS_GAINED || winevent == SDL2.WINDOWEVENT_SHOWN)
            window_paused[] = 0
        end
        # Note that window events pause the game, so at the end of any window
        # event, restart the timer so it doesn't have a HUGE frame.
        start!(frame_timer)
    end
    return 0
end

function getWindowSize(win)
    w,h,w_highDPI,h_highDPI = Int32[0],Int32[0],Int32[0],Int32[0]
    SDL2.GetWindowSize(win, w, h)
    SDL2.GL_GetDrawableSize(win, w_highDPI, h_highDPI)
    return w[],h[],w_highDPI[],h_highDPI[]
end

# Having a QuitException is useful for testing, since an exception will simply
# pause the interpreter. For release builds, the catch() block will call quitSDL().
struct QuitException <: Exception end

function quitSDL(win)
    # Need to close the callback before quitting SDL to prevent it from hanging
    # https://github.com/n0name/2D_Engine/issues/3
    SDL2.DelEventWatch(window_event_watcher_cfunc[], win);
    SDL2.Mix_CloseAudio()
    SDL2.TTF_Quit()
    SDL2.Quit()
end


text:
mutable struct JuliaProgram <: Program
	name::String
	id_file::String
	info_before::String
	info_after::String
    cmd_dependencies::Vector{CmdDependency}
	inputs::Vector{String}
	input_types::Vector{DataType}
	default_inputs::Vector
	validate_inputs::Function
	prerequisites::Function
    main::Function
	infer_outputs::Function
	outputs::Vector{String}
	output_types::Vector{DataType}
	default_outputs::Vector
	validate_outputs::Function
	wrap_up::Function
end

"""
	JuliaProgram <: Program

	JuliaProgram(;
		name::String               = "Unnamed Command Program",
		id_file::String            = "",
		info_before::String        = "auto",
		info_after::String         = "auto",
		cmd_dependencies::Vector{CmdDependency} = Vector{CmdDependency}(),
		inputs                     = Vector{String}(),
		validate_inputs::Function  = do_nothing,  # positional arguments: inputs::Dict{String}
		prerequisites::Function    = do_nothing,  # positional arguments: inputs, outputs::Dict{String}
		main::Function             = do_nothing,  # positional arguments: inputs, outputs::Dict{String}
		infer_outputs::Function    = do_nothing,  # positional arguments: inputs::Dict{String}
		outputs                    = Vector{String}(),
		validate_outputs::Function = do_nothing  # positional arguments: outputs::Dict{String},
		wrap_up::Function          = do_nothing  # positional arguments: inputs, outputs::Dict{String}
	) -> JuliaProgram

Julia program template. To run a `JuliaProgram`, use `run(::JuliaProgram; kwargs...).`

# Arguments

- `name::String`: Program name.

- `id_file::String`: The prefix of *run ID file*. To prevent from running the program with the same inputs and outputs twice, it will generate a unique *run ID file* after a successful run.

- `info_before::String`: Print it when starting the program.

- `info_after::String`: Print it when finishing the program.

- `cmd_dependencies::Vector{CmdDependency}`: Any command dependencies used in the program.

- `inputs` and `outputs`: Elements (or vectors containing elements) in the following format: (1) `keyword` (2) `keyword => data_type` (3) `keyword => default_value` (4) `keyword => default_value => data_type`.

  `keyword` is an argument name, can be `String` or `Symbol`.

  `default_value` is optional. If set, users may not provide this argument when running. Elsewise, users have to provide it. Caution: `nothing` is preserved and means default value not set. If `String`, it can contain other keywords, but need to quote using '<>', such as `"<arg>.txt"`

  `data_type` is optional. If set, the value provided have to be this data type, or an error will throw.

  *HOW DOES THIS WORK?*

  > `JuliaProgram` stores a Julia function, with two arguments `inputs::Dict{String}, outputs::Dict{String}`. The keys of the two arguments should be set in `inputs::Vector{String}` and `outputs::Vector{String}`.
  > Caution: the returned value of `p.main` will be assigned to new `outputs`. Please ensure the returned value is `Dict{String}` with proper keys.

- `validate_inputs::Function`: A function to validate inputs. It takes *one* argument `Dict{String}` whose keys are the same as `inputs`. If validation fail, throw error or return false.

- `prerequisites`: A function to run just before the main command. It prepares necessary things, such as creating directories. It takes *two* arguments `Dict{String}` whose keys are the same as `inputs` and `outputs`, respectively.

- `main::Function`: The main julia function. It takes *two* arguments `Dict{String}` whose keys are the same as `inputs` and `outputs`, respectively.

  > Caution: the returned value of `p.main` will be assigned to new `outputs`. Please ensure the returned value is `Dict{String}` with proper keys.

- `infer_outputs::Function`: A function to infer outputs from inputs. It takes *one* argument `Dict{String}` whose keys are the same as `inputs`.

- `validate_outputs::Function`: A function to validate outputs. It takes *one* argument `Dict{String}` whose keys are the same as `outputs`. If validation fail, throw error or return false.

- `wrap_up::Function`: the last function to run. It takes *two* arguments `Dict{String}` whose keys are the same as `inputs` and `outputs`, respectively.

# Example

	p = JuliaProgram(
		id_file = "id_file",
		inputs = [
			"a",
			"b" => Int
		],
		outputs =
			"c" => "<a>.<b>"
		,
		main = (inputs, outputs) -> begin
			a = inputs["a"]
			b = inputs["b"]
			println("inputs are ", a, " and ", b)
			println("You can also use info in outputs: ", outputs["c"])
	        println("The returned value will be assigned to a new outputs")

	        return Dict{String,Any}("c" => b^2)
		end
	)

	inputs = Dict(
		"a" => `in1`,
		"b" => 2
	)

	outputs = Dict(
		"c" => "out"
	)

	success, outputs = run(p, inputs, outputs;
		touch_run_id_file = false
	) # outputs will be refreshed
"""
function JuliaProgram(;
	name::String               = "Julia Program",
	id_file::String            = "",
	info_before::String        = "auto",
	info_after::String         = "auto",
	cmd_dependencies::Vector{CmdDependency} = Vector{CmdDependency}(),
	inputs                     = Vector{String}(),
	validate_inputs::Function  = do_nothing,  # positional arguments: inputs::Dict{String}
	prerequisites::Function    = do_nothing,  # positional arguments: inputs, outputs::Dict{String}
	main::Function             = do_nothing,  # positional arguments: inputs, outputs::Dict{String},
	infer_outputs::Function    = do_nothing,  # positional arguments: inputs::Dict{String}
	outputs                    = Vector{String}(),
	validate_outputs::Function = do_nothing,  # positional arguments: outputs::Dict{String}
	wrap_up::Function          = do_nothing  # positional arguments: inputs, outputs::Dict{String}
)
	inputs, input_types, default_inputs = parse_default(inputs)
	outputs, output_types, default_outputs = parse_default(outputs)

	JuliaProgram(
		name,
		id_file,
		info_before,
		info_after,
		cmd_dependencies,
		inputs,
		input_types,
		default_inputs,
		validate_inputs,
		prerequisites,
		main,
		infer_outputs,
		outputs,
		output_types,
		default_outputs,
		validate_outputs,
		wrap_up
	)
end

"""
------------

	run(
		p::JuliaProgram;
		inputs=Dict{String}(),
		outputs=Dict{String}(),
		dir::AbstractString="",
		check_dependencies::Bool=true,
		skip_when_done::Bool=true,
		touch_run_id_file::Bool=true,
		verbose::Bool=true,
		dry_run::Bool=false,
		stdout=nothing,
		stderr=nothing,
		stdlog=nothing,
		append::Bool=false
	) -> (success::Bool, outputs::Dict{String})

	run(p::JuliaProgram, inputs, outputs; kwargs...)

	run(p::JuliaProgram, inputs; kwargs...)
	)  # only usable when `p.infer_outputs` is defined, or default outputs are set in `p`.

Run Julia Program (JuliaProgram) using specified `inputs` and `outputs`.

Return `(success::Bool, outputs::Dict{String})`

- `p::JuliaProgram`: the command program template.

- `inputs::Dict{String}` and `outputs::Dict{String}`: `JuliaProgram` stores a Julia function, with two arguments `inputs::Dict{String}, outputs::Dict{String}`. The keys of the two arguments should be the same as `p.inputs::Vector{String}` and `p.outputs::Vector{String}`.

  > If data types of `inputs` and `outputs` are not `Dict{String}`, they will be converted as far as possible. If the conversion fails, program will throw an error.
  > *Caution:* the returned value of `p.main` will be assigned to new `outputs`. Please ensure the returned value of `p.main` is `Dict{String}` with proper keys.

- `dir::AbstractString = ""`: working directory to run the program.

- `check_dependencies::Bool = true`: check dependencies for `p` (`p.cmd_dependencies`).

- `skip_when_done::Bool = true`: Skip running the program and return `true` if it has been done before (the `run_id_file` exists and `p.validate_outputs(outputs)` passes.)

- `touch_run_id_file::Bool = true`: If `true`, touch a unique run ID file, which indicate the program is successfully run with given inputs and outputs. If `false`, the next time running the program, `skip_when_done=true` will not take effect.

- `verbose::Bool = true`: If `true`, print info and error messages. If `false`, print error messages only.

- `dry_run::Bool = false`: do not run the command, return `(fake_outputs::Dict{String}, run_id_file::String)`.

- `stdout`, `stderr`, `stdlog` and `append`: Redirect the program outputs to files. `stdlog` is the Julia logging of `@info`, `@warn`, `@error`, etc. Caution: If the original function (`p.main`) has redirection, arguments defined here might not be effective for it.

!!! note
    Redirecting in Julia are not thread safe, so unexpected redirection might be happen if you are running programs in different `Tasks` or multi-thread mode.

### Workflow

1. Go to the working directory. Establish redirection. (`dir`, `stdout`, `stderr`, `stdlog`, `append`).

2. Validate compatibility between `p` and `inputs/outputs`.

3. Check whether the program has run before. (`skip_when_done`, `p.validate_outputs(outputs)`)

4. Check command dependencies. (`check_dependencies`, `p.cmd_dependencies`)

5. Validate `inputs`. (`p.validate_inputs(inputs)`)

6. Preparing before running main command. (`p.prerequisites(inputs, outputs)`)

7. Run main function and ***the returned value will be assigned to new `outputs`***. (`outputs = p.main(inputs, outputs)`)

8. Validate `outputs`. (`p.validate_outputs(outputs)`)

9. Wrap up. (`p.wrap_up(inputs, outputs)`)

10. Success, touch run id file, and return `(success::Bool, outputs::Dict{String})`. (`touch_run_id_file::Bool`)

# Example

	p = JuliaProgram(
		id_file = "id_file",
		inputs = ["a", "b"],
		outputs = ["out"],
		main = (inputs, outputs) -> begin
			a = inputs["a"]
			b = inputs["b"]
			println("inputs are ", a, " and ", b)
			println("You can also use info in outputs: ", outputs["out"])

			println("The returned value will be assigned to a new outputs")
			return Dict{String,Any}("out" => b^2)
		end
	)

	inputs = Dict(
		"a" => `in1`,
		"b" => 2
	)

	outputs = Dict(
		"out" => "will_be_replaced"
	)

	success, outputs = run(p, inputs, outputs;
		touch_run_id_file = false
	) # outputs will be refreshed
"""
function _run(
	p::JuliaProgram;
	inputs=Dict{String, Any}(),
	outputs=Dict{String, Any}(),
	skip_when_done::Bool=true,
	check_dependencies::Bool=true,
	verbose::Bool=true,
	touch_run_id_file::Bool=true,
	dry_run::Bool=false
)
	# input/output completion
	inputs, outputs = xxputs_completion_and_check(p, inputs, outputs)

	# run id based on inputs and outputs
	run_id = generate_run_uuid(inputs, outputs)
	run_id_file = p.id_file * "." * string(run_id)

	# whether dry run
	if dry_run
		@goto dry_run_start
	end

	if verbose
		if p.info_before == "auto" || p.info_before == ""
			@info timestamp() * "Starting program: $(p.name)" run_id inputs outputs
		else
			@info timestamp() * p.info_before run_id inputs outputs
		end
	end

	# skip when done
	if skip_when_done && isfile(run_id_file) && isok(p.validate_outputs(outputs))
		@warn timestamp() * "Skipping finished program: $(p.name)" run_id inputs outputs
		return true, outputs
	end

	# check dependencies
	if check_dependencies
		foreach(check_dependency, p.cmd_dependencies)
	end

	# preparation: remove run id file
	isfile(run_id_file) && rm(run_id_file)

	# preparation: validate inputs
	try
		isok(p.validate_inputs(inputs)) || error("ProgramInputValidationError: $(p.name): the prerequisites function returns false.")
	catch e
		@error timestamp() * "ProgramInputValidationError: $(p.name): fail to validate inputs (before running the main command). See error messages above." validation_function=p.validate_inputs run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramInputValidationError")
		return false, outputs
	end

	@label dry_run_start

	if dry_run
		return (outputs, run_id_file)
		# dry_run stops
	end

	# preparation: run function prerequisites before the main command
	try
		isok(p.prerequisites(inputs, outputs)) || error("ProgramPrerequisitesError: $(p.name): the prerequisites function returns false.")
	catch e
		@error timestamp() * "ProgramPrerequisitesError: $(p.name): fail to run the prerequisites function (before running the main command). See error messages above." prerequisites=p.prerequisites run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramPrerequisitesError")
		return false, outputs
	end

	# run the main command
	outputs = try
		p.main(inputs, outputs)
	catch e
		@error timestamp() * "ProgramRunningError: $(p.name): fail to run the main command. See error messages above." run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramRunningError")
		return false, outputs
	end

	# check type of outputs
	if !isa(outputs, Dict{String})
		@error timestamp() * "ProgramMainReturnValueError: $(p.name): the returned value is not a Dict{String}" run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramMainReturnValueError")
		return false, outputs
	end

	# check keys in outputs::Dict{String}
	check_outputs_keywords(p, outputs)

	# confirmation: validate outputs
	try
		isok(p.validate_outputs(outputs)) || error("ProgramOutputValidationError: $(p.name): the validation function returns false.")
	catch e
		@error timestamp() * "ProgramOutputValidationError: $(p.name): fail to validate outputs (after running the main command). See error messages above." validation_function=p.validate_outputs run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramOutputValidationError")
		return false, outputs
	end

	try
		isok(p.wrap_up(inputs, outputs)) || error("ProgramWrapUpError: $(p.name): the wrap_up function returns false.")
	catch e
		@error timestamp() * "ProgramWrapUpError: $(p.name): fail to run the wrap_up function. See error messages above." wrap_up=p.wrap_up run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramWrapUpError")
		return false, outputs
	end

	# touch the run_id_file
	touch_run_id_file && touch(run_id_file)

	if verbose
		if p.info_after == "auto" || p.info_after == ""
			@info timestamp() * "Finishing program: $(p.name)" run_id inputs outputs
		else
			@info timestamp() * p.info_after run_id inputs outputs
		end
	end
	return true, outputs
end


text:
include("struct/distance.jl")

"""
Essaie de regrouper des données en commençant par celles qui sont les plus proches.
Deux clusters de données peuvent être fusionnés en un cluster C s'il n'existe aucune données x_i pour aucune caractéristique j qui intersecte l'intervalle représenté par les bornes minimale et maximale de C pour j (x_i,j n'appartient pas à [min_{x_k dans C} x_k,j ; max_{k dans C} x_k,j]).

Entrées :
- x : caractéristiques des données d'entraînement
- y : classe des données d'entraînement
- percentage : le nombre de clusters obtenu sera égal à n * percentage

Sorties :
- un tableau de Cluster constituant une partition de x
"""
function exactMerge(x, y)

    n = length(y)
    m = length(x[1,:])

    # Liste de tous les clusters obtenus
    # (les données non regroupées seront seules dans un cluster)
    clusters = Vector{Cluster}([])

    # Initialement, chaque donnée est dans un cluster
    for dataId in 1:size(x, 1)
        push!(clusters, Cluster(dataId, x, y))
    end

    # Id du cluster de chaque donnée dans clusters
    # (clusters[clusterId[i]] est le cluster contenant i)
    # (clusterId[i] = i, initialement)
    clusterId = collect(1:n)

    # Distances entre des couples de données de même classe
    distances = Vector{Distance}([])

    # Pour chaque couple de données de même classe
    for id1 in 1:n-1
        for id2 in id1+1:n
            if y[id1] == y[id2]
                # Ajoute leur distance
                push!(distances, Distance(id1, id2, x))
            end
        end
    end

    # Trie des distances par ordre croissant
    sort!(distances, by = v -> v.distance)

    # Pour chaque distance
    for distance in distances

        # Si les deux données associées ne sont pas déjà dans le même cluster
        cId1 = clusterId[distance.ids[1]]
        cId2 = clusterId[distance.ids[2]]

        if cId1 != cId2
            c1 = clusters[cId1]
            c2 = clusters[cId2]

            # Si leurs clusters satisfont les conditions de fusion
            if canMerge(c1, c2, x, y)

                # Les fusionner
                merge!(c1, c2)
                for id in c2.dataIds
                    clusterId[id]= cId1
                end

                # Vider le second cluster
                empty!(clusters[cId2].dataIds)
            end
        end
    end

    # Retourner tous les clusters non vides
    return filter(x -> length(x.dataIds) > 0, clusters)
end

"""
Regroupe des données en commençant par celles qui sont les plus proches jusqu'à ce qu'un certain pourcentage de clusters soit atteint

Entrées :
- x : caractéristiques des données
- y : classe des données
- gamma : le regroupement se termine quand il reste un nombre de clusters < n * gamma ou que plus aucun regroupement n'est possible

Sorties :
- un tableau de Cluster constituant une partition de x
"""
function simpleMerge(x, y, gamma)

    n = length(y)
    m = length(x[1,:])

    # Liste de tous les clusters obtenus
    # (les données non regroupées seront seules dans un cluster)
    clusters = Vector{Cluster}([])

    # Initialement, chaque donnée est dans un cluster
    for dataId in 1:size(x, 1)
        push!(clusters, Cluster(dataId, x, y))
    end

    # Id du cluster de chaque donnée dans clusters
    # (clusters[clusterId[i]] est le cluster contenant i)
    # (clusterId[i] = i, initialement)
    clusterId = collect(1:n)

    # Distances entre des couples de données de même classe
    distances = Vector{Distance}([])

    # Pour chaque couple de données de même classe
    for id1 in 1:n-1
        for id2 in id1+1:n
            if y[id1] == y[id2]
                # Ajoute leur distance
                push!(distances, Distance(id1, id2, x))
            end
        end
    end

    # Trie des distances par ordre croissant
    sort!(distances, by = v -> v.distance)

    remainingClusters = n
    distanceId = 1

    # Pour chaque distance et tant que le nombre de cluster souhaité n'est pas atteint
    while distanceId <= length(distances) && remainingClusters > n * gamma

        distance = distances[distanceId]
        cId1 = clusterId[distance.ids[1]]
        cId2 = clusterId[distance.ids[2]]

        # Si les deux données associées ne sont pas déjà dans le même cluster
        if cId1 != cId2
            remainingClusters -= 1

            # Fusionner leurs clusters
            c1 = clusters[cId1]
            c2 = clusters[cId2]
            merge!(c1, c2)
            for id in c2.dataIds
                clusterId[id]= cId1
            end

            # Vider le second cluster
            empty!(clusters[cId2].dataIds)
        end
        distanceId += 1
    end

    # Retourner tous les clusters non vides
    return filter(x -> length(x.dataIds) > 0, clusters)
end

"""
Test si deux clusters peuvent être fusionnés tout en garantissant l'optimalité

Entrées :
- c1 : premier cluster
- c2 : second cluster
- x  : caractéristiques des données d'entraînement
- y  : classe des données d'entraînement

Sorties :
- vrai si la fusion est possible ; faux sinon.
"""
function canMerge(c1::Cluster, c2::Cluster, x::Matrix{Float64}, y::Vector{Int})

    # Calcul des bornes inférieures si c1 et c2 étaient fusionnés
    mergedLBounds = min.(c1.lBounds, c2.lBounds)

    # Calcul des bornes supérieures si c1 et c2 étaient fusionnés
    mergedUBounds = max.(c1.uBounds, c2.uBounds)

    n = size(x, 1)
    id = 1
    canMerge = true

    # Tant que l'ont a pas vérifié que toutes les données n'intersectent la fusion de c1 et c2 sur aucune feature
    while id <= n && canMerge

        data = x[id, :]

        # Si la donnée n'est pas dans c1 ou c2 mais intersecte la fusion de c1 et c2 sur au moins une feature
        if !(id in c1.dataIds) && !(id in c2.dataIds) && isInABound(data, mergedLBounds, mergedUBounds)
            canMerge = false
        end

        id += 1
    end

    return canMerge
end

"""
Test si une donnée intersecte des bornes pour au moins une caractéristique

Entrées :
- v : les caractéristique de la donnée
- lowerBounds : bornes inférieures pour chaque caractéristique
- upperBounds : bornes supérieures pour chaque caractéristique

Sorties :
- vrai s'il y a intersection ; faux sinon.
"""
function isInABound(v::Vector{Float64}, lowerBounds::Vector{Float64}, upperBounds::Vector{Float64})
    isInBound = false

    featureId = 1

    # Tant que toutes les features n'ont pas été testées et qu'aucune intersection n'a été trouvée
    while !isInBound && featureId <= length(v)

        # S'il y a intersection
        if v[featureId] >= lowerBounds[featureId] && v[featureId] <= upperBounds[featureId]
            isInBound = true
        end
        featureId += 1
    end

    return isInBound
end

"""
Fusionne deux clusters

Entrées :
- c1 : premier cluster
- c2 : second cluster

Sorties :
- aucune, c'est le cluster en premier argument qui contiendra le second
"""
function merge!(c1::Cluster, c2::Cluster)

    append!(c1.dataIds, c2.dataIds)
    c1.x = vcat(c1.x, c2.x)
    c1.lBounds = min.(c1.lBounds, c2.lBounds)
    c1.uBounds = max.(c1.uBounds, c2.uBounds)
end

max_stars_repo_name: boumpteryx/RODM

max_issues_repo_name: boumpteryx/RODM

max_forks_repo_name: boumpteryx/RODM


text:
using ConvexTests, ECOS

@info "Starting ECOS tests"

do_tests("ECOS", () -> ECOS.Optimizer(verbose=0); exclude = [r"mip", r"sdp"])


text:
######### Stan program example  ###########

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoulli = "
  data {
    int<lower=0> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
      y ~ bernoulli(theta);
  }
  "

  bernoullidata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  global stanmodel, rc, dim
  stanmodel = Stanmodel(Variational(), name="bernoulli", model=bernoulli);

  rc, sim1 = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    monitor = ["theta"]
    sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
    Mamba.describe(sim)
    println()

    ## Plotting
    p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
    draw(p, ncol=1, nrow=4, filename="$(stanmodel.name)-variationalplot", fmt=:svg)
    draw(p, ncol=1, nrow=4, filename="$(stanmodel.name)-variationalplot", fmt=:pdf)
  end
end # cd


text:
function plot!(plot::Text)
    # attach a function to any text that calculates the glyph layout and stores it
    glyphcollection = lift(plot[1], plot.textsize, plot.font, plot.align,
            plot.rotation, plot.justification, plot.lineheight,
            plot.color, plot.strokecolor, plot.strokewidth) do str,
                ts, f, al, rot, jus, lh, col, scol, swi
        ts = to_textsize(ts)
        f = to_font(f)
        rot = to_rotation(rot)
        col = to_color(col)
        scol = to_color(scol)

        layout_text(str, ts, f, al, rot, jus, lh, col, scol, swi)
    end

    text!(plot, glyphcollection; plot.attributes...)

    plot
end

# TODO: is this necessary? there seems to be a recursive loop with the above
# function without these two interceptions, but I didn't need it before merging
# everything into the monorepo...
plot!(plot::Text{<:Tuple{<:GlyphCollection}}) = plot
plot!(plot::Text{<:Tuple{<:AbstractArray{<:GlyphCollection}}}) = plot

function plot!(plot::Text{<:Tuple{<:AbstractArray{<:AbstractString}}})

    glyphcollections = Observable(GlyphCollection[])
    position = Observable{Any}(nothing)
    rotation = Observable{Any}(nothing)

    onany(plot[1], plot.textsize, plot.position,
            plot.font, plot.align, plot.rotation, plot.justification,
            plot.lineheight, plot.color, plot.strokecolor, plot.strokewidth) do str,
                    ts, pos, f, al, rot, jus, lh, col, scol, swi

        ts = to_textsize(ts)
        f = to_font(f)
        rot = to_rotation(rot)
        col = to_color(col)
        scol = to_color(scol)

        gcs = GlyphCollection[]
        broadcast_foreach(str, ts, f, al, rot, jus, lh, col, scol, swi) do str,
                ts, f, al, rot, jus, lh, col, scol, swi
            subgl = layout_text(str, ts, f, al, rot, jus, lh, col, scol, swi)
            push!(gcs, subgl)
        end
        position.val = pos
        rotation.val = rot
        glyphcollections[] = gcs
    end

    # run onany once to initialize
    notify(plot[1])

    text!(plot, glyphcollections; position = plot.position, rotation = rotation,
        model = plot.model, offset = plot.offset, markerspace = plot.markerspace,
        visible=plot.visible, space = plot.space)

    plot
end

# overload text plotting for a vector of tuples of a string and a point each
function plot!(plot::Text{<:Tuple{<:AbstractArray{<:Tuple{<:AbstractString, <:Point}}}})
    strings_and_positions = plot[1]

    strings = Observable(first.(strings_and_positions[]))
    positions = Observable(to_ndim.(Ref(Point3f), last.(strings_and_positions[]), 0))

    attrs = plot.attributes
    pop!(attrs, :position)

    t = text!(plot, strings; position = positions, attrs...)

    # update both text and positions together
    on(strings_and_positions) do str_pos
        strs = first.(str_pos)
        poss = to_ndim.(Ref(Point3f), last.(str_pos), 0)
        # first mutate strings without triggering redraw
        t[1].val = strs
        # then update positions with trigger
        positions[] = poss
    end
    plot
end

function plot!(plot::Text{<:Tuple{<:Union{LaTeXString, AbstractVector{<:LaTeXString}}}})

    # attach a function to any text that calculates the glyph layout and stores it
    lineels_glyphcollection_offset = lift(plot[1], plot.textsize, plot.align, plot.rotation,
            plot.model, plot.color, plot.strokecolor, plot.strokewidth, plot.position) do latexstring,
                ts, al, rot, mo, color, scolor, swidth, _

        ts = to_textsize(ts)
        rot = to_rotation(rot)

        if latexstring isa AbstractVector
            tex_elements = []
            glyphcollections = GlyphCollection[]
            offsets = Point2f[]
            broadcast_foreach(latexstring, ts, al, rot, color, scolor, swidth) do latexstring,
                ts, al, rot, color, scolor, swidth

                te, gc, offs = texelems_and_glyph_collection(latexstring, ts,
                    al[1], al[2], rot, color, scolor, swidth)
                push!(tex_elements, te)
                push!(glyphcollections, gc)
                push!(offsets, offs)
            end
            tex_elements, glyphcollections, offsets
        else
            tex_elements, glyphcollections, offsets = texelems_and_glyph_collection(latexstring, ts,
                al[1], al[2], rot, color, scolor, swidth)
        end
    end

    glyphcollection = @lift($lineels_glyphcollection_offset[2])

    linepairs = Observable(Tuple{Point2f, Point2f}[])
    linewidths = Observable(Float32[])

    scene = Makie.parent_scene(plot)

    onany(lineels_glyphcollection_offset, scene.camera.projectionview) do (allels, gcs, offs), projview

        inv_projview = inv(projview)
        pos = plot.position[]
        ts = plot.textsize[]
        rot = plot.rotation[]

        ts = to_textsize(ts)
        rot = convert_attribute(rot, key"rotation"())

        empty!(linepairs.val)
        empty!(linewidths.val)

        # for the vector case, allels is a vector of vectors
        # so for broadcasting the single vector needs to be wrapped in Ref
        if gcs isa GlyphCollection
            allels = [allels]
        end
        broadcast_foreach(allels, offs, pos, ts, rot) do allels, offs, pos, ts, rot
            offset = Point2f(pos)

            els = map(allels) do el
                el[1] isa VLine || el[1] isa HLine || return nothing

                t = el[1].thickness * ts
                p = el[2]

                ps = if el[1] isa VLine
                    h = el[1].height
                    (Point2f(p[1], p[2]) .* ts, Point2f(p[1], p[2] + h) .* ts) .- Ref(offs)
                else
                    w = el[1].width
                    (Point2f(p[1], p[2]) .* ts, Point2f(p[1] + w, p[2]) .* ts) .- Ref(offs)
                end
                ps = Ref(rot) .* to_ndim.(Point3f, ps, 0)
                # TODO the points need to be projected to work inside Axis
                # ps = project ps with projview somehow

                ps = Point2f.(ps) .+ Ref(offset)
                ps, t
            end
            pairs = filter(!isnothing, els)
            append!(linewidths.val, repeat(last.(pairs), inner = 2))
            append!(linepairs.val, first.(pairs))
        end
        notify(linepairs)
    end

    notify(plot.position)

    text!(plot, glyphcollection; plot.attributes...)
    linesegments!(plot, linepairs, linewidth = linewidths, color = plot.color)

    plot
end

function texelems_and_glyph_collection(str::LaTeXString, fontscale_px, halign, valign,
        rotation, color, strokecolor, strokewidth)

    rot = convert_attribute(rotation, key"rotation"())

    all_els = generate_tex_elements(str)
    els = filter(x -> x[1] isa TeXChar, all_els)

    # hacky, but attr per char needs to be fixed
    fs = Vec2f(first(fontscale_px))

    scales_2d = [Vec2f(x[3] * Vec2f(fs)) for x in els]

    chars = [x[1].char for x in els]
    fonts = [x[1].font for x in els]

    extents = [FreeTypeAbstraction.get_extent(f, c) for (f, c) in zip(fonts, chars)]

    bboxes = map(extents, fonts, scales_2d) do ext, font, scale
        unscaled_hi_bb = FreeTypeAbstraction.height_insensitive_boundingbox(ext, font)
        return Rect2f(
            origin(unscaled_hi_bb) * scale,
            widths(unscaled_hi_bb) * scale
        )
    end

    basepositions = [to_ndim(Vec3f, fs, 0) .* to_ndim(Point3f, x[2], 0)
        for x in els]

    bb = isempty(bboxes) ? BBox(0, 0, 0, 0) : begin
        mapreduce(union, zip(bboxes, basepositions)) do (b, pos)
            Rect2f(Rect3f(b) + pos)
        end
    end

    xshift = if halign == :center
        width(bb) ./ 2
    elseif halign == :left
        minimum(bb)[1]
    elseif halign == :right
        maximum(bb)[1]
    end

    yshift = if valign == :center
        maximum(bb)[2] - (height(bb) / 2)
    elseif valign == :top
        maximum(bb)[2]
    else
        minimum(bb)[2]
    end

    positions = basepositions .- Ref(Vec3f(xshift, yshift, 0))
    positions .= Ref(rot) .* positions

    pre_align_gl = GlyphCollection(
        chars,
        fonts,
        Point3f.(positions),
        extents,
        scales_2d,
        rot,
        color,
        strokecolor,
        strokewidth,
    )

    all_els, pre_align_gl, Point2f(xshift, yshift)
end

MakieLayout.iswhitespace(l::LaTeXString) = MakieLayout.iswhitespace(replace(l.s, '$' => ""))


text:
function inex(digits)
#
#  test polynomial suggested by Goedecker
#
    p = poly([(10/11)*[1;1;1;1;1]; (20/11)*[1;1;1]; (30/11)*[1;1]])
    p = round.(10^digits*reverse(p.a)) ./ 10^digits
    z = reverse([10/11 5; 20/11 3; 30/11 2]; dims=1)
    p, PolyZeros(z)
end


text:
# restrict
function restrict(img::Union{WarpedView, InvWarpedView}, dims::Integer)
    # preserve axes information but eagerly evaluate the transformation
    # TODO: this eager evaluation destroy the laziness of `WarpedView`
    #       and `InvWarpedView`.
    restrict(OffsetArray(collect(img), axes(img)), dims)
end

# imresize
imresize(original::AbstractArray, dim1::T, dimN::T...; kwargs...) where T<:Union{Integer,AbstractUnitRange} = imresize(original, (dim1,dimN...); kwargs...)
function imresize(original::AbstractArray; ratio, kwargs...)
    all(ratio .> 0) || throw(ArgumentError("ratio $ratio should be positive"))
    new_size = ceil.(Int, size(original) .* ratio) # use ceil to avoid 0
    imresize(original, new_size; kwargs...)
end

function imresize(original::AbstractArray, short_size::Union{Indices{M},Dims{M}}; kwargs...) where M
    len_short = length(short_size)
    len_short > ndims(original) && throw(DimensionMismatch("$short_size has too many dimensions for a $(ndims(original))-dimensional array"))
    new_size = ntuple(i -> (i > len_short ? odims(original, i, short_size) : short_size[i]), ndims(original))
    imresize(original, new_size; kwargs...)
end

function imresize(original::AbstractArray, itp::Union{Interpolations.Degree,Interpolations.InterpolationType}, args...; kwargs...)
    throw(ArgumentError("Positional parameter for interpolation method is not supported; use keyword `method` instead."))
end

odims(original, i, short_size::Tuple{Integer,Vararg{Integer}}) = size(original, i)
odims(original, i, short_size::Tuple{}) = axes(original, i)
odims(original, i, short_size) = oftype(first(short_size), axes(original, i))

"""
    imresize(img, sz; [method]) -> imgr
    imresize(img, inds; [method]) -> imgr
    imresize(img; ratio, [method]) -> imgr

upsample/downsample the image `img` to a given size `sz` or axes `inds` using interpolations. If
`ratio` is provided, the output size is then `ceil(Int, size(img).*ratio)`.

!!! tip
    This interpolates the values at sub-pixel locations. If you are shrinking the image, you risk
    aliasing unless you low-pass filter `img` first.

# Arguments

- `img`: the input image array
- `sz`: the size of output array
- `inds`: the axes of output array
  If `inds` is passed, the output array `imgr` will be `OffsetArray`.

# Parameters

!!! info
    To construct `method`, you may need to load `Interpolations` package first.

- `ratio`: the upsample/downsample ratio used.
  The output size is `ceil(Int, size(img).*ratio)`. If `ratio` is larger than `1`, it is
  an upsample operation. Otherwise it is a downsample operation. `ratio` can also be a tuple,
  in which case `ratio[i]` specifies the resize ratio at dimension `i`.
- `method::InterpolationType`:
  specify the interpolation method used for reconstruction. conveniently, `methold` can
  also be a `Degree` type, in which case a `BSpline` object will be created.
  For example, `method = Linear()` is equivalent to `method = BSpline(Linear())`.

# Examples

```julia
using ImageTransformations, TestImages, Interpolations

img = testimage("lighthouse") # 512*768

# pass integers as size
imresize(img, 256, 384) # 256*384
imresize(img, (256, 384)) # 256*384
imresize(img, 256) # 256*768

# pass indices as axes
imresize(img, 1:256, 1:384) # 256*384
imresize(img, (1:256, 1:384)) # 256*384
imresize(img, (1:256, )) # 256*768

# pass resize ratio
imresize(img, ratio = 0.5) #256*384
imresize(img, ratio = (2, 1)) # 1024*768

# pass `CenterPoint`
imresize(img, (256, 384), CenterPoint(10, 15)) # 256*384, img[10,15] ≈ imgr[10,15]
imresize(img, CenterPoint(10, 15), ratio = 0.5) # 256*384, img[10,15] ≈ imgr[10,15]

# use different interpolation method
imresize(img, (256, 384), method=Linear()) # 256*384 bilinear interpolation
imresize(img, (256, 384), method=Lanczos4OpenCV()) # 256*384 OpenCV-compatible Lanczos 4 interpolation
```

For downsample with `ratio=0.5`, [`restrict`](@ref) is a much faster two-fold implementation that
you can use.
"""
function imresize(original::AbstractArray{T,0}, new_inds::Tuple{}; kwargs...) where T
    Tnew = imresize_type(first(original))
    copyto!(similar(original, Tnew), original)
end

function imresize(original::AbstractArray{T,N}, new_size::Dims{N}; kwargs...) where {T,N}
    Tnew = imresize_type(first(original))
    inds = axes(original)
    if map(length, inds) == new_size
        dest = similar(original, Tnew, new_size)
        if axes(dest) == inds
            copyto!(dest, original)
        else
            copyto!(dest, CartesianIndices(axes(dest)), original, CartesianIndices(inds))
        end
    else
        imresize!(similar(original, Tnew, new_size), original; kwargs...)
    end
end

function imresize(original::AbstractArray{T,N}, new_inds::Indices{N}; kwargs...) where {T,N}
    Tnew = imresize_type(first(original))
    if axes(original) == new_inds
        copyto!(similar(original, Tnew), original)
    else
        imresize!(similar(original, Tnew, new_inds), original; kwargs...)
    end
end

function imresize(original::AbstractArray{T,N}, new_size::Dims{N}, center_point::CenterPoint{N}; kwargs...) where {T,N}
    Tnew = imresize_type(first(original))
    cp = center_point.p
    checkbounds(original, cp)
    topleft = firstindex.(Ref(original), Tuple(1:N))
    offset = @. cp.I - new_size * (cp.I - topleft) ÷ $size(original) - 1
    newimage = OffsetArray(similar(original, Tnew, new_size), offset)
    imresize!(newimage, original; kwargs...)
end

function imresize(original::AbstractArray{T,N}, center_point::CenterPoint{N}; ratio, kwargs...) where {T,N}
    all(ratio .> 0) || throw(ArgumentError("ratio $ratio should be positive"))
    new_size = ceil.(Int, size(original) .* ratio) # use ceil to avoid 0
    imresize(original, new_size, center_point; kwargs...)
end
# To choose the output type, rather than forcing everything to
# Float64 by multiplying by 1.0, we exploit the fact that the scale
# changes correspond to integer ratios.  We mimic ratio arithmetic
# without actually using Rational (which risks promoting to a
# Rational type, too slow for image processing).
function imresize_type(c::Colorant)
    CT = base_colorant_type(c)
    isconcretetype(CT) && return CT # special 0-parameter colorant types: ARGB32, etc
    CT{eltype(imresize_type(Gray(c)))}
end
imresize_type(c::Gray) = Gray{imresize_type(gray(c))}
imresize_type(c::FixedPoint) = typeof(c)
imresize_type(c) = typeof((c*1)/1)

function imresize!(resized::AbstractArray{T,N}, original::AbstractArray{S,N}; method::Union{Interpolations.Degree,Interpolations.InterpolationType}=BSpline(Linear()), kwargs...) where {T,S,N}
    # FIXME: avoid allocation for interpolation
    itp = interpolate(original, wrap_BSpline(method))
    imresize!(resized, itp)
end

function imresize!(resized::AbstractArray{T,N}, original::AbstractInterpolation{S,N}) where {T,S,N}
    # Define the equivalent of an affine transformation for mapping
    # locations in `resized` to the corresponding position in
    # `original`. We take the viewpoint that a pixel at `i, j` is a
    # sensor that *integrates* the intensity over an area spanning
    # `i±0.5, j±0.5` (this is a good model of how a camera pixel
    # actually works). We then map the *outer corners* of the two
    # images to each other, i.e., in typical cases
    #     (0.5, 0.5) -> (0.5, 0.5)  (outer corner, top left)
    #     size(resized)+0.5 -> size(original)+0.5  (outer corner, lower right)
    # This ensures that both images cover exactly the same area.
    Ro, Rr = CartesianIndices(axes(original)), CartesianIndices(axes(resized))
    sf = map(/, (last(Ro)-first(Ro)).I .+ 1, (last(Rr)-first(Rr)).I .+ 1) # +1 for outer corners
    offset = map((io,ir,s)->io - 0.5 - s*(ir-0.5), first(Ro).I, first(Rr).I, sf)
    if all(x->x >= 1, sf)
        @inbounds for I in Rr
            I_o = map3((i,s,off)->s*i+off, I.I, sf, offset)
            resized[I] = original(I_o...)
        end
    else
        @inbounds for I in Rr
            I_o = clampR(map3((i,s,off)->s*i+off, I.I, sf, offset), Ro)
            resized[I] = original(I_o...)
        end
    end
    resized
end

# map isn't optimized for 3 tuple-arguments, so do it here
@inline map3(f, a, b, c) = (f(a[1], b[1], c[1]), map3(f, tail(a), tail(b), tail(c))...)
@inline map3(f, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()

function clampR(I::NTuple{N}, R::CartesianIndices{N}) where N
    map3(clamp, I, first(R).I, last(R).I)
end