text:
mutable struct HRU_Input
    #inputs (alphabetic order)
    Area_Elevations::Array{Float64,1}
    Area_HRU:: Float64
    Area_Glacier::Array{Float64,1} # smaller than 1
    Elevation_Count::Array{Int64}
    Nr_Elevationbands:: Int8
    Catchment_Elevation::Tuple
    Snow_Redistribution::Tuple
    #Potential_Evaporation::Array{Float64,1} #muss später auch Array werden!!! average Epot for soiL!!!
    Potential_Evaporation_Mean:: Float64
    Precipitation::Array{Float64,1}
    Riparian_Discharge:: Float64 #only necessary for riparian HRU
    Temp_Elevation::Array{Float64,1}
    Total_Effective_Precipitation::Float64
    Total_Interception_Evaporation::Float64
end

mutable struct Parameters
    # parameters (alphabetic order)
    beta:: Float64
    Ce:: Float64
    Drainagecapacity:: Float64 #only necessary for riparian HRU
    Interceptionstoragecapacity:: Float64
    Kf:: Float64
    Meltfactor:: Float64
    Mm:: Float64
    #Percolationcapacity:: Float64 #only necessary for hillslope HRU
    Ratio_Pref:: Float64 #only necessary for hillslope HRU
    Soilstoragecapacity:: Float64
    Temp_Thresh:: Float64
end

mutable struct Slow_Paramters
    Ks:: Float64
    Ratio_Riparian:: Float64
end

mutable struct Storages
    Fast:: Float64
    Interception::Array{Float64,1}
    Snow::Array{Float64,1}
    Snow_Cover::Array{Float64,1}
    Soil:: Float64
end

mutable struct Outflows
    Fast_Discharge:: Float64
    GWflow:: Float64 #only necessary for hillslope HRU
    Soil_Evaporation:: Float64
    Interception_Evaporation:: Float64
end

mutable struct Elevations
    Thickness_Band:: Float64
    Min_elevation:: Float64
    Max_elevation:: Float64
    Measured_Prec_Elevation:: Float64
    Measured_Temp_Elevation:: Float64
end

mutable struct Drought
    #inputs (alphabetic order)
    Nr_Drought_Days_Past::Array{Float64,1}
    Nr_Drought_Days_Future::Array{Float64,1}
    Nr_Drought_Events_Past::Array{Float64,1}
    Nr_Drought_Events_Future::Array{Float64,1}
    Max_Drought_Length_Past::Array{Float64,1}
    Max_Drought_Length_Future::Array{Float64,1}
    Mean_Drought_Length_Past::Array{Float64,1}
    Mean_Drought_Length_Future::Array{Float64,1}
    Max_Deficit_Past::Array{Float64,1}
    Max_Deficit_Future::Array{Float64,1}
    Mean_Deficit_Past::Array{Float64,1}
    Mean_Deficit_Future::Array{Float64,1}
    Total_Deficit_Past::Array{Float64,1}
    Total_Deficit_Future::Array{Float64,1}
    Max_Intensity_Past::Array{Float64,1}
    Max_Intensity_Future::Array{Float64,1}
    Mean_Intensity_Past::Array{Float64,1}
    Mean_Intensity_Future::Array{Float64,1}
end

mutable struct Drought_Extremes
    Longest_Drought_Length_Past::Array{Float64,1}
    Longest_Drought_Length_Future ::Array{Float64,1}
    Longest_Drought_Deficit_Past::Array{Float64,1}
    Longest_Drought_Deficit_Future::Array{Float64,1}
    Longest_Drought_Start_Past::Array{Float64,1}
    Longest_Drought_Start_Future::Array{Float64,1}
    Longest_Drought_End_Past::Array{Float64,1}
    Longest_Drought_End_Future::Array{Float64,1}
    Severest_Drought_Length_Past::Array{Float64,1}
    Severest_Drought_Length_Future ::Array{Float64,1}
    Severest_Drought_Deficit_Past::Array{Float64,1}
    Severest_Drought_Deficit_Future::Array{Float64,1}
    Severest_Drought_Start_Past::Array{Float64,1}
    Severest_Drought_Start_Future::Array{Float64,1}
    Severest_Drought_End_Past::Array{Float64,1}
    Severest_Drought_End_Future::Array{Float64,1}
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
using CompileBot

snoop_bench(
    BotConfig(
        "RecipesPipeline",
        yml_path= "SnoopCompile.yml",
        else_os = "linux",
        else_version = "1.5",
    )
)

ext: jl


text:
if basename(pwd()) == "aoc"
    cd("2015/12")
end

function part1(filename::AbstractString)
    sum(parse(Int, m.match) for m in eachmatch(r"-?\d+", readchomp(filename)))
end

using JSON
part2(filename::AbstractString) = stripred(JSON.parse(readchomp(filename)))

function stripred(d::AbstractDict)
    xs = values(d)
    if "red" in xs
        0
    else
        stripred(xs)
    end
end
stripred(xs::Union{AbstractArray,Base.ValueIterator}) = sum(stripred(x) for x in xs)
stripred(x::Integer) = x
stripred(::Any) = 0

@assert stripred(JSON.parse("[1,2,3]")) == 6
@assert stripred(JSON.parse("""[1,{"c":"red","b":2},3]""")) == 4
@assert stripred(JSON.parse("""{"d":"red","e":[1,2,3,4],"f":5}""")) == 0
@assert stripred(JSON.parse("""[1,"red",5]""")) == 6

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
struct SimData
    mechanism::Mechanism
    env::Environment
    state_cache::Vector{StateCache}
    envj_cache::Vector{EnvironmentJacobianCache}
    Δt::Real
    vs::VariableSelector
    cs::ConstraintSelector
    generate_solver_fn::Symbol
    extract_sol::Symbol
    lower_vs::Vector{VariableSelector}
    lower_cs::Vector{ConstraintSelector}
    lower_options::Vector{Dict}
    normal_vs::Vector{VariableSelector}
    normal_cs::Vector{ConstraintSelector}
    normal_options::Vector{Dict}
    fric_vs::Vector{VariableSelector}
    fric_cs::Vector{ConstraintSelector}
    fric_options::Vector{Dict}
    N::Int
    con_fns::Vector{Tuple{Symbol,Any}}
    obj_fns::Vector{Tuple{Symbol,Any}}
    fric_osqp_models::Vector{OSQP.Model}
    normal_osqp_models::Vector{OSQP.Model}
end

function add_eq!(sim_data::SimData, name::Symbol, size::Int, fun)
    add_eq!(sim_data.cs, name, size)
    push!(sim_data.con_fns, (name, fun))

    sim_data.con_fns
end

function add_ineq!(sim_data::SimData, name::Symbol, size::Int, fun)
    add_ineq!(sim_data.cs, name, size)
    push!(sim_data.con_fns, (name, fun))

    sim_data.con_fns
end

function add_obj!(sim_data::SimData, name::Symbol, fun)
    push!(sim_data.obj_fns, (name, fun))

    sim_data.obj_fns
end

function add_box_con!(sim_data::SimData, name::Symbol, var_name::Symbol, min::AbstractArray{T}, max::AbstractArray{T}) where T
    min_name = Symbol(name, "_min")
    add_ineq!(sim_data.cs, min_name, length(min))
    push!(sim_data.con_fns, (min_name, x -> min - sim_data.vs(x, var_name)))

    max_name = Symbol(name, "_max")
    add_ineq!(sim_data.cs, max_name, length(max))
    push!(sim_data.con_fns, (max_name, x -> sim_data.vs(x, var_name) - max))

    sim_data.con_fns
end

function simulate(sim_data::SimData,control!,state0::MechanismState,N::Int;
                  opt_tol=1e-6,major_feas=1e-6,minor_feas=1e-6,verbose=0)

    results = zeros(sim_data.vs.num_vars, N)
    results[sim_data.vs(:qnext), 1] = configuration(state0)
    results[sim_data.vs(:vnext), 1] = velocity(state0)

    x_ctrl = MechanismState(sim_data.mechanism)
    u0 = zeros(num_velocities(sim_data.mechanism))

    for i in 2:N
        x = results[:,i-1]
        q0 = sim_data.vs(x, :qnext)
        v0 = sim_data.vs(x, :vnext)

        set_configuration!(x_ctrl, q0)
        set_velocity!(x_ctrl, v0)
        control!(u0, (i-2)*sim_data.Δt, x_ctrl)

        solver_fn = eval(sim_data.generate_solver_fn)(sim_data, q0, v0, u0)

        options = Dict{String, Any}()
        options["Derivative option"] = 1
        options["Verify level"] = -1 # -1 => 0ff, 0 => cheap
        options["Major optimality tolerance"] = opt_tol
        options["Major feasibility tolerance"] = major_feas
        options["Minor feasibility tolerance"] = minor_feas

        xopt, info = snopt(solver_fn, sim_data.cs.num_eqs, sim_data.cs.num_ineqs, x, options)

        if verbose >= 1
            println(info)
        end

        results[:,i] .= xopt
    end

    sol = eval(sim_data.extract_sol)(sim_data, results)

    sol
end


text:
using LinearAlgebra: Tridiagonal

using Oceananigans:
    CPU, RegularRectilinearGrid, ZFaceField, ComputedField, set!, compute!, interior,
    FluxBoundaryCondition, GradientBoundaryCondition, TracerBoundaryConditions, fill_halo_regions!,
    Forcing, IncompressibleModel, Simulation, run!

using Oceananigans.Grids: Periodic, Bounded
using Oceananigans.OutputWriters: NetCDFOutputWriter, TimeInterval
using Oceananigans.AbstractOperations: @at, ∂z

function convective_adjustment!(model, Δt, K)
    Nz, Δz = model.grid.Nz, model.grid.Δz
    T = model.tracers.T

    ∂T∂z = ComputedField(@at (Center, Center, Center) ∂z(T))
    compute!(∂T∂z)

    κ = zeros(Nz)
    for i in 1:Nz
        κ[i] = ∂T∂z[1, 1, i] < 0 ? K : 0
    end

    ld = [-Δt/Δz^2 * κ[i]   for i in 2:Nz]
    ud = [-Δt/Δz^2 * κ[i+1] for i in 1:Nz-1]

    d = zeros(Nz)
    for i in 1:Nz-1
        d[i] = 1 + Δt/Δz^2 * (κ[i] + κ[i+1])
    end
    d[Nz] = 1 + Δt/Δz^2 * κ[Nz]

    𝓛 = Tridiagonal(ld, d, ud)

    T′ = 𝓛 \ interior(T)[:]
    set!(model, T=reshape(T′, (1, 1, Nz)))

    return nothing
end

function oceananigans_convective_adjustment_nn(ds; output_dir, nn_filepath, Δt=600)
    ρ₀ = 1027.0
    cₚ = 4000.0
    f  = ds.metadata["coriolis_parameter"]
    α  = ds.metadata["thermal_expansion_coefficient"]
    β  = 0.0
    g  = ds.metadata["gravitational_acceleration"]

    heat_flux = ds.metadata["temperature_flux"]
    ∂T₀∂z = ds.metadata["dθdz_deep"]

    T = ds["T"]
    wT = ds["wT"]
    Nz = size(T, 3)
    zc = znodes(T)
    zf = znodes(wT)
    Lz = abs(zf[1])

    Nt = size(T, 4)
    times = T.times
    stop_time = times[end]

    ## Grid setup

    topo = (Periodic, Periodic, Bounded)
    grid = RegularRectilinearGrid(topology=topo, size=(1, 1, Nz), extent=(1, 1, Lz))

    ## Boundary conditions

    T_bc_top = FluxBoundaryCondition(heat_flux)
    T_bc_bottom = GradientBoundaryCondition(∂T₀∂z)
    T_bcs = TracerBoundaryConditions(grid, top=T_bc_top, bottom=T_bc_bottom)

    ## Neural network forcing

    final_nn = jldopen(nn_filepath, "r")
    neural_network = final_nn["neural_network"]
    T_scaling = final_nn["T_scaling"]
    wT_scaling = final_nn["wT_scaling"]
    close(final_nn)

    μ_T, σ_T = T_scaling.μ, T_scaling.σ
    μ_wT, σ_wT = wT_scaling.μ, wT_scaling.σ

    function ∂z_wT(wT)
        wT_field = ZFaceField(CPU(), grid)
        set!(wT_field, reshape(wT, (1, 1, Nz+1)))
        fill_halo_regions!(wT_field, CPU(), nothing, nothing)
        ∂z_wT_field = ComputedField(@at (Center, Center, Center) ∂z(wT_field))
        compute!(∂z_wT_field)
        return interior(∂z_wT_field)[:]
    end

    enforce_fluxes(wT) = cat(0, wT, heat_flux, dims=1)

    # convective adjustment diffusivity
    K = wT_scaling.σ / T_scaling.σ * stop_time / Lz * 10

    function diagnose_wT_NN(model)
        T = interior(model.tracers.T)[:]
        T = T_scaling.(T)
        wT_NN_interior = neural_network(T)
        wT_NN_interior = inv(wT_scaling).(wT_NN_interior)
        wT_NN = enforce_fluxes(wT_NN_interior)

        ∂T∂z = ComputedField(@at (Center, Center, Face) ∂z(model.tracers.T))
        compute!(∂T∂z)

        κ = zeros(Nz+1)
        for i in 1:Nz+1
            κ[i] = ∂T∂z[1, 1, i] < 0 ? K : 0
        end

        K∂T∂z = κ .* interior(∂T∂z)[:]

        return wT_NN .- K∂T∂z
    end

    neural_network_forcing = Chain(
        T -> T_scaling.(T),
        neural_network,
        wT -> inv(wT_scaling).(wT),
        enforce_fluxes,
        ∂z_wT
    )

    ## TODO: Benchmark NN performance.

    ∂z_wT_NN = zeros(Nz)
    forcing_params = (∂z_wT_NN=∂z_wT_NN,)
    @inline neural_network_∂z_wT(i, j, k, grid, clock, model_fields, p) = - p.∂z_wT_NN[k]
    T_forcing = Forcing(neural_network_∂z_wT, discrete_form=true, parameters=forcing_params)

    ## Model setup

    model_convective_adjustment = IncompressibleModel(grid=grid, boundary_conditions=(T=T_bcs,))
    model_neural_network = IncompressibleModel(grid=grid, boundary_conditions=(T=T_bcs,), forcing=(T=T_forcing,))

    T₀ = reshape(Array(interior(T)[1, 1, :, 1]), size(grid)...)
    set!(model_convective_adjustment, T=T₀)
    set!(model_neural_network, T=T₀)

    ## Simulation setup

    function progress_convective_adjustment(simulation)
        clock = simulation.model.clock
        # @info "Convective adjustment: iteration = $(clock.iteration), time = $(prettytime(clock.time))"
        convective_adjustment!(simulation.model, simulation.Δt, K)
        return nothing
    end

    function progress_neural_network(simulation)
        model = simulation.model
        clock = simulation.model.clock

        # @info "Neural network: iteration = $(clock.iteration), time = $(prettytime(clock.time))"

        T = interior(model.tracers.T)[:]
        ∂z_wT_NN .= neural_network_forcing(T)

        convective_adjustment!(model, simulation.Δt, K)

        return nothing
    end

    simulation_convective_adjustment = Simulation(model_convective_adjustment, Δt=Δt, iteration_interval=1,
                                                  stop_time=stop_time, progress=progress_convective_adjustment)
    simulation_neural_network = Simulation(model_neural_network, Δt=Δt, iteration_interval=1,
                                           stop_time=stop_time, progress=progress_neural_network)

    ## Output writing

    filepath_CA = joinpath(output_dir, "oceananigans_convective_adjustment.nc")
    outputs_CA = (T  = model_convective_adjustment.tracers.T,)

    simulation_convective_adjustment.output_writers[:solution] =
        NetCDFOutputWriter(model_convective_adjustment, outputs_CA,
                           schedule = TimeInterval(Δt),
                           filepath = filepath_CA,
                           mode = "c")

    filepath_NN = joinpath(output_dir, "oceananigans_neural_network.nc")
    outputs_NN = (T  = model_neural_network.tracers.T,
                  wT = model -> diagnose_wT_NN(model))

    simulation_neural_network.output_writers[:solution] =
        NetCDFOutputWriter(model_neural_network, outputs_NN,
                           schedule = TimeInterval(Δt),
                           filepath = filepath_NN,
                           mode = "c",
                           dimensions = (wT=("zF",),))

    @info "Running convective adjustment simulation..."
    run!(simulation_convective_adjustment)

    @info "Running convective adjustment simulation + neural network..."
    run!(simulation_neural_network)

    ds_ca = NCDstack(filepath_CA)
    ds_nn = NCDstack(filepath_NN)

    T_ca = dropdims(Array(ds_ca[:T]), dims=(1, 2))
    T_nn = dropdims(Array(ds_nn[:T]), dims=(1, 2))
    wT_nn = Array(ds_nn[:wT])

    convective_adjustment_solution = (T=T_ca, wT=nothing)
    neural_network_solution = (T=T_nn, wT=wT_nn)

    return convective_adjustment_solution, neural_network_solution
end


text:
# Sieve of Eratosthenes based algorithm to find prime numbers up to a limit:
print("Enter an upper limit to find prime numbers upto, starting from one: ")
upperLimit = parse(Int, readline())
isPrime = [(i % 2) > 0 for i in 1:upperLimit]
isPrime[2] = true  # 2 is an exception, being the only even prime.
for num in 3:2:floor(Int, sqrt(upperLimit))
    isPrime[num] && foreach(i -> isPrime[i] = false, num * num:2 * num:upperLimit)
end
# Print all prime numbers from 1 to upperLimit:
foreach(i -> isPrime[i] && println(i), 2:upperLimit)

max_stars_repo_path: Utility Functions/Algorithms/SOE.jl

max_issues_repo_path: Utility Functions/Algorithms/SOE.jl

max_forks_repo_path: Utility Functions/Algorithms/SOE.jl


text:
using LaTeXStrings

struct TakeSnapshots{S <: AbstractVector{Array{Float64, 3}}} <: ProcessingMethod
    problem::FluidFlowProblem
    every_t::Union{Int, Vector{Int}}
    snapshots::S
    timesteps::Vector{Int}
end
TakeSnapshots(problem, every_t) =
    TakeSnapshots(problem, every_t, Array{Float64, 3}[], Int[])

function next!(pm::TakeSnapshots, q, f_in, t::Int64)
    if pm.every_t isa Int && mod(t, pm.every_t) != 0
        return false
    end
    if pm.every_t isa Vector{Int} && t ∉ pm.every_t
        return false
    end

    # snapshot = take_snapshot(pm, q, f_in, t)
    snapshot = copy(f_in)

    push!(pm.snapshots, snapshot)
    push!(pm.timesteps, t)

    visualize(pm, q)

    return false
end

function visualize(pm::TakeSnapshots, q::Quadrature)
    velocity_profile_x = plot(xlabel = L"x", ylabel = L"u_y")
    velocity_profile_y = plot(
        xlabel = "x",
        ylabel = L"u_y",
        legend = :bottomright,
        title = string("Velocity profile at ", latexstring("y = \\pi")),
    )
    pressure_profile = plot(xlabel = "x", ylabel = L"p")
    temperature_profile = plot(xlabel = "x", ylabel = L"T")
    sigma_xx_profile = plot(xlabel = "x", ylabel = L"\sigma_{xx}")
    sigma_xy_profile = plot(xlabel = "x", ylabel = L"\sigma_{xy}")

    problem = pm.problem
    for (timestep, f_in) in zip(pm.timesteps, pm.snapshots)
        time = delta_t(problem) * timestep
        # Pre-allocate Macroscopic Variables
        ρ = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
        u = Array{Float64}(undef, size(f_in, 1), size(f_in, 2), dimension(q))
        p = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
        T = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
        σ_xx = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))
        σ_xy = Array{Float64}(undef, size(f_in, 1), size(f_in, 2))

        Nx = size(f_in, 1)
        Ny = size(f_in, 2)

        x_range, y_range = range(problem)

        f = Array{Float64}(undef, size(f_in, 3))
        u_ = zeros(dimension(q))
        @inbounds for x_idx in 1:Nx, y_idx in 1:Ny
            @inbounds for f_idx in 1:size(f_in, 3)
                f[f_idx] = f_in[x_idx, y_idx, f_idx]
            end
            x = x_range[x_idx]
            y = y_range[y_idx]

            ρ_ = density(q, f)
            velocity!(q, f, ρ_, u_)
            T_ = temperature(q, f, ρ_, u_)
            p_ = pressure(q, f, ρ_, u_)

            τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem)
            σ_ = LatticeBoltzmann.deviatoric_tensor(q, τ, f, ρ_, u_)

            ρ_ = LatticeBoltzmann.dimensionless_density(problem, ρ_)
            u_ = LatticeBoltzmann.dimensionless_velocity(problem, u_)
            T_ = LatticeBoltzmann.dimensionless_temperature(q, problem, T_)
            p_ = LatticeBoltzmann.dimensionless_pressure(q, problem, p_)
            ρ[x_idx, y_idx] = ρ_
            u[x_idx, y_idx, :] = u_
            p[x_idx, y_idx] = p_
            T[x_idx, y_idx] = T_

            # u[x_idx, y_idx, :] .-= LatticeBoltzmann.velocity(problem, x, y, time)
            σ_ = LatticeBoltzmann.dimensionless_stress(problem, σ_)
            σ_xx[x_idx, y_idx] = σ_[1, 1]
            σ_xy[x_idx, y_idx] = σ_[1, 2]
        end

        x_pos = max(round(Int, problem.NX / 2), 1)
        y_pos = max(round(Int, problem.NY / 2), 1)
        domain = x_range[1:Nx]

        x_range, y_range = range(problem)
        velocity = (x, y, t) -> LatticeBoltzmann.velocity(problem, x, y, t)
        σ = (x, y, t) -> LatticeBoltzmann.deviatoric_tensor(q, problem, x, y, t)
        pr = (x, y, t) -> LatticeBoltzmann.pressure(q, problem, x, y, t)

        exact_range = range(0.0, length = 1000, stop = problem.domain_size[1])

        plot!(
            velocity_profile_x,
            domain,
            u[:, y_pos, 2],
            label = latexstring("t = ", time),
            linecolor = :gray,
        )
        y_exact = map(exact_range) do y
            velocity(y, y_range[y_pos], time)[2]
        end
        # plot!(velocity_profile_x, exact_range, y_exact,  label = "", linecolor = :gray, linealpha = 0.2, linestyle = :dash)

        # annotation_x = 0.5
        # annotation_y = u[x_pos, end - 3, 1]
        # annotation_y = u[x_pos, round(Int, problem.NY /2) + 1, 1] + 0.005
        # annotate!(velocity_profile_x, annotation_y, annotation_x, text(L"\nu / t = 0.5", rotation = 00))

        # plot!(velocity_profile_x, domain, u[x_pos, :, 1], label = latexstring("t = ", time), linecolor = :gray)
        # plot!(velocity_profile_x, (y) -> velocity(x_range[x_pos], y, time)[1], exact_range, label = "", linecolor = :gray, linealpha = 0.2, linestyle = :dash)

        # annotation_x = 0.5
        # annotation_y = u[x_pos, end - 3, 1]
        # annotation_y = u[x_pos, round(Int, problem.NY /2) + 1, 1] + 0.005
        # annotate!(velocity_profile_x, annotation_x, annotation_y, text(L"\nu / t = 0.5", rotation = 00))

        # annotate!(velocity_profile_x, [(5, y[5], Plots.text("this is #5", 16, :red, :center)), (10, y[10], Plots.text("this is #10", :right, 20, "courier"))])

        # scatter!(velocity_profile_x, domain, u[x_pos, :, 2], label = latexstring("t = ", timestep))

        scatter!(
            velocity_profile_y,
            domain,
            u[x_pos, :, 2],
            label = latexstring("t = ", timestep),
            markershape = :auto,
            markersize = 6,
        )
        plot!(
            (y) -> velocity(x_range[x_pos], y, time)[2],
            exact_range,
            label = "",
            linecolor = :gray,
            linealpha = 0.2,
            linestyle = :dash,
        )

        scatter!(
            pressure_profile,
            domain,
            p[x_pos, :],
            label = latexstring("t = ", timestep),
        )
        plot!(
            pressure_profile,
            (y) -> pr(x_range[x_pos], y, time),
            exact_range,
            label = "",
            linecolor = :gray,
            linealpha = 0.2,
            linestyle = :dash,
        )

        plot!(
            temperature_profile,
            domain,
            T[x_pos, :],
            label = latexstring("t = ", timestep),
        )

        scatter!(
            sigma_xx_profile,
            domain,
            σ_xx[x_pos, :],
            label = latexstring("t = ", timestep),
        )
        plot!(
            sigma_xx_profile,
            (y) -> σ(x_range[x_pos], y, time)[1, 1],
            exact_range,
            label = "",
            linecolor = :gray,
            linealpha = 0.2,
            linestyle = :dash,
        )

        scatter!(
            sigma_xy_profile,
            domain,
            σ_xy[x_pos, :],
            label = latexstring("t = ", timestep),
        )
        plot!(
            sigma_xy_profile,
            (y) -> σ(x_range[x_pos], y, time)[1, 2],
            exact_range,
            label = "",
            linecolor = :gray,
            linealpha = 0.2,
            linestyle = :dash,
        )
    end

    ps = (
        velocity_profile_x = velocity_profile_x,
        velocity_profile_y = velocity_profile_y,
        temperature_profile = temperature_profile,
        pressure_profile = pressure_profile,
        sigma_xx_profile = sigma_xx_profile,
        sigma_xy_profile = sigma_xy_profile,
    )

    return ps
end


text:
# Matern 3/2 isotropic covariance function

"""
    Mat32Iso <: MaternIso

Isotropic Matern 3/2 kernel (covariance)
```math
k(x,x') = σ²(1 + √3|x-x'|/ℓ)\\exp(-√3|x-x'|/ℓ)
```
with length scale ``ℓ`` and signal standard deviation ``σ``.
"""
mutable struct Mat32Iso <: MaternIso
    "Length scale"
    ℓ::Float64
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Mat32Iso(ll::Float64, lσ::Float64)

    Create `Mat32Iso` with length scale `exp(ll)` and signal standard deviation `exp(lσ)`.
    """
    Mat32Iso(ll::Float64, lσ::Float64) = new(exp(ll), exp(2 * lσ), [])
end

function set_params!(mat::Mat32Iso, hyp::VecF64)
    length(hyp) == 2 || throw(ArgumentError("Matern 3/2 only has two parameters"))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2 * hyp[2])
end

get_params(mat::Mat32Iso) = Float64[log(mat.ℓ), log(mat.σ2) / 2 ]
get_param_names(mat::Mat32Iso) = [:ll, :lσ]
num_params(mat::Mat32Iso) = 2

Statistics.cov(mat::Mat32Iso, r::Float64) =
    (s = √3 * r / mat.ℓ; mat.σ2 * (1 + s) * exp(-s))

@inline dk_dll(mat::Mat32Iso, r::Float64) =
    (s = √3 * r / mat.ℓ; mat.σ2 * s^2 * exp(-s))


text:
using Images

using PCG.Types
using PCG.Geometry
using PCG.Universes
using PCG.Topologies.SquareGreedTopologies
using PCG.Neighborhoods
using PCG.Operations
using PCG.Spaces
using PCG.Recorders.GreedImages
using PCG.Recorders.SquareGreedImages
using PCG.Recorders.TurnsLogger

const WIDTH = 80
const HEIGHT = 80
const CELL_SIZE = Size(10, 10)
const DEBUG = false

struct Terrain <: Checkable
    value::Int64
end

function Operations.check(element::Element, parameters::Terrain)
    return element.properties.terrain == parameters
end

const GRASS = Terrain(1)
const WATER = Terrain(2)
const SAND = Terrain(3)
const FOREST = Terrain(4)

struct Properties <: AbstractProperties
    terrain::Terrain
end

const SPRITE_GRASS = SquareSprite(RGB(0, 1, 0), CELL_SIZE)
const SPRITE_WATER = SquareSprite(RGB(0, 0, 1), CELL_SIZE)
const SPRITE_SAND = SquareSprite(RGB(1, 1, 0), CELL_SIZE)
const SPRITE_FOREST = SquareSprite(RGB(0, 0.5, 0), CELL_SIZE)
const SPRITE_ERROR = SquareSprite(RGB(1, 0, 0), CELL_SIZE)

function GreedImages.choose_sprite(recorder::GreedImageRecorder, element)
    if element |> GRASS |> exists
        return SPRITE_GRASS
    end

    if element |> WATER |> exists
        return SPRITE_WATER
    end

    if element |> SAND |> exists
        return SPRITE_SAND
    end

    if element |> FOREST |> exists
        return SPRITE_FOREST
    end

    return SPRITE_ERROR
end

function process(universe::Universe)

    ring = Neighborhood(universe.topology, ring_distance)
    euclidean = Neighborhood(universe.topology, euclidean_distance)

    complete_turn!(universe)

    universe() do element
        if element |> Fraction(0.01) |> exists
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> Fraction(0.8) |> GRASS |> euclidean(1, 3) |> WATER |> exists
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> GRASS |> ring() |> WATER |> exists
            element << (terrain=SAND,)
        end
    end

    universe(turns=3) do element
        if element |> Fraction(0.1) |> GRASS |> ring() |> SAND |> exists
            element << (terrain=SAND,)
        end
    end

    universe(turns=4) do element
        if element |> SAND |> ring() |> WATER |> count >= 6
            element << (terrain=WATER,)
        end
    end

    universe() do element
        if element |> GRASS |> Fraction(0.03) |> exists
            element << (terrain=FOREST,)
        end
    end

    universe() do element
        if (element |> GRASS |> Fraction(0.1) |> exists &&
            element |> ring(2, 2) |> FOREST |> exists &&
            element |> ring() |> new |> FOREST |> not_exists)

            element << (terrain=FOREST,)
        end
    end

end

topology = SquareGreedTopology(WIDTH, HEIGHT)

universe = initialize(topology,
                      Properties(GRASS),
                      DEBUG ? Recorder[] : [TurnsLoggerRecorder(),
                                            GreedImageRecorder(CELL_SIZE, 1000, "output.gif")])

@time process(universe)

finish_recording!(universe)

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
"""
    generate_ta_func_with_arrays(d, symb)

Generate function code (level 1) with arrays as input / output

# Examples
```julia
julia> generate_ta_func_with_arrays(d, :MAC)
"..."
```
"""
function generate_ta_func_with_arrays(d::OrderedDict{Symbol,Any}, symb_func::Symbol)

    func_info = d[symb_func]

    s_doc_RequiredInputArguments = ""
    s_doc_OptionalInputArguments = ""
    s_doc_OutputArguments = ""

    params_lv0 = String[]
    params_lv1 = String[]
    params_RequiredInputArguments = String[]
    params_OptionalInputArguments = String[]
    params_OutputArguments = String[]

    jltypes = String[]

    for arg = ["startIdx", "endIdx"]
        varname = uncamel(arg)
        push!(params_lv0, varname)
    end

    for arg = func_info["RequiredInputArguments"]
        varname = replace_var(arg["Name"])
        vartyp = d_typ_to_jl[arg["Type"]]
        push!(params_lv0, string(varname))
        push!(params_lv1, string(varname))
        push!(params_RequiredInputArguments, string(varname))
        push!(jltypes, string(vartyp))
        s_doc_RequiredInputArguments *= "\n        - $varname::$vartyp"
    end

    for arg = func_info["OptionalInputArguments"]
        varname = fix_varname(arg["Name"])
        vartyp = d_typ_to_jl[arg["Type"]]
        def_val = arg["DefaultValue"]
        push!(params_lv0, varname)
        push!(params_lv1, varname)
        push!(params_OptionalInputArguments, varname)
        push!(jltypes, string(vartyp))
        s_doc_OptionalInputArguments *= "\n        - $varname=$vartyp($def_val)"
    end

    for arg = ["outBegIdx", "outNbElement"]
        varname = arg #uncamel(arg)
        push!(params_lv0, varname)
        push!(jltypes, "Ref{Cint}")
    end

    for arg = func_info["OutputArguments"]
        varname = arg["Name"]
        vartyp = d_typ_to_jl[arg["Type"]]
        push!(params_lv0, varname)
        push!(params_OutputArguments, varname)
        push!(jltypes, string(vartyp))
        s_doc_OutputArguments *= "\n        - $varname::$vartyp"
    end

    params_lv0 = join(params_lv0, ", ")
    params_lv1 = join(params_lv1, ", ")

    params_lv1_with_types = ""
    for (index, arg) in enumerate(func_info["RequiredInputArguments"])
        if index != 1
            params_lv1_with_types *= ", "
        end
        varname = replace_var(arg["Name"])
        vartyp = d_typ_to_jl[arg["Type"]]
        params_lv1_with_types *= "$varname::$vartyp"
    end
    for (index, arg) in enumerate(func_info["OptionalInputArguments"])
        if index != 1
            params_lv1_with_types *= ", "
        else
            params_lv1_with_types *= "; "
        end
        varname = fix_varname(arg["Name"])
        vartyp = d_typ_to_jl[arg["Type"]]
        def_val = arg["DefaultValue"]
        params_lv1_with_types *= "$varname=$vartyp($def_val)"
    end

    jltypes = join(jltypes, ", ")

    funcname = string(symb_func)
    first_req_input_array = replace_var(func_info["RequiredInputArguments"][1]["Name"])

    s = "

\"\"\"
    $funcname($params_lv1_with_types)

$(func_info["ShortDescription"]) ($(func_info["CamelCaseName"]))

    $(func_info["GroupId"])

    Level: 1 - Arrays

Arguments:

    RequiredInputArguments:$(s_doc_RequiredInputArguments)

    OptionalInputArguments:$(s_doc_OptionalInputArguments)

Returns:$(s_doc_OutputArguments)

\"\"\"
function $funcname($params_lv1_with_types)
    N = length($first_req_input_array)
    start_idx = 0
    end_idx = N - 1
    outBegIdx = Ref{Cint}(0)
    outNbElement = Ref{Cint}(0)"

for arg = params_OutputArguments
    varname = arg
    s *= "\n"
    s *= INDENT * "$(varname) = fill(NaN, N)"
end

s *= "\n"
s *= INDENT * "ret_code = _TA_$funcname($params_lv0)"
s *= "\n"
s *= INDENT * "_ta_check_success(\"$funcname\", ret_code)"

for arg = params_OutputArguments
    s *= "\n"
    s *= INDENT * "$arg = circshift($arg, outBegIdx[])"
end

s *= "\n"
if length(params_OutputArguments) == 1
    s *= INDENT * params_OutputArguments[1]
else
    s *= INDENT * "[" * join(params_OutputArguments, " ") * "]"  # output as Array{Any,2}
    #s *= INDENT * join(params_OutputArguments, ", ")  # output as tuple of array
end
s *= "\n"
s *= "end"

s

end


text:
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

using Printf

"""
Approximate |z*exp(1-z)|==1 with z = r(ϕ)*exp(iϕ) for 101 values in [0,π]
```
Piecewise[{{1/E, 2*ϕ == Pi || 2*ϕ == 3*Pi}},
  (-ProductLog[(Abs[Cos[ϕ]]*
        Piecewise[{{1, Pi/2 <= ϕ <= (3*Pi)/2}}, -1])/E])* Sec[ϕ]]
```
"""
const curve_r_of_phi = (
    1.00000000000000000, 0.96938799932339478, 0.94031280913960984,
    0.91267447197920273, 0.88638114186065965, 0.86134829907186178,
    0.83749805315932746, 0.81475852287976240, 0.79306328346805633,
    0.77235087292502438, 0.75256435016959500, 0.73365089886863247,
    0.71556147158168498, 0.69825046956111954, 0.68167545414968103,
    0.66579688623349311, 0.65057789065218041, 0.63598404285018716,
    0.62198317538372242, 0.60854520218380338, 0.59564195872407880,
    0.58324705645795586, 0.57133575007763025, 0.55988481631183171,
    0.54887244312274892, 0.53827812828850114, 0.52808258646807237,
    0.51826766394286338, 0.50881626031469696, 0.49971225651573406,
    0.49094044855261233, 0.48248648646631564, 0.47433681804178095,
    0.46647863684787701, 0.45889983422986405, 0.45158895491338878,
    0.44453515591202428, 0.43772816845980179, 0.43115826271651143,
    0.42481621501712971, 0.41869327745787858, 0.41278114963040691,
    0.40707195233265375, 0.40155820310031468, 0.39623279341667475,
    0.39108896747105881, 0.38612030234742765, 0.38132068953484406,
    0.37668431766076138, 0.37220565635644773, 0.36787944117144232,
    0.36370065946082444, 0.35966453717533229, 0.35576652649006006,
    0.35200229421264355, 0.34836771091656781, 0.34485884074953856,
    0.34147193187079335, 0.33820340747482355, 0.33504985736226709,
    0.33200803002174285, 0.32907482518915661, 0.32624728685353880,
    0.32352259668079724, 0.32089806782890166, 0.31837113912997935,
    0.31593936961660664, 0.31360043337124438, 0.31135211467929895,
    0.30919230346770490, 0.30711899101223196, 0.30513026589792689,
    0.30322431021821778, 0.30139939599924305, 0.29965388183692661,
    0.29798620973521020, 0.29639490213468041, 0.29487855912159646,
    0.29343585580803966, 0.29206553987457226, 0.29076642926741500,
    0.28953741004273363, 0.28837743435116807, 0.28728551855624701,
    0.28626074148080829, 0.28530224277599438, 0.28440922140781493,
    0.28358093425666686, 0.28281669482557941, 0.28211587205330827,
    0.28147788922874187, 0.28090222300340517, 0.28038840249915426,
    0.27993600850844928, 0.27954467278487599, 0.27921407742185853,
    0.27894395431776862, 0.27873408472589181, 0.27858429888795920,
    0.27849447575019571, 0.27846454276107380)

function get_curve()
    ϕ1 = LinRange(0, π, 101)
    ϕ2 = -ϕ1[end-1:-1:2]
    r1, r2 = curve_r_of_phi, curve_r_of_phi[end-1:-1:2]
    line_seg = [collect(r1).*exp.(1im*collect(ϕ1))...,
                collect(r2).*exp.(1im*collect(ϕ2))...]
    return CV_LinePainter([line_seg], true)
end

const fontface = cv_fontface("sans-serif")

function show_comparison(n, save=false)
    param_n = CV_TranslateByOffset(Int)
    param_n.value = n

    function trafo2(z)
        n = param_n.value
        z = z * n
        result = big(1.0)
        for k in 1:n
          result += big(z)^k/factorial(big(k))
        end
        return ComplexF64(result)
    end

    trafo1 = z -> exp(z*param_n.value)

    codomain1 = CV_Math2DCanvas(-0.5 + 1.0im, 1.5 - 1.0im, 200)
    codomain2 = CV_Math2DCanvas(-0.5 + 1.0im, 1.5 - 1.0im, 200)

    label_style = cv_black → fontface → cv_fontsize(20)
    nstyle = cv_black → fontface → cv_fontsize(40)
    nrulers=(
        CV_TickLabelAppearance(; label_style) ↦ ("%.0f" ⇒ 0.0:1.0),
        CV_TickLabelAppearance(; label_style=nstyle, tick_length=0) ↦
            (("n = " * string(n)) ⇒ 1.2,))
    rulers=(
        CV_TickLabelAppearance(; label_style) ↦ ("%.0f" ⇒ 0.0:1.0),)

    curve = (cv_black → cv_linewidth(3)) ↦ get_curve()

    scene = cv_scene_comp_codomains_std(tuple(), trafo1, trafo2,
        codomain1, codomain2; codomain1_re_rulers=rulers,
        codomain1_im_rulers=rulers,
        codomain2_re_rulers=rulers,
        codomain2_im_rulers=nrulers,
        painter1=CV_PortraitPainter(trafo1) → curve,
        painter2=CV_PortraitPainter(trafo2) → curve)
    cv_get_redraw_func(scene)()

    cv_save_image(cv_get_can_layout(scene),
        @sprintf("ExpSeriesRescaled_%03i.png", n))
    cv_destroy(scene)
end

for n in (1, 5, 10, 30, 60, 80, 100)
    println("n = ", n)
    show_comparison(n)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
###############################
#
# FGL from Danaher et al.
#
# currently implements only two groups
###############################

function fusedGraphicalLasso!(
    Θx::StridedMatrix{T}, Θy::StridedMatrix{T},
    Zx::StridedMatrix{T}, Zy::StridedMatrix{T},
    Ux::StridedMatrix{T}, Uy::StridedMatrix{T},
    Σx::StridedMatrix{T}, nx::Int,
    Σy::StridedMatrix{T}, ny::Int,
    λ1::Real, λ2::Real;
    options::ADMMOptions = ADMMOptions(),
    penalize_diag::Bool=true) where {T<:AbstractFloat}

    maxiter = options.maxiter
    ρ = options.ρ
    α = options.α
    abstol = options.abstol
    reltol = options.reltol

    n = nx + ny
    γ = 1. / ρ
    γx = nx * γ / n
    γy = ny * γ / n

    p = size(Σx, 1)

    gx = ProxGaussLikelihood(Symmetric(Σx))
    gy = ProxGaussLikelihood(Symmetric(Σy))

    Tx = zeros(p, p)  # temp storage
    Ty = zeros(p, p)  # temp storage

    Zxold = copy(Zx)
    Zyold = copy(Zy)

    # ADMM
    for it=1:maxiter
        # update Θx and Θy
        @. Tx = Zx - Ux
        prox!(gx, Θx, Tx, γx)
        @. Ty = Zy - Uy
        prox!(gy, Θy, Ty, γy)

        copyto!(Zxold, Zx)
        copyto!(Zyold, Zy)
        # update Zx and Zy with relaxation
        # α⋅X + (1-α)⋅Z
        @. Tx = α * Θx + (1. - α)*Zx + Ux
        @. Ty = α * Θy + (1. - α)*Zy + Uy
        for i=1:p
            for j=i:p
                if i==j
                    z1, z2 = proxL1Fused(Tx[i,i], Ty[i,i], penalize_diag ? λ1 : 0., λ2, γ)
                    Zx[i,i] = z1
                    Zy[i,i] = z2
                else
                    z1, z2 = proxL1Fused(Tx[i,j], Ty[i,j], λ1, λ2, γ)
                    Zx[i,j] = Zx[j,i] = z1
                    Zy[i,j] = Zy[j,i] = z2
                end
            end
        end

        # update U
        @. Ux = Tx - Zx
        @. Uy = Ty - Zy

        # check convergence
        r_norm = norm_diff(Θx, Zx, 2) + norm_diff(Θy, Zy, 2)
        s_norm = (norm_diff(Zx, Zxold, 2) + norm_diff(Zy, Zyold, 2)) * sqrt(ρ)
        eps_pri = p*abstol + reltol * max( norm(Θx) + norm(Θy), norm(Zx) + norm(Zy))
        eps_dual = p*abstol + reltol * ρ * ( norm(Ux) + norm(Uy) )
        if r_norm < eps_pri && s_norm < eps_dual
            break
        end
    end
    (Zx, Zy)
end

function fusedGraphicalLasso(
    Σx::StridedMatrix{T}, nx::Int,
    Σy::StridedMatrix{T}, ny::Int,
    λ1::Real, λ2::Real;
    options::ADMMOptions = ADMMOptions(),
    penalize_diag::Bool=true) where {T<:AbstractFloat}

    p = size(Σx, 1)

    Θx = eye(p)
    Θy = eye(p)
    Zx = zeros(p, p)
    Zy = zeros(p, p)
    Ux = zeros(p, p)
    Uy = zeros(p, p)

    fusedGraphicalLasso!(
    Θx, Θy,
    Zx, Zy,
    Ux, Uy,
    Σx, nx, Σy, ny,
    λ1, λ2; options=options, penalize_diag=penalize_diag)
end

####################################
#
# Fused neighborhood selection
#
# ADMM implementation of NIPS paper
####################################

function fusedNeighborhoodSelection(
    A1, b1,
    A2, b2,
    λ1, λ2;
    options::ADMMOptions = ADMMOptions()
    )

    @assert size(A1, 2) == size(A2, 2)
    @assert size(A1, 1) == length(b1)
    @assert size(A2, 1) == length(b2)

    maxiter = options.maxiter
    ρ = options.ρ
    α = options.α
    abstol = options.abstol
    reltol = options.reltol
    γ = 1. / ρ

    n1 = size(A1, 1)
    n2 = size(A2, 1)
    p = size(A1, 2)

    AtA1 = A1'A1 / n1
    AtA2 = A2'A2 / n2
    for i=1:p
        AtA1[i,i] += ρ
        AtA2[i,i] += ρ
    end
    cAtA1 = cholesky!(AtA1, Val(false))
    cAtA2 = cholesky!(AtA2, Val(false))

    Atb1 = A1'*b1 / n1
    Atb2 = A2'*b2 / n2

    x1 = zeros(p)
    x2 = zeros(p)
    z1 = zeros(p)
    z2 = zeros(p)
    u1 = zeros(p)
    u2 = zeros(p)

    T1 = zeros(p)  # temp storage
    T2 = zeros(p)

    z1old = copy(z1)
    z2old = copy(z2)

    # ADMM
    for it=1:maxiter
        # update x1 and x2
        @. T1 = Atb1 + (z1 - u1) * ρ
        x1 = cAtA1 \ T1
        @. T2 = Atb2 + (z2 - u2) * ρ
        x2 = cAtA2 \ T2

        copyto!(z1old, z1)
        copyto!(z2old, z2)
        # update z1 and z2
        @. T1 = α*x1 + (1. - α)*z1 + u1
        @. T2 = α*x2 + (1. - α)*z2 + u2
        for i=1:p
            z1[i], z2[i] = proxL1Fused(T1[i], T2[i], λ1, λ2, γ)
        end

        # update U
        @. u1 = T1 - z1
        @. u2 = T2 - z2

        # check convergence
        r_norm = norm_diff(x1, z1, 2) + norm_diff(x2, z2, 2)
        s_norm = (norm_diff(z1, z1old, 2) + norm_diff(z2, z2old, 2)) * sqrt(ρ)
        eps_pri = sqrt(p)*abstol + reltol * max( norm(x1) + norm(x2), norm(z1) + norm(z2))
        eps_dual = sqrt(p)*abstol + reltol * ρ *( norm(u1) + norm(u2) )
        if r_norm < eps_pri && s_norm < eps_dual
            break
        end
    end
    (z1, z2)
end


text:
export SaintVenant

"""
    SaintVenant(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
Saint-Venant (or shallow water) model.

# Argument
`param` is of type `NamedTuple` and must contain
- the dimensionless parameter `ϵ` (nonlinearity);
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Saint-Venant"`);

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `SaintVenant.f!` to be called in explicit time-integration solvers;
2. a function `SaintVenant.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `SaintVenant.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct SaintVenant <: AbstractModel

	label 	:: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info	:: String

    function SaintVenant(param::NamedTuple;
						mesh = Mesh(param),
						dealias=0,ktol=0,
						label="Saint-Venant"
						)

		m=WhithamBoussinesq(merge(param,(μ=1,));Boussinesq=true,
							mesh=mesh,
							a=0,b=0,
							dealias=dealias,ktol=ktol,
							label=label)

		info = "Saint-Venant model.\n"
		info *= "├─Nonlinearity parameter ϵ=$(param.ϵ).\n"
		if dealias == 0
			info *= "└─No dealiasing. "
		else
			info *= "└─Dealiasing with Orszag's rule adapted to power $(dealias + 1) nonlinearity. "
		end
		if ktol == 0
			info *= "No Krasny filter. "
		else
			info *= "Krasny filter with tolerance $ktol."
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		new(m.label, m.f!, m.mapto, m.mapfro, info)
    end
end


text:
"""
    AssemblyModule

Module for assemblers  of system matrices and vectors.
"""
module AssemblyModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import SparseArrays: sparse
import LinearAlgebra: diag

"""
    AbstractSysmatAssembler

Abstract type of system-matrix assembler.
"""
abstract type AbstractSysmatAssembler end;

"""
    SysmatAssemblerSparse{T<:Number} <: AbstractSysmatAssembler

Type for assembling a sparse global matrix from elementwise matrices.
"""
mutable struct SysmatAssemblerSparse{T<:Number} <: AbstractSysmatAssembler
    buffer_length::FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer::FInt;
    ndofs_row::FInt; ndofs_col::FInt;
end

"""
    SysmatAssemblerSparse(zero::T=0.0) where {T<:Number}

Construct blank system matrix assembler. The matrix entries are of type `T`.

# Example

This is how a sparse matrix is assembled from two rectangular dense matrices.
```
	a = SysmatAssemblerSparse(0.0)
	startassembly!(a, 5, 5, 3, 7, 7)
	m = [0.24406   0.599773    0.833404  0.0420141
		0.786024  0.00206713  0.995379  0.780298
		0.845816  0.198459    0.355149  0.224996]
	assemble!(a, m, [1 7 5], [5 2 1 4])
	m = [0.146618  0.53471   0.614342    0.737833
		 0.479719  0.41354   0.00760941  0.836455
		 0.254868  0.476189  0.460794    0.00919633
		 0.159064  0.261821  0.317078    0.77646
		 0.643538  0.429817  0.59788     0.958909]
	assemble!(a, m, [2 3 1 7 5], [6 7 3 4])
	A = makematrix!(a)
```
"""
function SysmatAssemblerSparse(zero::T=0.0) where {T<:Number}
    return SysmatAssemblerSparse{T}(0,[zero],[0],[0],0,0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparse{T},
      elem_mat_nrows::FInt, elem_mat_ncols::FInt, elem_mat_nmatrices::FInt,
      ndofs_row::FInt, ndofs_col::FInt) where {T<:Number}

Start the assembly of a global matrix.
"""
function startassembly!(self::SysmatAssemblerSparse{T},
  elem_mat_nrows::FInt, elem_mat_ncols::FInt, elem_mat_nmatrices::FInt,
  ndofs_row::FInt, ndofs_col::FInt) where {T<:Number}
    # Start the assembly of a global matrix.
    # The method makes buffers for matrix assembly. It must be called before
    # the first call to the method assemble.
    # elem_mat_nrows= number of rows in typical element matrix,
    # elem_mat_ncols= number of columns in a typical element matrix,
    # elem_mat_nmatrices= number of element matrices,
    # ndofs_row= Total number of equations in the row direction,
    # ndofs_col= Total number of equations in the column direction.
    self.buffer_length = elem_mat_nmatrices*elem_mat_nrows*elem_mat_ncols;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs_row = ndofs_row;
    self.ndofs_col = ndofs_col;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T},
      dofnums_row::FIntMat, dofnums_col::FIntMat) where {T<:Number}

Assemble a rectangular matrix.
"""
function assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T}, dofnums_row::FIntVec, dofnums_col::FIntVec) where {T<:Number}
    # Assembly of a rectangular matrix.
    # The method assembles a rectangular matrix using the two vectors of
    # equation numbers for the rows and columns.
    nrows=length(dofnums_row); ncolumns=length(dofnums_col);
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @inbounds for j=1:ncolumns
        @inbounds for i=1:nrows
            self.matbuffer[p] = mat[i,j] # serialized matrix
            self.rowbuffer[p] = dofnums_row[i];
            self.colbuffer[p] = dofnums_col[j];
            p=p+1
        end
    end
    self.buffer_pointer=p;
    return self
end

function assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T}, dofnums_row::FIntMat, dofnums_col::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums_row), vec(dofnums_col))
end

"""
    makematrix!(self::SysmatAssemblerSparse)

Make a sparse matrix.
"""
function makematrix!(self::SysmatAssemblerSparse)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs_row) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j]=self.ndofs_row+1;
        end
        if (self.colbuffer[j] > self.ndofs_col) || (self.colbuffer[j] <= 0)
            self.colbuffer[j]=self.ndofs_col+1;
        end
    end
    S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs_row+1, self.ndofs_col+1);
    self=SysmatAssemblerSparse(0.0*self.matbuffer[1])# get rid of the buffers
    return S[1:end-1,1:end-1]
end

"""
    SysmatAssemblerSparseSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a **symmetric square** matrix  assembled from symmetric square
matrices.
"""
mutable struct SysmatAssemblerSparseSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    buffer_length:: FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer:: FInt;
    ndofs:: FInt;
end

"""
    SysmatAssemblerSparseSymm(zero::T=0.0) where {T<:Number}

Construct blank system matrix assembler for symmetric matrices. The matrix
entries are of type `T`.

# Example

This is how a symmetric sparse matrix is assembled from two square dense matrices.
```
	a = SysmatAssemblerSparseSymm(0.0)
	startassembly!(a, 5, 5, 3, 7, 7)
	m = [0.24406   0.599773    0.833404  0.0420141
		0.786024  0.00206713  0.995379  0.780298
		0.845816  0.198459    0.355149  0.224996]
	assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
	m = [0.146618  0.53471   0.614342    0.737833
		 0.479719  0.41354   0.00760941  0.836455
		 0.254868  0.476189  0.460794    0.00919633
		 0.159064  0.261821  0.317078    0.77646
		 0.643538  0.429817  0.59788     0.958909]
	assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
	A = makematrix!(a)
```
# See also
[]
"""
function SysmatAssemblerSparseSymm(zero::T=0.0) where {T<:Number}
    return SysmatAssemblerSparseSymm{T}(0,[zero],[0],[0],0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparseSymm{T},
      elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
      ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric square global matrix.
"""
function startassembly!(self::SysmatAssemblerSparseSymm{T},
  elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
  ndofs::FInt, ignore2::FInt) where {T<:Number}
    # Start the assembly of a global matrix.
    # The method makes buffers for matrix assembly. It must be called before
    # the first call to the method assemble.
    # elem_mat_nrows= number of rows in typical element matrix,
    # elem_mat_ncols= number of columns in a typical element matrix,
    # elem_mat_nmatrices= number of element matrices,
    # ndofs_row= Total number of equations in the row direction,
    # ndofs_col= Total number of equations in the column direction.
    self.buffer_length = elem_mat_nmatrices*elem_mat_dim^2;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are the same).
"""
function assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}
    # Assembly of a square symmetric matrix.
    # The method assembles the lower triangle of the square symmetric matrix using the two vectors of
    # equation numbers for the rows and columns.
    nrows=length(dofnums); ncolumns=nrows;
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @inbounds for j=1:ncolumns
        @inbounds for i=j:nrows
            self.matbuffer[p] = mat[i,j] # serialized matrix
            self.rowbuffer[p] = dofnums[i];
            self.colbuffer[p] = dofnums[j];
            p=p+1
        end
    end
    self.buffer_pointer=p;
    return self
end

"""
	assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},
	  dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are the same).
"""
function assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T}, dofnums::FIntMat, ignore::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums), vec(ignore))
end

"""
    makematrix!(self::SysmatAssemblerSparseSymm)

Make a sparse symmetric square matrix.
"""
function makematrix!(self::SysmatAssemblerSparseSymm)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j]=self.ndofs+1;
        end
        if (self.colbuffer[j] > self.ndofs) || (self.colbuffer[j] <= 0)
            self.colbuffer[j]=self.ndofs+1;
        end
    end
    S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs+1, self.ndofs+1);
    S = S+copy(transpose(S)); # This is to address the lack of a function for adding together a sparse matrix with a transpose of the sparse matrix; at the moment (January 2018), this defaults to the addition of two dense matrices and the result is DENSE. Hence a copy of the transpose needs to be made.
    #   S = S+transpose(S);    # construct the other triangle
    @inbounds for j=1:size(S,1)
        S[j,j]=S[j,j]/2.0;      # the diagonal is there twice; fix it;
    end
    self=SysmatAssemblerSparse(0.0*self.matbuffer[1])# get rid of the buffers
    S = S[1:end-1,1:end-1]
    return S
end

"""
    AbstractSysvecAssembler

Abstract type of system vector assembler.
"""
abstract type AbstractSysvecAssembler end;

"""
    SysvecAssembler

Assembler for the system vector.
"""
mutable struct SysvecAssembler{T<:Number} <: AbstractSysvecAssembler
    F_buffer::Vector{T};
    ndofs::FInt
end

"""
    SysvecAssembler(zero::T=0.0) where {T<:Number}

Construct blank system vector assembler. The vector entries are of type `T`.
"""
function SysvecAssembler(zero::T=0.0) where {T<:Number}
    return SysvecAssembler([zero], 1)
end

"""
    startassembly!(self::SysvecAssembler{T},
      ndofs_row::FInt) where {T<:Number}

Start assembly.

The method makes the buffer for the vector assembly. It must be called before
the first call to the method assemble.

`ndofs_row`= Total number of degrees of freedom.
"""
function startassembly!(self::SysvecAssembler{T},  ndofs_row::FInt) where {T<:Number}
    self.ndofs = ndofs_row
    self.F_buffer = zeros(T,self.ndofs);
end

"""
    assemble!(self::SysvecAssembler{T}, vec::MV,
      dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}

Assemble an elementwise vector.

The method assembles a column element vector using the vector of degree of
freedom numbers for the rows.
"""
function assemble!(self::SysvecAssembler{T}, vec::MV,
  dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}
    for i = 1:length(dofnums)
        gi = dofnums[i];
        if (0 < gi <= self.ndofs)
            self.F_buffer[gi] = self.F_buffer[gi] + vec[i];
        end
    end
end

"""
    makevector!(self::SysvecAssembler)

Make the global vector.
"""
function makevector!(self::SysvecAssembler)
  return deepcopy(self.F_buffer);
end

"""
    SysmatAssemblerSparseHRZLumpingSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a **symmetric lumped square** matrix  assembled from  **symmetric square**
matrices.

Reference: A note on mass lumping and related processes in the finite element method,
E. Hinton, T. Rock, O. C. Zienkiewicz, Earthquake Engineering & Structural Dynamics,
volume 4, number 3, 245--249, 1976.
}

"""
mutable struct SysmatAssemblerSparseHRZLumpingSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    buffer_length:: FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer:: FInt;
    ndofs:: FInt;
end

"""
    SysmatAssemblerSparseHRZLumpingSymm(zer::T=0.0) where {T<:Number}

Construct blank system matrix assembler. The matrix entries are of type `T`.
"""
function SysmatAssemblerSparseHRZLumpingSymm(zer::T=0.0) where {T<:Number}
    return SysmatAssemblerSparseHRZLumpingSymm{T}(0,[zer],[0],[0],0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparseHRZLumpingSymm{T},
      elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
      ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric lumped square global matrix.
"""
function startassembly!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt,ignore2::FInt) where {T<:Number}
    # Start the assembly of a global matrix.
    # The method makes buffers for matrix assembly. It must be called before
    # the first call to the method assemble.
    # elem_mat_nrows= number of rows in typical element matrix,
    # elem_mat_ncols= number of columns in a typical element matrix,
    # elem_mat_nmatrices= number of element matrices,
    # ndofs_row= Total number of equations in the row direction,
    # ndofs_col= Total number of equations in the column direction.
    self.buffer_length = elem_mat_nmatrices*elem_mat_dim^2;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},
      dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble a HRZ-lumped square symmetric matrix.
"""
function assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}
    # Assembly of a HRZ-lumped square symmetric matrix.
    # The method assembles the scaled diagonal of the square symmetric matrix using the two vectors of equation numbers for the rows and columns.
    nrows=length(dofnums); ncolumns=nrows;
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    # Now comes the lumping procedure
    em2 = sum(sum(mat, dims = 1));
    dem2 = sum(diag(mat));
    ffactor = em2/dem2 # total-element-mass compensation factor
    @inbounds for j=1:ncolumns
        @inbounds for i=j:nrows
            if i == j # only the diagonal elements are assembled
                self.matbuffer[p] = mat[i,j]*ffactor # serialized matrix
                self.rowbuffer[p] = dofnums[i];
                self.colbuffer[p] = dofnums[j];
                p=p+1
            end
        end
    end
    self.buffer_pointer=p;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},
        dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble an HRZ-lumped square symmetric matrix.
"""
function assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T}, dofnums::FIntMat, ignore::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums), vec(ignore))
end

"""
    makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)

Make a sparse HRZ-lumped **symmetric square**  matrix.
"""
function makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j]=self.ndofs+1;
        end
        if (self.colbuffer[j] > self.ndofs) || (self.colbuffer[j] <= 0)
            self.colbuffer[j]=self.ndofs+1;
        end
    end
    S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs+1, self.ndofs+1);
    S = S+copy(transpose(S)); # This is to address the lack of a function for adding together a sparse matrix with a transpose of the sparse matrix; at the moment (January 2018), this defaults to the addition of two dense matrices and the result is DENSE. Hence a copy of the transpose needs to be made.
    #   S = S+transpose(S);    # construct the other triangle
    @inbounds for j=1:size(S,1)
        S[j,j]=S[j,j]/2.0;      # the diagonal is there twice; fix it;
    end
    self=SysmatAssemblerSparse(0.0*self.matbuffer[1])# get rid of the buffers
    S = S[1:end-1,1:end-1]
    return S
end

end


text:
struct l1Manifold{n} <: Manifold{ℝ}
    nnz_coords::BitArray{1}
end
function l1Manifold(nnz_coords)
    return l1Manifold{length(nnz_coords)}(convert(BitArray, nnz_coords))
end

function get_nnz_indices(M::l1Manifold)
    inds = Vector{Int64}(undef, sum(M.nnz_coords))

    i = 1
    for (ind::Int64, val::Bool) in enumerate(M.nnz_coords)
        if val != false
            inds[i] = ind
            i+=1
        end
    end
    return inds
end

==(M::l1Manifold, N::l1Manifold) = (M.nnz_coords == N.nnz_coords)

function <(M::l1Manifold, N::l1Manifold)
    return M.nnz_coords < N.nnz_coords
end

function check_manifold_point(M::l1Manifold{n}, x; kwargs...) where {n}
    if size(x) != (n,)
        return DomainError(size(x), "x should be vector of size $n.")
    end
    norm_zerocoords = norm(x .* (1 .- M.nnz_coords))
    if norm_zerocoords > 0
        return DomainError(
            norm_zerocoords,
            "Coordinates of x supposed to be 0 have norm $norm_zerocoords.",
        )
    end
    return nothing
end

function check_tangent_vector(M::l1Manifold, x, ξ; check_base_point = true, kwargs...)
    if check_base_point
        mpe = check_manifold_point(M, x; kwargs...)
        mpe === nothing || return mpe
    end
    check_manifold_point(M, ξ; kwargs...)
    return nothing
end

copy(M::l1Manifold) = l1Manifold(copy(M.nnz_coords))

# decorated_manifold(M::l1Manifold) = Euclidean(representation_size(M)...; field = ℝ)

distance(::l1Manifold, p, q) = norm(p - q)

# Euclidean to riemannian gradients, hessians at vectors.
function ehess_to_rhess!(M::l1Manifold, Hessf_xξ, x, ∇f_x, ∇²f_ξ, ξ)
    return project!(M, Hessf_xξ, x, ∇²f_ξ)
end

# embed!(::l1Manifold, X, x) = copyto!(X, x)
# embed!(::l1Manifold, Xi, x, ξ) = copyto!(Xi, ξ)

embed(::l1Manifold, x) = x
embed(::l1Manifold, x, ξ) = ξ

embedding_dimension(M::l1Manifold) = length(M.nnz_coords)

exp!(::l1Manifold, y, x, ξ) = (@. y = x + ξ)

@inline inner(::l1Manifold, x, ξ, η) = dot(ξ, η)

@inline log!(M::l1Manifold, ξ, x, y) = (ξ .= y .- x)

manifold_dimension(M::l1Manifold{n}) where {n} = sum(M.nnz_coords)

function name(M::l1Manifold; short = true)
    return short ? "l1-$(sum(M.nnz_coords))/$(length(M.nnz_coords))" :
           "l1Manifold with $(sum(M.nnz_coords)) nnz"
end

norm(::l1Manifold, x, ξ) = norm(ξ)

function project!(M::l1Manifold, ξ, x, X)
    (@. ξ = X * M.nnz_coords)
    return ξ
end
project!(M::l1Manifold, x, X) = (@. x = X * M.nnz_coords)

function project(M::l1Manifold{n}, X) where {n}
    nnz_inds = get_nnz_indices(M)
    return sparsevec(nnz_inds, X[nnz_inds], n)
end

function randomMPoint(M::l1Manifold{n}) where {n}
    nnz = manifold_dimension(M)
    return sparsevec(get_nnz_indices(M), 2 * rand(nnz) .- 1, n)
end

function randomTVector(M::l1Manifold, x)
    pt = randomMPoint(M)
    pt /= norm(pt)
    return pt
end

representation_size(::l1Manifold{n}) where {n} = (n)

retract!(::l1Manifold, y, x, ξ) = (@. y = x + ξ)

function Base.show(io::IO, M::l1Manifold{n}) where {n}
    return print(io, name(M))
end

zero_tangent_vector(M::l1Manifold{n}, ::Any) where {n} = sparsevec(get_nnz_indices(M), zeros(manifold_dimension(M)), n)
zero_tangent_vector!(M::l1Manifold, v, p) = fill!(v, 0)

# const l1MPoint = Vector
# const l1TVector = Vector

# copy(m::l1Manifold) = l1Manifold(m.nnz_coords)

# ==(m1::l1Manifold, m2::l1Manifold) = m1.nnz_coords == m2.nnz_coords

# # Functions
# # ---
# manifold_dimension(M::l1Manifold) = sum(M.nnz_coords)

# # metric
# inner(M::l1Manifold, x::l1MPoint, ξ::l1TVector, ν::l1TVector) = dot(ξ, ν)
# norm(M::l1Manifold, x::l1MPoint, ξ::l1TVector) = inner(M, x, ξ, ξ)

# # Tangent space
# function project!(M::l1Manifold, ξ::l1TVector, x::l1MPoint, q::AbstractArray)
#     @. ξ = q * M.nnz_coords
#     return ξ
# end

# function project(M::l1Manifold, x::l1MPoint, q::AbstractArray)                  # ! This should factor out, with `zro_tangent_vector`
#     ξ = zeros(size(x))
#     return project!(M, ξ, x, q)
# end

# function retract!(M::l1Manifold, y::l1MPoint, x::l1MPoint, ξ::l1TVector)
#     return @. y = x + ξ
# end
# function retract(M::l1Manifold, x::l1MPoint, ξ::l1TVector)
#     y = zeros(size(x))
#     return retract!(M, y, x, ξ)
# end

# # Display
# # ---
# show(io::IO, M::l1Manifold) = print(io, M.abbreviation)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7952ee6-45b5-11eb-1158-5bb2ff274ce9
begin
	using DrWatson

		md"""
	*mtg_cards loaded from JSON*

	## EDH Deck: Vannifar's Circus
	"""
end

# ╔═╡ 0c145632-3927-11eb-19b9-877e05c1bcdc
begin
	@quickactivate

	using GameZero
	using JSON
	using HTTP
	using Images: load, RGB
	using PlutoUI
	using Serialization
	using PlaymatSimulator
	#using DataStructures
	using ImageTransformations: imresize

	AC = PlaymatSimulator.Actors

	pd = projectdir()

	mtg_cards = JSON.parsefile("$pd/Base/json/oracle-cards-20201224220555.json")

	md"**Found $(length(mtg_cards)) cards in JSON db!**"
end

# ╔═╡ 621b08a4-384e-11eb-0109-61e9b9ecf125
if !(@isdefined deck)
	deck = Dict{Symbol,Any}(
	:name => split(@__DIR__, "/")[end],
    :commander_names => [
        "Prime Speaker Vannifar",
    	],
    :card_names => [
		"Arixmethes, Slumbering Isle",
		"Alchemist's Refuge",
		"Birds of Paradise",
        "Botanical Sanctum",
        "Brainstorm",
        "Breeding Pool",
        "City of Brass",
        "Coiling Oracle",
        "Command Tower",
        "Counterspell",
        "Crop Rotation",
        "Cryptic Command",
        "Cultivate",
        "Deadeye Navigator",
        "Devoted Druid",
        "Dryad Arbor",
        "Dryad of the Ilysian Grove",
        "Elvish Mystic",
        "Elvish Reclaimer",
        "Eternal Witness",
        "Experiment Kraj",
        "Fae of Wishes // Granted",
        "Faerie Conclave",
        "Fblthp, the Lost",
        "Flooded Grove",
        "Flooded Strand",
        "Forbidden Orchard",
        "Forest",
		"Forest",
        "Glen Elendra Archmage",
        "Grand Architect",
		"Green Sun's Zenith",
		"Growth Spiral",
		"Gyre Engineer",
        "Hinterland Harbor",
		"Incubation Druid",
		"Island",
		"Island",
		"Jwari Disruption // Jwari Ruins",
        "Kinnan, Bonder Prodigy",
        "Kiora's Follower",
        "Kodama's Reach",
        "Leech Bonder",
        "Lesser Masticore",
        "Ley Weaver",
        "Lightning Greaves",
        "Llanowar Elves",
        "Llanowar Reborn",
        "Lore Weaver",
        "Magosi, the Waterveil",
        "Maze of Ith",
        "Meekstone",
        "Minamo, School at Water's Edge",
        "Misty Rainforest",
        "Murkfiend Liege",
        "Mystic Sanctuary",
        "Nykthos, Shrine to Nyx",
        "Paradise Mantle",
		"Parcelbeast",
		"Pemmin's Aura",
		"Phyrexian Metamorph",
		"Pili-Pala",
		"Prime Speaker Vannifar",
        "Ramunap Excavator",
        "Reclamation Sage",
        "Reflecting Pool",
        "Regrowth",
        "Remand",
        "Rishkar, Peema Renegade",
        "Safe Haven",
        "Sapseep Forest",
        "Seedborn Muse",
        "Sensei's Divining Top",
        "Simic Growth Chamber",
        "Skullclamp",
        "Snapcaster Mage",
        "Solemn Simulacrum",
        "Sol Ring",
        "Spellseeker",
        "Spore Frog",
        "Strip Mine",
        "Teferi, Mage of Zhalfir",
        "Temple of Mystery",
        "Temporal Mastery",
        "Thousand-Year Elixir",
        "Tolaria West",
        "Trinket Mage",
        "Tropical Island",
		"Vastwood Fortification // Vastwood Thicket",
        "Venser, Shaper Savant",
        "Scryb Ranger",
        "Vizier of Tumbling Sands",
		"Walking Ballista",
		"Waterlogged Grove",
        "Willbreaker",
        "Wirewood Symbiote",
        "Worldly Tutor",
        "Yavimaya Coast",
        "Yavimaya Elder",
        "Young Wolf",
    ]
)
else
	nothing
end

# ╔═╡ fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
length(deck[:card_names])

# ╔═╡ c61fa79e-4583-11eb-3b71-2d334aca843d
md"""
Look OK? Keep in mind that images that do not require in-game scaling will suffer less distortion.
"""

# ╔═╡ 7420cf10-45cc-11eb-2780-4f320bd8a2cf
begin
	cards = []

	for n in sort(deck[:card_names])

		for c in mtg_cards

			if n == c["name"]
				push!(cards, c)
			end
		end
	end

	md"""
	**Found $(length(cards)) matching cards in mtg_cards!**
	"""
end

# ╔═╡ c31dd202-50c6-11eb-0631-13c70535635e
missing_cards = filter!(x->!(x in [ c["name"] for c in cards ]), vcat(deck[:card_names], deck[:commander_names]))

# ╔═╡ c5374766-4ef1-11eb-2555-c159dba953f0
md"""Card #: $(@bind i Slider(1:length(cards), default=50, show_value=true))
Shrink card by $(@bind card_ratio Slider(0.1:0.05:1.25, default=0.5, show_value=true))

or: $(@bind customw CheckBox(default=true)) specify the default card width: $(@bind card_width NumberField(1:10; default=270)) in pixels
"""

# ╔═╡ 0899b1d0-5972-11eb-0470-4d480cf95d53
md"""
##### Download card imgs? $(@bind download_imgs CheckBox())
"""

# ╔═╡ be2af776-5971-11eb-13ec-3d7982a01ea3
md"""
##### Save data to disk? $(@bind save_data CheckBox())
"""

# ╔═╡ 43cdcbdc-5f68-11eb-11bd-51fbbf610333
deck

# ╔═╡ ce216c54-468a-11eb-13b8-7f3dac7af44a
function get_face_img(img_uri::String)
	img_resp = HTTP.get(img_uri)
	card_img = img_resp.body |> IOBuffer |> load
end

# ╔═╡ 2ab53d00-50cd-11eb-1cd4-5bf94ce53692
function get_card_preview_img(c)
	if haskey(all_cards[i], "card_faces") && haskey(all_cards[i]["card_faces"][1], "image_uris")
		hcat([get_face_img(f["image_uris"]["border_crop"])
			for f in all_cards[i]["card_faces"] ]...)
	else
		get_face_img(all_cards[i]["image_uris"]["border_crop"])
	end
end

# ╔═╡ 97bc3768-50ce-11eb-3f74-95d4fefe3792
function get_card_faces(c, faces=[])
	if haskey(c, "image_uris")
		push!(faces, get_face_img(c["image_uris"]["border_crop"]))

	elseif haskey(c, "card_faces")
		for f in c["card_faces"]
			push!(faces, get_face_img(f["image_uris"]["border_crop"]))
		end
	end

	return faces
end

# ╔═╡ 73d1cd18-4647-11eb-3994-7d4eb92eddca
if length(cards) > 0 && (@isdefined cards)

	if customw
		fs = get_card_faces(cards[i])
		sz = size(fs[end])
		imresize.(fs, ceil(Int, card_width * sz[1]/sz[2]), card_width)[end]
	else
		imresize.(get_card_faces(cards[i]), ratio=card_ratio)[end]
	end
else
	nothing
end

# ╔═╡ dfc9b56e-50ce-11eb-0e7f-83ec6e831901
if download_imgs
	CARD_FACES = []

	for c in cards
		faces = get_card_faces(c)

		if customw
			sz = size(faces[end])
			faces = imresize.(faces, ceil(Int, card_width * sz[1]/sz[2]), card_width)
		else
			faces = imresize.(faces, ratio=card_ratio)
		end

		push!(CARD_FACES, c["name"] => faces)
		sleep(0.1)
	end
end

# ╔═╡ 3a17fc4c-5f69-11eb-29a1-c9bb584fee32
if @isdefined CARD_FACES
	CARD_FACES[i][end][end]
end

# ╔═╡ d654ad1e-468a-11eb-2348-695621b7b9b0
function search_mtg_cards_by_keyword(q::String, mtg_cards::Array)
	[ n for n in [ c["name"] for c in mtg_cards ] if occursin(q, n) ]
end

# ╔═╡ 7cbdfda6-5eb5-11eb-2be6-0f93c3374a47
function hbox(x, y, gap=16; sy=size(y), sx=size(x))
	w,h = (max(sx[1], sy[1]),
		   gap + sx[2] + sy[2])

	slate = fill(RGB(1,1,1), w,h)
	slate[1:size(x,1), 1:size(x,2)] .= RGB.(x)
	slate[1:size(y,1), size(x,2) + gap .+ (1:size(y,2))] .= RGB.(y)
	slate
end

# ╔═╡ 5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
begin
	CARD_BACK_PATH = "$pd/Base/ui/cards/card_back.png"
	CARD_BACK_IMG = imresize(load(CARD_BACK_PATH), ratio=card_ratio)

	hbox(CARD_BACK_IMG, CARD_BACK_IMG)
end

# ╔═╡ 5e6f7046-4da5-11eb-0122-bd82397aab4f
if save_data
	deck[:CARD_BACK_IMG] = CARD_BACK_IMG
	deck[:CARD_FACES] = CARD_FACES

	fn = "$(projectdir())/games/MtG/EDH/decks/$(deck[:name])/$(deck[:name]).jls"
	serialize(fn, deck)

	md"Deck data saved to $fn"
end

# ╔═╡ 85103516-5eb5-11eb-3abc-dfc9ae9caf96
vbox(x,y, gap=16) = hbox(x', y')'

# ╔═╡ Cell order:
# ╟─c7952ee6-45b5-11eb-1158-5bb2ff274ce9
# ╟─0c145632-3927-11eb-19b9-877e05c1bcdc
# ╟─621b08a4-384e-11eb-0109-61e9b9ecf125
# ╟─fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
# ╟─5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
# ╟─c61fa79e-4583-11eb-3b71-2d334aca843d
# ╟─7420cf10-45cc-11eb-2780-4f320bd8a2cf
# ╟─c31dd202-50c6-11eb-0631-13c70535635e
# ╟─c5374766-4ef1-11eb-2555-c159dba953f0
# ╟─73d1cd18-4647-11eb-3994-7d4eb92eddca
# ╟─0899b1d0-5972-11eb-0470-4d480cf95d53
# ╟─3a17fc4c-5f69-11eb-29a1-c9bb584fee32
# ╟─dfc9b56e-50ce-11eb-0e7f-83ec6e831901
# ╟─be2af776-5971-11eb-13ec-3d7982a01ea3
# ╟─5e6f7046-4da5-11eb-0122-bd82397aab4f
# ╟─43cdcbdc-5f68-11eb-11bd-51fbbf610333
# ╟─ce216c54-468a-11eb-13b8-7f3dac7af44a
# ╟─2ab53d00-50cd-11eb-1cd4-5bf94ce53692
# ╟─97bc3768-50ce-11eb-3f74-95d4fefe3792
# ╟─d654ad1e-468a-11eb-2348-695621b7b9b0
# ╟─7cbdfda6-5eb5-11eb-2be6-0f93c3374a47
# ╟─85103516-5eb5-11eb-3abc-dfc9ae9caf96


text:
import ..Topologies
using ..RecursiveApply

function dss_transform(arg, local_geometry)
    RecursiveApply.rmap(arg) do x
        dss_transform(x, local_geometry)
    end
end
dss_transform(arg::Number, local_geometry) = arg
dss_transform(
    arg::Geometry.AxisTensor{T, N, <:Tuple{Vararg{Geometry.CartesianAxis}}},
    local_geometry::Geometry.LocalGeometry,
) where {T, N} = arg
dss_transform(
    arg::Geometry.CartesianVector,
    local_geometry::Geometry.LocalGeometry,
) where {T, N} = arg
dss_transform(
    arg::Geometry.AxisTensor{T, N, <:Tuple{Vararg{Geometry.LocalAxis}}},
    local_geometry::Geometry.LocalGeometry,
) where {T, N} = arg
dss_transform(
    arg::Geometry.LocalVector,
    local_geometry::Geometry.LocalGeometry,
) where {T, N} = arg

function dss_transform(
    arg::Geometry.AxisVector,
    local_geometry::Geometry.LocalGeometry,
)
    ax = axes(local_geometry.∂x∂ξ, 1)
    axfrom = axes(arg, 1)
    # TODO: make this consistent
    if ax isa Geometry.UVAxis && (
        axfrom isa Geometry.Covariant3Axis ||
        axfrom isa Geometry.Contravariant3Axis
    )
        return arg
    end
    if ax isa Geometry.UWAxis && (
        axfrom isa Geometry.Covariant2Axis ||
        axfrom isa Geometry.Contravariant2Axis
    )
        return arg
    end
    Geometry.transform(ax, arg, local_geometry)
end

function dss_untransform(refarg, targ, local_geometry)
    RecursiveApply.rmap(refarg, targ) do rx, tx
        dss_untransform(rx, tx, local_geometry)
    end
end
dss_untransform(refarg::Number, targ, local_geometry) = targ
dss_untransform(
    refarg::T,
    targ::T,
    local_geometry::Geometry.LocalGeometry,
) where {T <: Geometry.AxisTensor} = targ
function dss_untransform(
    refarg::Geometry.AxisTensor,
    targ::Geometry.AxisTensor,
    local_geometry::Geometry.LocalGeometry,
)
    ax = axes(refarg, 1)
    Geometry.transform(ax, targ, local_geometry)
end

function dss_1d!(
    dest,
    src,
    local_geometry_data,
    htopology::Topologies.AbstractTopology,
    Nq::Int,
    Nv::Int = 1,
)
    idx1 = CartesianIndex(1, 1, 1, 1, 1)
    idx2 = CartesianIndex(Nq, 1, 1, 1, 1)
    for (elem1, face1, elem2, face2, reversed) in
        Topologies.interior_faces(htopology)
        for level in 1:Nv
            @assert face1 == 1 && face2 == 2 && !reversed
            local_geometry_slab1 = slab(local_geometry_data, level, elem1)
            src_slab1 = slab(src, level, elem1)
            local_geometry_slab2 = slab(local_geometry_data, level, elem2)
            src_slab2 = slab(src, level, elem2)
            val =
                dss_transform(src_slab1[idx1], local_geometry_slab1[idx1]) ⊞
                dss_transform(src_slab2[idx2], local_geometry_slab2[idx2])

            dest_slab1 = slab(dest, level, elem1)
            dest_slab2 = slab(dest, level, elem2)

            dest_slab1[idx1] = dss_untransform(
                dest_slab1[idx1],
                val,
                local_geometry_slab1[idx1],
            )
            dest_slab2[idx2] = dss_untransform(
                dest_slab2[idx2],
                val,
                local_geometry_slab2[idx2],
            )
        end
    end
    return dest
end

"""
    horizontal_dss!(dest, src, topology, Nq)

Apply horizontal direct stiffness summation (DSS) to `src`, storing the result in `dest`.
"""
function dss_2d!(
    dest,
    src,
    local_geometry_data,
    topology::Topologies.AbstractTopology,
    Nq::Int,
    Nv::Int = 1,
)
    # TODO: generalize to extruded domains by returning a cartesian index?
    # iterate over the interior faces for each element of the mesh
    for (elem1, face1, elem2, face2, reversed) in
        Topologies.interior_faces(topology)
        for level in 1:Nv
            # iterate over non-vertex nodes
            src_slab1 = slab(src, level, elem1)
            src_slab2 = slab(src, level, elem2)
            local_geometry_slab1 = slab(local_geometry_data, level, elem1)
            local_geometry_slab2 = slab(local_geometry_data, level, elem2)
            dest_slab1 = slab(dest, level, elem1)
            dest_slab2 = slab(dest, level, elem2)
            for q in 2:(Nq - 1)
                i1, j1 = Topologies.face_node_index(face1, Nq, q, false)
                i2, j2 = Topologies.face_node_index(face2, Nq, q, reversed)
                val =
                    dss_transform(
                        src_slab1[i1, j1],
                        local_geometry_slab1[i1, j1],
                    ) ⊞ dss_transform(
                        src_slab2[i2, j2],
                        local_geometry_slab2[i2, j2],
                    )

                dest_slab1[i1, j1] = dss_untransform(
                    dest_slab1[i1, j1],
                    val,
                    local_geometry_slab1[i1, j1],
                )
                dest_slab2[i2, j2] = dss_untransform(
                    dest_slab2[i2, j2],
                    val,
                    local_geometry_slab2[i2, j2],
                )
            end
        end
    end

    # iterate over all vertices
    for vertex in Topologies.vertices(topology)
        # for each level
        for level in 1:Nv
            # gather: compute sum over shared vertices
            sum_data = mapreduce(⊞, vertex) do (elem, vertex_num)
                i, j = Topologies.vertex_node_index(vertex_num, Nq)
                local_geometry_slab = slab(local_geometry_data, level, elem)
                src_slab = slab(src, level, elem)
                dss_transform(src_slab[i, j], local_geometry_slab[i, j])
            end

            # scatter: assign sum to shared vertices
            for (elem, vertex_num) in vertex
                dest_slab = slab(dest, level, elem)
                i, j = Topologies.vertex_node_index(vertex_num, Nq)
                local_geometry_slab = slab(local_geometry_data, level, elem)
                dest_slab[i, j] = dss_untransform(
                    dest_slab[i, j],
                    sum_data,
                    local_geometry_slab[i, j],
                )
            end
        end
    end
    return dest
end

function horizontal_dss!(dest, src, space::AbstractSpace)
    if space isa ExtrudedFiniteDifferenceSpace
        Nv = nlevels(space)
        hspace = space.horizontal_space
    else
        Nv = 1
        hspace = space
    end
    htopology = hspace.topology
    Nq = Quadratures.degrees_of_freedom(hspace.quadrature_style)::Int
    if hspace isa SpectralElementSpace1D
        dss_1d!(dest, src, local_geometry_data(space), htopology, Nq, Nv)
    elseif hspace isa SpectralElementSpace2D
        dss_2d!(dest, src, local_geometry_data(space), htopology, Nq, Nv)
    end
end

horizontal_dss!(data, space::AbstractSpace) =
    horizontal_dss!(data, data, space::AbstractSpace)

weighted_dss!(dest, src, space::AbstractSpace) = horizontal_dss!(
    dest,
    Base.Broadcast.broadcasted(⊠, src, space.dss_weights),
    space,
)

weighted_dss!(dest, src, space::ExtrudedFiniteDifferenceSpace) =
    horizontal_dss!(
        dest,
        Base.Broadcast.broadcasted(⊠, src, space.horizontal_space.dss_weights),
        space,
    )

weighted_dss!(data, space::AbstractSpace) = weighted_dss!(data, data, space)


text:
import RowEchelon
import LinearAlgebra

#=
Homework 8
Section 4.2
Problems 5
=#

A = [1 -4 0 2 0 ;
     0 0 1 -5 0 ;
     0 0 0 0 2]

nulA = LinearAlgebra.nullspace(A)
println("Original Matrix, A:")
println(A)
println("NulA:")
println(nulA)


text:
configs = Dict()
setup = [
    (test_prob, behavior_activated)
    for test_prob in 0:0.1:1
    for behavior_activated in [true, false]
]
for (i, (test_prob, behavior_activated)) in enumerate(setup)
    configs[i] = EpidemicModel.Config(
        n_nodes=300,
        n0=3,
        k=3,
        agentcount=1000,
        transmission_prob=0.05,
        test_prob=test_prob,
        quarantine_duration=140,
        behavior_reduction_factor=0.25,
        behavior_activated=behavior_activated,
        tick_authority_recommendation=100,
        tick_authority_policy=200,
        n_public_spaces=100,
        initial_infectives=5,
        n_iter=1000
    )
end


text:
"Utilities"
module Util
using Spec

include("box.jl")
include("rng.jl")
include("namedtuple.jl") # Named tuple utilities
include("misc.jl")
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
export get_conjunctive_normal_form

"""
    get_conjunctive_normal_form(formula)

Return `formula` in the conjunctive normal form, expressed as clauses.

`formula` must not contain quantifiers. The returned formula is logically
equivalent to `formula`.
"""
function get_conjunctive_normal_form(formula)
# Ideally, we would use a single @match clause to handle both
# negations and junctions. The negation case is tricky to implement
# in such a unified approach, however, so we generate the CNF in two
# steps.
    function collapse_negations(formula)
        @match formula begin
        f::Negation => begin
            @match f.formula begin
            ff::Negation => collapse_negations(ff.formula)
            ff::Conjunction || ff::Disjunction => begin
                opposite_type = typeof(ff) == Conjunction ? Disjunction : Conjunction
                opposite_type(
                    collapse_negations(Negation(ff.formula1)),
                    collapse_negations(Negation(ff.formula2)),
                    )
                end
            ff::AtomicFormula => f
            end
        end
        f::Conjunction || f::Disjunction => typeof(f)(
            collapse_negations(f.formula1),
            collapse_negations(f.formula2),
        )
        f::AtomicFormula => f
        _ => throw(TypeError(:collapse_negations, "", Formula, typeof(formula)))
        end
    end
    function expand_disjunctions(formula)
        obj = @match formula begin
        f::Disjunction => begin
            cnf1 = expand_disjunctions(f.formula1)
            cnf2 = expand_disjunctions(f.formula2)
            # Concatenate the clauses in each pair of the cartesian product and
            # return the whole thing as a vector of clauses (that is, in CNF).
            resultcnf = CNF()
            for clause1 in cnf1
                for clause2 in cnf2
                    newclause = union(clause1, clause2)
                    push!(resultcnf, newclause)
                end
            end
            resultcnf
            end
        f::Conjunction => begin
            clauses = union(expand_disjunctions(f.formula1),
                            expand_disjunctions(f.formula2))
            CNF(clauses)
        end
        f::AtomicFormula => CNF([Clause([Literal(false, f)])])
        f::Negation => CNF([Clause([Literal(true, f.formula)])])
        end
        return obj
    end
    collapsed_formula = collapse_negations(formula)
    expand_disjunctions(collapsed_formula)
end


text:
module SimpleRoots

export Bisection, Secant, Brent

export Lower, Upper, Both

export findzero, quad

abstract type AbstractRootFinder end

abstract type AbstractBisection <: AbstractRootFinder end

bracket(method::AbstractBisection) = method.bracket

struct Bisection{B} <: AbstractBisection
    bracket::B
end
Bisection(arg1, arg2) = Bisection((arg1, arg2))

"""
    Secant(bracket)
    Secant(bracket...)

Secant method of bracketed root-finding.

Bracket can be a length 2 vector, a tuple or 2 argument
"""
struct Secant{B} <: AbstractBisection
    bracket::B
end
Secant(arg1, arg2) = Secant((arg1, arg2))

struct Brent{B} <: AbstractBisection
    bracket::B
end
Brent(arg1, arg2) = Brent((arg1, arg2))

"""
    findzero(f, bracket; atol=1e-7, max_iter=100)
    findzero(f, method::AbstractBisection; atol=1e-7, max_iter=100)

Find roots using the passed in method. If a bracket is passed in without
specifying the method, use the [`Secant`](@ref) method.
"""
findzero(f, bracket; kwargs...) = findzero(f, Secant(bracket); kwargs...)
findzero(f, method::AbstractRootFinder; atol=1e-7, max_iter=100) =
    findzero(f, method, atol, max_iter)

"""
    findzero(f, method::Secant, atol, max_iter)
    findzero(f, method::Secant; atol=1e-7, max_iter=100)

Find root using the secant method.
"""
function findzero(f, method::Secant, atol, max_iter)
    local x
    x0, x1 = bracket(method)
    y0, y1 = f(x0), f(x1)

    for _ in 1:max_iter
        x = x1 - y1 * (x1 - x0) / (y1 - y0)
        if abs(x - x1) < atol
            return x, true
        end
        x0 = x1
        y0 = y1
        x1 = x
        y1 = f(x1)
    end
    return x, false
end

"""
    findzero(f, method::Bisection, atol, max_iter)
    findzero(f, method::Bisection; atol=1e-7, max_iter=100)

Find root using the bisection method.
"""
function findzero(f, method::Bisection, atol, max_iter)
    a, b = bracket(method)
    fa = f(a)
    local c
    for _ in 1:max_iter
        if (b - a) <= atol
            return b, true
        end

        c = (a + b)/2
        fc = f(c)
        if fa * fc > zero(fa * fc)
            a = c # Root is in the right half of [a, b].
            fa = fc
        elseif fc == zero(fc)
            return c, true
        else
            b = c # Root is in the left half of [a, b].
        end
    end

    return b, false
end

"""
    findzero(f, method::Brent, atol, maxiter)
    findzero(f, method::Brent; atol=1e-7, max_iter=100)

Find root using Brents method. Returns a tuple of the found value
and a Bool specifying wether the root was found within the tolerance, or not.
"""
function findzero(f, method::Brent, atol, maxiter)
    local e, d
    epsval = eps()
    a, b = bracket(method)
    fa, fb = f(a), f(b)
    c = b
    fc = fb
    for iter = 1:maxiter
        if (fb > zero(fb) && fc > zero(fc)) || (fb < zero(fb) && fc < zero(fc))
            c = a
            fc = fa
            d = b - a
            e = d
        end
        if abs(fc) < abs(fb)
            a = b
            b = c
            c = a
            fa = fb
            fb = fc
            fc = fa
        end
        tol1 = 2epsval * abs(b) + atol/2
        xm = (c - b)/2
        if (abs(xm) <= tol1 || fb == zero(fb))
            return b, true
        end
        if abs(e) >= tol1 && abs(fa) > abs(fb)
            s = fb / fa
            if (a == c)
                p = 2xm * s
                q = 1 - s
            else
                q = fa / fc
                r = fb / fc
                p = s * (2xm * q * (q - r) - (b - a) * (r - 1))
                q = (q - 1) * (r - 1) * (s - 1)
            end
            if p > zero(p)
                q = -q
            end
            p = abs(p)
            if (2p < min(3xm * q - abs(tol1 * q), abs(e * q)))
                e = d
                d = p / q
            else
                d = xm
                e = d
            end
        else
            d = xm
            e = d
        end
        a = b
        fa = fb
        if (abs(d) > tol1)
            b = b + d
        else
            b = b + tol1 * sign(xm)
        end
        fb = f(b)
    end

    return b, false
end

"""
Used to specifiy the result of quad, either
[`Lower`](@ref), [`Upper`](@ref) or [`Both`](@ref)
"""
abstract type QuadraticResult end

"""
Return the lower result value
"""
struct Lower <: QuadraticResult end

"""
Return the upper result value
"""
struct Upper <: QuadraticResult end

"""
Return both lower and upper result values as a tuple
"""
struct Both <: QuadraticResult end

struct QuadraticError <: Exception
   var::String
end

Base.showerror(io::IO, e::QuadraticError) = print(io, e.var, "!")

"""
    function quad(res::QuadraticResult, a, b, c)

Simple quadratic equation solver.

Returns the result specified with `res`
[`Lower`](@ref), [`Upper`](@ref) or [`Both`](@ref)
"""
function quad(res::QuadraticResult, a, b, c)
    x = (b^2 - 4 * a * c)
    x < zero(x) && throw(QuadraticError("Imaginary roots in quadratic"))
    if a == zero(a)
        if b == zero(b)
            c != zero(c) && throw(QuadraticError("Can't solve quadratic"))
            # return zero with correct units
            return out(res, zero(-c / b))
        else
            return out(res, -c / b)
        end
    else
        return side(res, x, a, b)
    end
end
quad(a, b, c) = quad(Both(), a, b, c)

@inline side(::Both, x, a, b) = side(Lower(), x, a, b), side(Upper(), x, a, b)
@inline side(::Upper, x, a, b) = (-b + sqrt(x)) / 2a
@inline side(::Lower, x, a, b) = (-b - sqrt(x)) / 2a

@inline out(::Both, x) = (x, x)
@inline out(::Upper, x) = x
@inline out(::Lower, x) = x

end # module


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function point_intersect(fname)
    Drawing(1000, 1000, fname)
    origin()
    setopacity(0.4)

    for i in 1:100
        randomhue()
        p1 = Point(rand(-450:450), rand(-450:450))
        p2 = Point(rand(-450:450), rand(-450:450))
        p3 = Point(rand(-450:450), rand(-450:450))
        p4 = Point(rand(-450:450), rand(-450:450))
        flag, intersection_point = intersectionlines(p1, p2, p3, p4, crossingonly=true)
        if flag
            circle(intersection_point, 5, :fill)
            move(p1)
            line(p2)
            strokepath()
            move(p3)
            line(p4)
            strokepath()
        end
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

point_intersect("line_intersection.pdf")


text:
using .DistributionsAD: TuringDirichlet, TuringWishart, TuringInverseWishart,
                        FillVectorOfUnivariate, FillMatrixOfUnivariate,
                        MatrixOfUnivariate, FillVectorOfMultivariate, VectorOfMultivariate,
                        TuringScalMvNormal, TuringDiagMvNormal, TuringDenseMvNormal
using Distributions: AbstractMvLogNormal

# Bijectors

bijector(::TuringDirichlet) = SimplexBijector()
bijector(::TuringWishart) = PDBijector()
bijector(::TuringInverseWishart) = PDBijector()
bijector(::TuringScalMvNormal) = Identity{1}()
bijector(::TuringDiagMvNormal) = Identity{1}()
bijector(::TuringDenseMvNormal) = Identity{1}()

bijector(d::FillVectorOfUnivariate{Continuous}) = up1(bijector(d.v.value))
bijector(d::FillMatrixOfUnivariate{Continuous}) = up1(up1(bijector(d.dists.value)))
bijector(d::MatrixOfUnivariate{Discrete}) = Identity{2}()
bijector(d::MatrixOfUnivariate{Continuous}) = TruncatedBijector{2}(_minmax(d.dists)...)
bijector(d::VectorOfMultivariate{Discrete}) = Identity{2}()
for T in (:VectorOfMultivariate, :FillVectorOfMultivariate)
    @eval begin
        bijector(d::$T{Continuous, <:MvNormal}) = Identity{2}()
        bijector(d::$T{Continuous, <:TuringScalMvNormal}) = Identity{2}()
        bijector(d::$T{Continuous, <:TuringDiagMvNormal}) = Identity{2}()
        bijector(d::$T{Continuous, <:TuringDenseMvNormal}) = Identity{2}()
        bijector(d::$T{Continuous, <:MvNormalCanon}) = Identity{2}()
        bijector(d::$T{Continuous, <:AbstractMvLogNormal}) = Log{2}()
        bijector(d::$T{Continuous, <:SimplexDistribution}) = SimplexBijector{2}()
        bijector(d::$T{Continuous, <:TuringDirichlet}) = SimplexBijector{2}()
    end
end
bijector(d::FillVectorOfMultivariate{Continuous}) = up1(bijector(d.dists.value))

isdirichlet(::VectorOfMultivariate{Continuous, <:Dirichlet}) = true
isdirichlet(::VectorOfMultivariate{Continuous, <:TuringDirichlet}) = true
isdirichlet(::TuringDirichlet) = true

function link(
    d::TuringDirichlet,
    x::AbstractVecOrMat{<:Real},
    proj::Bool = true,
)
    return SimplexBijector{proj}()(x)
end

function link_jacobian(
    d::TuringDirichlet,
    x::AbstractVector{T},
    proj::Bool = true,
) where {T<:Real}
    return jacobian(SimplexBijector{proj}(), x)
end

function invlink(
    d::TuringDirichlet,
    y::AbstractVecOrMat{<:Real},
    proj::Bool = true
)
    return inv(SimplexBijector{proj}())(y)
end
function invlink_jacobian(
    d::TuringDirichlet,
    y::AbstractVector{T},
    proj::Bool = true
) where {T<:Real}
    return jacobian(inv(SimplexBijector{proj}()), y)
end

ispd(::TuringWishart) = true
ispd(::TuringInverseWishart) = true
function getlogp(d::TuringWishart, Xcf, X)
    return 0.5 * ((d.df - (dim(d) + 1)) * logdet(Xcf) - tr(d.chol \ X)) + d.logc0
end
function getlogp(d::TuringInverseWishart, Xcf, X)
    Ψ = d.S
    return -0.5 * ((d.df + dim(d) + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) + d.logc0
end