text:
#!/usr/bin/env julia

Fract_Zener = RheoModelClass(
          # Model name
          name="frac_zener",
          # Model parameters,
          p = [ :cₐ, :a, :cᵦ, :β, :cᵧ, :γ],
          # Relaxation modulus
          G = quote
                 cᵦ*t^(-β)*mittleff(a - β, 1 - β, -cᵦ*t^(a - β)/cₐ) + cᵧ*t^(-γ)/gamma(1 - γ)
              end,
          # Creep modulus
          J = quote
                  Ĵ(s) = (1/s)*(cₐ*s^a + cᵦ*s^β)/(cₐ*s^a*cᵦ*s^β + cᵧ*s^γ*(cₐ*s^a + cᵦ*s^β))
                  InverseLaplace.talbot(s -> Ĵ(s), t)
              end,
          # Storage modulus
          Gp = quote
                   denominator = (cₐ*ω^a)^2 + (cᵦ*ω^β)^2 + 2*(cₐ*ω^a)*(cᵦ*ω^β)*cos((a - β)*π/2)
                   numerator = ((cᵦ*ω^β)^2)*(cₐ*ω^a)*cos(a*π/2) + ((cₐ*ω^a)^2)*(cᵦ*ω^β)*cos(β*π/2)
                   numerator/denominator + cᵧ*ω^γ*cos(γ*π/2)
               end,
         # Loss modulus
          Gpp = quote
                   denominator = (cₐ*ω^a)^2 + (cᵦ*ω^β)^2 + 2*(cₐ*ω^a)*(cᵦ*ω^β)*cos((a - β)*π/2)
                   numerator = ((cᵦ*ω^β)^2)*(cₐ*ω^a)*sin(a*π/2) + ((cₐ*ω^a)^2)*(cᵦ*ω^β)*sin(β*π/2)
                   numerator/denominator + cᵧ*ω^γ*sin(γ*π/2)
                end,
          # Constraints
          constraint = quote
                   all([   (a<1) & (a>0)
                           (β<1) & (β>0)
                            -a+β < 0
                           (γ<1) & (γ>0)] )
                  end,
          # Network
          info= "

                  ______╱╲__________╱╲______
                 |      ╲╱          ╲╱      |
          _______|      cₐ,a         cᵦ, β  |_______
                 |                          |
                 |____________╱╲____________|
                              ╲╱
                              cᵧ, γ
                     "
          )

FractSLS_Zener = RheoModelClass(
        # Model name
        name="fracsls_Zener",
        # Model parameters,
        p = [ :cₐ, :a, :kᵦ, :kᵧ],
        # Relaxation modulus
        G = quote
                kᵦ*mittleff(a, -(kᵦ/cₐ)*t^a) + kᵧ
            end,
        # Creep modulus
        J = quote
                Ĵ(s) = (1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ))
                 InverseLaplace.talbot(s -> Ĵ(s), t)
            end,
        # Storage modulus
        Gp = quote
                denominator = (cₐ*ω^a)^2 + kᵦ^2 + 2*(cₐ*ω^a)*kᵦ*cos(a*π/2)
                numerator = kᵦ^2*(cₐ*ω^a)*cos(a*π/2) + ((cₐ*ω^a)^2)*kᵦ
                numerator/denominator + kᵧ
             end,
       # Loss modulus
        Gpp = quote
                denominator = (cₐ*ω^a)^2 + kᵦ^2 + 2*(cₐ*ω^a)*kᵦ*cos(a*π/2)
                numerator = kᵦ^2*(cₐ*ω^a)*sin(a*π/2)
                numerator/denominator
              end,
        # Constraints
        constraint = quote
                 (a<1) & (a>0)
                end,
        # Network
        info= "

                _____╱╲_______╱╲  ╱╲  ╱╲  _____
               |     ╲╱         ╲╱  ╲╱  ╲╱     |
        _______|       cₐ,a              kᵦ    |_______
               |                               |
               |__________╱╲  ╱╲  ╱╲  _________|
                            ╲╱  ╲╱  ╲╱
                                 kᵧ
                   "
        )

SLS_Zener = RheoModelClass(
          # Model name
          name="SLS_Zener",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = quote
                 kᵧ + kᵦ*exp(-t*kᵦ/η)
              end,
          # Creep modulus
          J = quote
                  c₀ = 1/kᵧ
                  c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
                  τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
                  c₀ - c₁*exp(-t/τᵣ)
              end,
          # Storage modulus
          Gp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω^2*τ^2*kᵦ
                  numerator/denominator + kᵧ
               end,
         # Loss modulus
          Gpp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω*τ*kᵦ
                  numerator/denominator
                end,
          # Network
          info= "
                      ___
                  _____| |________╱╲  ╱╲  ╱╲  ___
                 |    _|_|          ╲╱  ╲╱  ╲╱   |
          _______|      η                  kᵦ    |_______
                 |                               |
                 |__________╱╲  ╱╲  ╱╲  _________|
                              ╲╱  ╲╱  ╲╱
                                   kᵧ
                     "
          )

FractJeffreys_Zener = RheoModelClass(
                    # Model name
                    name="fjeff_Zener",
                    # Model parameters,
                    p = [ :ηₐ, :cᵦ, :β, :ηᵧ],
                    # Relaxation modulus
                    G = quote
                            diracterm = t!=0.0 ? 0.0 : Inf
                            cᵦ*t^(-β)*mittleff(1-β, 1-β, -cᵦ*t^(1-β)/ηₐ) + ηᵧ*diracterm
                        end,
                    # Creep modulus
                    J = quote
                            Ĵ(s) = (1/s)*(ηₐ*s + cᵦ*s^β)/(ηₐ*s*cᵦ*s^β + ηᵧ*s*(ηₐ*s + cᵦ*s^β))
                            InverseLaplace.talbot(s -> Ĵ(s), t)
                        end,
                    # Storage modulus
                    Gp = quote
                            denominator = (ηₐ*ω)^2 + (cᵦ*ω^β)^2 + 2*(ηₐ*ω)*(cᵦ*ω^β)*cos((1 - β)*π/2)
                            numerator = ((ηₐ*ω)^2)*(cᵦ*ω^β)*cos(β*π/2)
                            # cosine floating point error work-around
                            if β==1.0
                                0.0
                            else
                                numerator/denominator
                            end
                         end,
                   # Loss modulus
                    Gpp = quote
                            denominator = (ηₐ*ω)^2 + (cᵦ*ω^β)^2 + 2*(ηₐ*ω)*(cᵦ*ω^β)*cos((1 - β)*π/2)
                            numerator = ((cᵦ*ω^β)^2)*(ηₐ*ω) + ((ηₐ*ω)^2)*(cᵦ*ω^β)*sin(β*π/2)
                            numerator/denominator + ηᵧ*ω
                          end,
                    # Constraints
                    constraint = quote
                             (β<1) & (β>0)
                            end,
                    # Network
                    info= "

                                    ___
                            _________| |_________╱╲________
                           |        _|_|         ╲╱        |
                    _______|          ηₐ            cᵦ, β  |_______
                           |              ___              |
                           |_______________| |_____________|
                                          _|_|
                                             ηᵧ
                               "
                    )

Jeffreys_Zener = RheoModelClass(
                # Model name
                name="jeffreys_Zener",
                # Model parameters,
                p = [:ηₐ, :k, :ηᵧ],
                # Relaxation modulus
                G = quote
                        diracterm = t!=0.0 ? 0.0 : Inf
                        k*exp(-k*t/ηₐ) + ηᵧ*diracterm
                    end,
                # Creep modulus
                J = quote
                        Ĵ(s) = (1/s)*(ηₐ*s + k)/(ηₐ*s*k + ηᵧ*s*(ηₐ*s + k))
                        InverseLaplace.talbot(s -> Ĵ(s), t)
                    end,
                # Storage modulus
                Gp = quote
                        denominator = (ηₐ*ω)^2 + k^2
                        numerator = ((ηₐ*ω)^2)*k
                        numerator/denominator
                     end,
               # Loss modulus
                Gpp = quote
                        denominator = (ηₐ*ω)^2 + k^2
                        numerator = (k^2)*(ηₐ*ω)
                        numerator/denominator + ηᵧ*ω
                      end,
                # Network
                info= "

                              ___
                        _______| |_______╱╲  ╱╲  ╱╲  ___
                       |      _|_|         ╲╱  ╲╱  ╲╱  |
                _______|          ηₐ            k      |_______
                       |              ___              |
                       |_______________| |_____________|
                                      _|_|
                                         ηᵧ
                           "
                )

FractSolid = RheoModelClass(
        # Model name
        name="fractsolid",
        # Model parameters,
        p = [:η, :cᵦ, :β, :k],
        # Relaxation modulus
        G = quote
              k + cᵦ*t^(-β)*mittleff(1 - β, 1 - β, -cᵦ*(t^(1 - β))/η)
            end,
        # Creep modulus
        J = quote
              a = η/cᵦ
              b = k*η/cᵦ
              Ĵ(s) = (1/s^2)*(1+a*s^(1-β))/(η+k/s+b*s^(-β))
              InverseLaplace.talbot(s -> Ĵ(s), t)
            end,
        # Storage modulus
        Gp = quote
               denominator = (η*ω)^2 + (cᵦ*ω^β)^2
               numerator = ((η*ω)^2)*(cᵦ*ω^β)*cos(β*π/2)
               numerator/denominator + k
             end,
       # Loss modulus
        Gpp = quote
                denominator = (η*ω)^2 + (cᵦ*ω^β)^2
                numerator = ((cᵦ*ω^β)^2)*(η*ω) + ((η*ω)^2)*(cᵦ*ω^β)*sin(β*π/2)
                numerator/denominator
              end,
        # Constraints
        constraint = quote
                 (β<1) & (β>0)
                end,
        # Network
        info= "
                      ___
                  _____| |__________╱╲__________
                 |    _|_|          ╲╱          |
             ___ |      η              cᵦ, β    |___
                 |                              |
                 |__________╱╲  ╱╲  ╱╲  ________|
                              ╲╱  ╲╱  ╲╱
                                k
               "
        )


text:
using Revise
using DrWatson
@quickactivate "SICMModelingToolkit"
DrWatson.greet()
using ModelingToolkit, Unitful, DifferentialEquations, Optim
using Latexify
using BenchmarkTools, LinearAlgebra, ForwardDiff, Zygote
using Quadrature, StaticArrays
using SICMModelingToolkit
using Plots
using Measures

#### section 1.4

@parameters t m k
D = Differential(t)

locals1 = generate_generic_localtuple(1)
locals2 = generate_generic_localtuple(2)
locals3 = generate_generic_localtuple(3)

velocity(locals3)

L_free_mass = L_free_particle(m)

Lagrange_equations(L_free_mass, locals3)

L_harmonic_spring_apply = L_harmonic_spring(m, k)
Lagrange_equations(L_harmonic_spring_apply, locals1)
Lagrange_equations(L_harmonic_spring_apply, locals3)

ode_generator = ode_problem_generator(L_free_particle, locals3, [m])
#Lagrangian_pre_param, locals, params = L_free_particle, locals3, [m]

#ode_func = ODEFunction(ode_generator[1])

#du = [0., 0., 0., 0., 0., 0.]
#u0 = @SVector [1., 0.5, 0.25, 0., 0., 0.]
#p0 = @SVector [1.0]
#ode_func(du, u0, p0, 0.0)

#@show du

#dyn_ode_func(du, u, p, t) = DynamicalODEFunction{true}(ode_func)

du = [0., 0., 0.]
q0 = [0., 0., 0.]
qv0 = [1.0, 0.5, 0.25]
#u0 = [1., 0.5, 0.25, 0., 0., 0.]
u0 = ArrayPartition(qv0, q0)
p0 = [1.0]
t0 = 0.0

tspan = (0.0, 2.0)

prob = ode_generator(qv0, q0, tspan, p0; saveat=0.05)

sol = solve(prob, McAte5(); dt=0.05)

plot(sol, vars=(4,5,6); dpi=350)

### kepler

kepler_ode_generator = ode_problem_generator(L_kepler, locals2, [m, k])

begin
    q0_kep = [0.0, 1.0]
    qv0_kep = [-1.0, 0.0]
	tspan_kep = (0.0, 5.0)
    p_kep = [3.0, 10.0]
    dt = 0.01

    kepler_prob = kepler_ode_generator(qv0_kep, q0_kep, tspan_kep, p_kep)
end

kepler_sol = solve(kepler_prob, McAte5(); dt=dt)
kepler_sol_nosymp = solve(kepler_prob, Tsit5(); dt=dt)

curxlims = (-1.1, 1.1)
curylims = (9.0/16) .* curxlims
cur_left_margin = Measures.Length{:mm, Float64}(12)
cur_bottom_margin = Measures.Length{:mm, Float64}(7)
gr()
ellipseplot = plot(kepler_sol, vars=(4,3); dpi=500, ylims=curylims, xlims=curxlims,
    label="orbiter", xlabel="x", ylabel="y", title="Kepler Orbital Problem Symplectic",
    left_margin=cur_left_margin, bottom_margin=cur_bottom_margin)
scatter!(ellipseplot, [0.0], [0.0]; label="planet")
save("kepler_symplectic.png", ellipseplot)
ellipseplot_nosymp = plot(kepler_sol_nosymp, vars=(4,3); dpi=500, ylims=curylims, xlims=curxlims,
    label="orbiter", xlabel="x", ylabel="y", title="Kepler Orbital Problem Nonsymplectic",
    left_margin=cur_left_margin, bottom_margin=cur_bottom_margin)
scatter!(ellipseplot_nosymp, [0.0], [0.0]; label="planet")
save("kepler_nosymplectic.png", ellipseplot_nosymp)
plot(kepler_sol, vars=(0,3,4); dpi=350)

#### harmonic

harmonic_ode_generator = ode_problem_generator(L_harmonic_spring, locals1, [m, k])

begin
    q0_harm = [1.0]
    qv0_harm = [0.0]
	tspan_harm = (0.0, 5.0)
    p_harm = [1.0, 1.0]
    dt = 0.01

    harmonic_prob = harmonic_ode_generator(qv0_harm, q0_harm, tspan_harm, p_harm)
end

harmonic_sol = solve(harmonic_prob, McAte5(); dt=dt)
plot(harmonic_sol; dpi=500, left_margin=cur_left_margin, bottom_margin=cur_bottom_margin)

ext: jl

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
include("required.jl")

## Conditional arguments
# outcome = "rf"
# outcome = "bart"
# save_figures = true

# Make sure "outcome" is defined
if !(@isdefined outcome)
    @warn "'outcome' not defined, must be either 'raw', 'bio', 'rf', or 'bart'"
elseif !(outcome in ["raw", "bio", "rf", "bart"])
    @warn "'outcome' invalid, must be either 'raw', 'bio', 'rf', or 'bart'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Prepare data
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name
distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "proc", "distributions_$(outcome).tif"), i
    ) for i in eachindex(spenames)
]

Y = Ymatrix(distributions)
richness = calculate_richness(Y, distributions[1])
lcbd = calculate_lcbd(Y, distributions[1]; transform=true)

plot(richness; c=:viridis)

## Prepare subareas

# Define subareas
coords_NE = (left=-80.0, right=-60.0, bottom=40.0, top=50.0)
coords_SW = (left=-120.0, right=-100.0, bottom=30.0, top=40.0)
distributions_NE = [d[coords_NE] for d in distributions]
# Define window size
wsize = size(distributions_NE[1])

## Moving windows

function get_windows_indices(index_mat, wsize, steps::Tuple{Int64,Int64})
    lat_step, lon_step = steps
    # Get all windows
    winds = []
    jrange = 1:lon_step:(size(index_mat, 2) - (wsize[2] - 1))
    irange = 1:lat_step:(size(index_mat, 1) - (wsize[1] - 1))
    for j in jrange, i in irange
        subrows, subcols = i:(i + (wsize[1] - 1)), j:(j + (wsize[2] - 1))
        subinds = index_mat[subrows, subcols]
        push!(winds, subinds)
    end
    return winds
end
function get_windows_indices(index_mat, wsize, steps::Int64)
    return get_windows_indices(index_mat, wsize, (steps, steps))
end

# Apply moving windows
@time begin
    # Get matrix Y
    Y = Ymatrix(distributions)
    # Create matrix of indices
    index_mat = reshape(eachindex(distributions[1].grid), size(distributions[1])) |> Array
    # Get windows indices
    steps = Int64.(round.(1.0 ./ stride(distributions[1])))
    windows_inds = get_windows_indices(index_mat, wsize, steps)
    # Extract windows from Y
    Ywindows = [Y[vec(winds), :] for winds in windows_inds]
    # Remove ones without observation
    allnothing = map(x -> all(isnothing, x), Ywindows)
    deleteat!(windows_inds, allnothing)
    deleteat!(Ywindows, allnothing)
    # Calculate LCBD values for Ywindows
    LCBDwindows = @showprogress [
        calculate_lcbd(Yobs, distributions_NE[1]; transform=true, relative=true) for
        Yobs in Ywindows
    ]

    # Expand LCBD windows to match full extent
    LCBDbatch = []
    for i in eachindex(LCBDwindows)
        mat = fill(nothing, size(distributions[1])) |> Array{Union{Nothing,Float32}}
        mat[windows_inds[i]] = LCBDwindows[i].grid
        push!(LCBDbatch, mat)
    end
    # Stack windows in single matrix
    LCBDmat = reduce(hcat, map(vec, LCBDbatch))
    # Get mean LCBD value per site
    LCBDmean = map(x -> filter(!isnothing, x), eachrow(LCBDmat))
    LCBDmean = map(x -> isempty(x) ? nothing : mean(x), LCBDmean)
    # Arrange mean values as layer
    LCBDwindow = similar(distributions[1])
    LCBDwindow.grid = reshape(LCBDmean, size(LCBDwindow.grid))
end;
plot(LCBDwindow; c=:viridis)

## Visualize result
if outcome == "raw"
    plotfct = :plot_layer
else
    plotfct = :plot
end
function plot_lcbd_relationship(richness, lcbd; maintitle="", kw...)
    p1 = eval(plotfct)(
        lcbd; c=:viridis, title="LCBD", colorbar_title="Relative LCBD score", clim=(0, 1)
    )
    p2 = histogram2d(
        richness,
        lcbd;
        c=:viridis,
        bins=40,
        title="Relationship",
        xlabel="Richness",
        ylabel="LCBD",
        colorbar_title="Number of sites",
        xlim=(1, 45),
        ylim=(0.0, 1.0),
        bottommargin=3.0mm,
    )
    if maintitle != ""
        l = @layout [t{0.01h}; grid(1, 2)]
        ptitle = plot(; annotation=(0.5, 0.5, "$maintitle"), framestyle=:none)
        p = plot(ptitle, p1, p2; layout=l, size=(900, 300), kw...)
    else
        l = @layout [a b]
        p = plot(p1, p2; layout=l, size=(900, 300), kw...)
    end
    return p
end
window_full = plot_lcbd_relationship(richness, LCBDwindow)
window_NE = plot_lcbd_relationship(
    richness[coords_NE], LCBDwindow[coords_NE]; maintitle="Northeast subarea"
)
window_SW = plot_lcbd_relationship(
    richness[coords_SW], LCBDwindow[coords_SW]; maintitle="Southwest subarea"
)
psubareas = plot(
    window_NE,
    window_SW;
    layout=grid(2, 1),
    size=(900, 600),
    bottommargin=0.0mm,
    title=["" "" "" ""],
)

# Export figures
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    savefig(
        plot(window_full; dpi=200),
        joinpath("archive", "fig", outcome, "06-0_$(outcome)_moving-windows_full.png"),
    )
    savefig(
        plot(psubareas; dpi=200),
        joinpath("archive", "fig", outcome, "06-1_$(outcome)_moving-windows_subareas.png"),
    )
end

## GIF
left = -71.0;
right = -64.0;
bottom = 46.5;
top = 50.0;
dim_ratio = (top - bottom) / (right - left)
asp_ratio = 92.60 / 60.75
coords_subarea = (left=left, right=right, bottom=bottom, top=top)
subarea_plots = []
nplots = 0
@time while left > -145.0 + asp_ratio && bottom > 20.0 + asp_ratio * dim_ratio
    global nplots += 1
    global left -= asp_ratio
    global bottom -= asp_ratio * dim_ratio
    local coords_subarea = (left=left, right=right, bottom=bottom, top=top)
    local p = plot_lcbd_relationship(
        richness[coords_subarea],
        LCBDwindow[coords_subarea];
        formatter=f -> "$(round(f, digits=1))",
        dpi=200,
    )
    push!(subarea_plots, p)
end
anim = @animate for p in subarea_plots[Not(1)]
    plot(p)
end
gif(anim; fps=3)
if (@isdefined save_figures) && save_figures == true
    gif(anim, joinpath("fig", outcome, "06-3_$(outcome)_moving-windows.gif"); fps=3)
end

#### 3 scales comparison

# Extract LCBD & relationship subplots for first, middle, last GIF plots
mid_ind = median(1:length(subarea_plots)) |> round |> Int64
ps = subarea_plots[[1, mid_ind, end]]

# Combine 3 scales
p = plot(
    ps...;
    dpi=200,
    layout=(3, 1),
    size=(900, 900),
    title=["LCBD" "Relationship" "" "" "" ""],
)

# Export figures
# save_figures = true
if (@isdefined save_figures) && save_figures == true
    savefig(p, joinpath("archive", "fig", outcome, "06-2_$(outcome)_moving-windows_3scales.png"))
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

__precompile__(true)
module Gaston

export closefigure, closeall, figure,
       plot, plot!, scatter, scatter!, stem, bar, histogram, imagesc,
       surf, surf!, contour, scatter3, scatter3!,
       printfigure, set

import Base.show, Base.isempty

using Random
using DelimitedFiles

const VERSION = v"0.11"

## Handle Unix/Windows differences
#
# Define gnuplot's end-of-plot delimiter. It is different in Windows
# than in Unix, thanks to different end-of-line conventions.
gmarker_start = "GastonBegin\n"
gmarker_done = "GastonDone\n"
if Sys.iswindows()
    gmarker_start = "GastonBegin\r\n"
    gmarker_done = "GastonDone\r\n"
end

# load files
include("gaston_types.jl")
include("gaston_aux.jl")
include("gaston_config.jl")
include("gaston_figures.jl")
include("gaston_llplot.jl")
include("gaston_2d.jl")
include("gaston_3d.jl")
include("gaston_histograms.jl")
include("gaston_images.jl")
include("gaston_print.jl")

# define function to determine if function is empty
Base.isempty(f::Figure) = (f.curves == nothing)

# initialize internal state
gnuplot_state = GnuplotState()

mutable struct Pipes
    gstdin :: Pipe
    gstdout :: Pipe
    gstderr :: Pipe
    Pipes() = new()
end

const P = Pipes()

# initialize gnuplot
function __init__()
    global P
    try
        success(`gnuplot --version`)
    catch
        error("Gaston cannot be loaded: gnuplot is not available on this system.")
    end
    gstdin = Pipe()
    gstdout = Pipe()
    gstderr = Pipe()
    gproc = run(pipeline(`gnuplot`,
                         stdin = gstdin, stdout = gstdout, stderr = gstderr),
                wait = false)
    process_running(gproc) || error("There was a problem starting up gnuplot.")
    close(gstdout.in)
    close(gstderr.in)
    close(gstdin.out)
    P.gstdin = gstdin
    P.gstdout = gstdout
    P.gstderr = gstderr

    global IsJupyterOrJuno = false
    if isdefined(Main, :IJulia) && Main.IJulia.inited
        global IsJupyterOrJuno = true
    elseif isdefined(Main, :Juno) && Main.Juno.isactive()
        global IsJupyterOrJuno = true
    end

    global config = default_config()

    return nothing
end

end


text: # include("convolution.jl")

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork


text:
include("test_mean.jl")
include("test_variance.jl")
include("test_check.jl")


text:
## Get user commands
# Populates `parsed_args`
include(joinpath(dirname(@__DIR__),"parsecommands.jl"))

# Assume `fec` will be defined by script generating the data, else load --previous
if length(ARGS) == 0 && !@isdefined(fec)
  @info "loading" parsed_args["previous"]
  parsed_args["folder_name"] = "labrun2"
  parsed_args["remoteprocs"] = 0
  parsed_args["localprocs"] = 1
  parsed_args["batch_resolve"] = false
  parsed_args["vis2d"] = true
  parsed_args["vis3d"] = false
  parsed_args["imshow"] = true
  parsed_args["msgloops"] = 3000
  parsed_args["usesimmodels"] = true
end

# bring up the multiprocess cluster
!@isdefined(fec) && include(joinpath(@__DIR__, "CarFeedCommonSetup.jl"))

dfg = if @isdefined(fec)
  fec.slam.dfg
else
  dfg = initfg()
  getSolverParams(dfg).logpath=parsed_args["previous"]
  loadDFG!(dfg, joinLogPath(dfg, "finalGraph"))
  dfg
end

## =========================================================================================
# interpose results
## =========================================================================================

allD = jsonResultsSLAM2D(dfg)

allStr = JSON2.write(allD)

fid = open(joinLogPath(dfg, "$(runfile)_results_prebatch.json"),"w")
println(fid, allStr)
close(fid)

## save the factor graph

if parsed_args["savedfg"]
  saveDFG(dfg, joinLogPath(dfg, "fg_final_prebatch"))
end

## draw results

if parsed_args["vis2d"]

pl = plotSLAM2D(dfg, drawhist=false)
pl |> PDF(joinLogPath(dfg,"fg_final_prebatch.pdf"))

if parsed_args["report_factors"]
  reportFactors(dfg, Pose2Pose2, show=false)
  reportFactors(dfg, Pose2Pose2, lsf(dfg, r"l"),prefix="LM_", show=false)
end

end

## batch resolve

if parsed_args["batch_resolve"]

# fg2 = deepcopy(dfg)

enableSolveAllNotDRT!(dfg)
# dontMarginalizeVariablesAll!(dfg)
# foreach(x->setSolvable!(dfg, x, 1), ls(dfg))
# foreach(x->setSolvable!(dfg, x, 1), lsf(dfg))

tree, smt, hist = solveTree!(dfg, storeOld=true)

if parsed_args["savedfg"]
  saveDFG(dfg, joinLogPath(dfg, "fg_final_resolve"))
end

# report a second time after the resolve
if parsed_args["report_factors"]
  reportFactors(dfg, Pose2Pose2, show=false)
  reportFactors(dfg, Pose2Pose2, lsf(dfg, r"l"),prefix="LM_", show=false)
end

end

## after resolve interpose results

allD = jsonResultsSLAM2D(dfg)

allStr = JSON2.write(allD)

fid = open(joinLogPath(dfg, "$(runfile)_results.json"),"w")
println(fid, allStr)
close(fid)

## draw a second time

if parsed_args["batch_resolve"] && parsed_args["vis2d"]

pl = plotSLAM2D(dfg, drawhist=false)
pl |> PDF(joinLogPath(dfg,"fg_final_resolve.pdf"))

end

@show getLogPath(dfg)


text:
#ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется
# ровно одна внутренняя перегородка в форме прямоугольника.
#Робот - в произвольной клетке поля между внешней и
#внутренней перегородками.
#РЕЗУЛЬТАТ: Робот - в исходном положении
#и по всему периметру внутренней перегородки
#поставлены маркеры.
function try_move!(robot, side)::Bool
    if isborder(robot, side) # - это вызов виртуальной функции (кторая реально будет определена в каком-то из производных типов)
        false
    else
        move!(robot, side) # - это вызов виртуальной функции (кторая реально будет определена в каком-то из производных типов)
        true
    end
end
function proverka(r)
    for side in (Ost,West,Sud,Nord)
        if try_move(r,side) == false
            return false
        else
            return true
        end
    end
end
next(side::HorizonSide)=HorizonSide(mod(Int(side)+1,4))
function find_border(r)
    num_steps_max=1
    side=Nord
    while isborder(r, side)==false
        for _ in 1:2
            fborder(r,side,num_steps_max)
            side=next(side)
        end
        num_steps_max+=1
    end
end

function fborder(r,side,num_steps_max)
    for _ in 1:num_steps_max
        if isborder(r,side) && proverka(r) == false
            return nothing
        end
        move!(r,side)
    end
end


text:
using BlackBoxOptim

# Let's try to solve some TSP problems with continuous optimization
# algorithms by mapping from the float candidate vector to a TSP
# solution.
include("read_tsplib_file.jl")

TspProblemFile = "dantzig42.tsp" # Best known solution is 699
TspProblemFile = "fri26.tsp" # Best known solution is 937
const TSP = read_tsplib_file(TspProblemFile)

# We get inspiration from the mapping methods in section 2.4 of the paper:
# https://www.sciencedirect.com/science/article/pii/S2210650219304468?casa_token=Pyu6JJAzBQoAAAAA:J-EYa_kdwfDQidRAHbxmvnAwsmttUWO-hDESjTQyLlr70RMOfYdOijJs8Reao1yr-pXrIgpuxsY

abstract type FloatVectorMapping end
abstract type PermutationMapping <: FloatVectorMapping end

# Some mappings will use the best so far as a guide but by default not.
updatefitness!(p::PermutationMapping, fitness::Float64) = nothing

# This has also been called SPV (Smallest Position Value) but they
# look the same.
struct RankedOrderValueMapping <: PermutationMapping end

function apply(m::RankedOrderValueMapping, v::Vector{Float64})
    rov = zeros(Int, length(v))
    rov[sortperm(v)] = 1:length(v)
    rov
end
# v = [2.54, 1.52, 2.13, 3.48, 2.37, 2.87]
# @assert apply(RankedOrderValueMapping(), v) == [4, 1, 2, 6, 3, 5]

# sortperm right away should be the same though but faster:
struct SortPermMapping <: PermutationMapping
    reverse::Bool
    SortPermMapping(r::Bool = false) = new(r)
end
apply(m::SortPermMapping, v::Vector{Float64}) = sortperm(v, rev=m.reverse)

# The literature on mappings is strange with some papers
# basically just reversing the order and claiming novelty.
# For example: Li, Xiangtao, and Minghao Yin. "A hybrid cuckoo
# search via Lévy flights for the permutation flow shop scheduling
# problem." International Journal of Production Research 51.16 (2013): 4732-4754.
# Proposes Largest Ranked Value which just seems to be the reverse
# order of SPV/ROV above. I can't imagine this can have any effect
# on optimization performance.
struct LargestRankedValueMapping <: PermutationMapping end
function apply(m::LargestRankedValueMapping, v::Vector{Float64})
    rov = zeros(Int, length(v))
    rov[sortperm(v, rev=true)] = 1:length(v)
    rov
end

# Let's try something of our own, i.e. truncate/floor countinuous
# values to integers, then handle the duplicates by assigning the
# remaining values based on which of the remaining ones it is closest
# to.
struct TruncateThenClosestRemainingMapping <: PermutationMapping
    N::Int
    s::Set{Int}
    TruncateThenClosestRemainingMapping(N::Int) = new(N, Set(1:N))
end
searchrange(m::TruncateThenClosestRemainingMapping) = (1.0, m.N + 1.0 - 1e-10)

function apply(m::TruncateThenClosestRemainingMapping, v::Vector{Float64})
    ints = floor.(Int, v)
    # Now repair by selecting the closest unselected
    selected = unique(ints)
    unselected = collect(setdiff(m.s, selected))
    counts = zeros(Int, m.N)
    for (idx, ival) in enumerate(ints)
        if counts[ival] > 0
            _, i = findmin(map(r -> abs(ival - r), unselected))
            r = unselected[i]
            counts[r] += 1
            ints[idx] = r
            deleteat!(unselected, i)
        else
            counts[ival] += 1
        end
    end
    ints
end

const ROV = RankedOrderValueMapping()
const SPM = SortPermMapping()
const LRV = LargestRankedValueMapping()
const TRC = TruncateThenClosestRemainingMapping(size(TSP))

function make_perm_mapping_fitness(tsp, pm::PermutationMapping)
    function tspfitnessfn(v::Vector{Float64})
        c = cost(tsp, apply(pm, v))
        updatefitness!(pm, c)
        c
    end
    tspfitnessfn
end

res_rov = bboptimize(make_perm_mapping_fitness(TSP, ROV);
    SearchRange = (0.0, 10.0), NumDimensions = size(TSP),
    PopulationSize = 1000,
    MaxTime = 30.0)
# dantzig42.tsp best runs: 836, 818
# fri26.tsp best runs: 977

res_spm = bboptimize(make_perm_mapping_fitness(TSP, SPM);
    SearchRange = (0.0, 10.0), NumDimensions = size(TSP),
    PopulationSize = 1000,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 749, 753
# fri26.tsp best runs: 937

res_spm = bboptimize(make_perm_mapping_fitness(TSP, SPM);
    SearchRange = (0.0, 10.0), NumDimensions = size(TSP),
    PopulationSize = 5000,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 753

res_lrv = bboptimize(make_perm_mapping_fitness(TSP, LRV);
    SearchRange = (0.0, 10.0), NumDimensions = size(TSP),
    PopulationSize = 5000,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 881

res_spm = bboptimize(make_perm_mapping_fitness(TSP, SPM);
    SearchRange = (0.0, 10.0), NumDimensions = size(TSP),
    Method = :xnes,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 853

res_trc1000 = bboptimize(make_perm_mapping_fitness(TSP, TRC);
    SearchRange = searchrange(TRC), NumDimensions = size(TSP),
    PopulationSize = 1000,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 699, 722, 699, 699
# fri26.tsp best runs: 937

res_trc100 = bboptimize(make_perm_mapping_fitness(TSP, TRC);
    SearchRange = searchrange(TRC), NumDimensions = size(TSP),
    PopulationSize = 100,
    MaxTime = 30.0)
# dantzig42.tsp Best runs: 737, 778

# The TRC mapping is good on this problem since it "repairs" duplicates
# by turning them into the closest remaining city in the natural encoding.

# An interesting representation might be to encode the probability
# to take the road from one city to another. This scales as O(n^2) but would
# enable global learning since the same position in the genome means the same
# thing for all individuals. We only add the starting city as a rounded float
# and then use the highest probability of the already non-visited nodes
# from then on to create a unique TSP traversal.
struct EdgeTraversalProbabilitiesMapping <: PermutationMapping
    N::Int
    EdgeTraversalProbabilitiesMapping(N::Int) = new(N)
end
function searchspace(m::EdgeTraversalProbabilitiesMapping)
    ss = [(1.0, float(m.N)+0.99999999)]
    for i in 1:(m.N^2)
        push!(ss, (0.0, 1.0))
    end
    ss
end

apply(m::EdgeTraversalProbabilitiesMapping, v::Vector{Float64}) =
    apply!(m, v, Array{Int}(undef, m.N))

function apply!(m::EdgeTraversalProbabilitiesMapping,
    v::Vector{Float64}, permutation::Vector{Int})

    i = 1
    permutation[i] = currentcity = floor(Int, v[1])
    visited = BitArray{1}(zeros(m.N))
    visited[currentcity] = 1
    while i < m.N
        i += 1
        startidx = 2 + (currentcity-1)*m.N
        endidx = startidx + m.N - 1
        perm = sortperm(view(v, startidx:endidx), rev=true)
        for candidatecity in perm
            if candidatecity != currentcity && !visited[candidatecity]
                permutation[i] = candidatecity
                visited[Int(candidatecity)] = 1
                currentcity = candidatecity
                break
            end
        end
    end
    permutation
end

const ETP = EdgeTraversalProbabilitiesMapping(size(TSP))

# Doesn't seem to work well, I guess the scaling makes it costly.
res_etp1000 = bboptimize(make_perm_mapping_fitness(TSP, ETP);
    SearchSpace = searchspace(ETP),
    PopulationSize = 1000,
    MaxTime = 60.0)

# We could try TRC with backups i.e. each also has a backup which is checked
# and only if that also is not available do we search for the closest remaining one.


text:
using REPL
using Terminals

ccall(:jl_install_sigint_handler, Void, ())

(quiet,repl,startup,color_set,history) = Base.process_options(ARGS)

if repl
    if !isa(STDIN,Base.TTY)
        if !color_set
            eval(Base, :(have_color = false))
        end
        # note: currently IOStream is used for file STDIN
        if isa(STDIN,File) || isa(STDIN,IOStream)
            # reading from a file, behave like include
            eval(Base, :(is_interactive = false))
            eval(parse_input_line(readall(STDIN)))
        else
            # otherwise behave repl-like
            eval(Base, :(is_interactive = true))
            while !eof(STDIN)
                eval_user_input(parse_input_line(STDIN), true)
            end
        end
    else
        t = Terminals.Unix.UnixTerminal(get(ENV,"TERM",""),STDIN,STDOUT,STDERR)

        if !color_set
            eval(Base,:(have_color = $(hascolor(t))))
        end

        eval(Base, :(is_interactive = true))
        REPL.run_repl(t)
    end
end
if Base.have_color
    print(Base.color_normal)
end
quit()


text:
using Revise
using Random
using Printf
using StatsBase
using StaticArrays
using Distributions
using LinearAlgebra
using POMDPs
using AutoViz
using AutomotiveDrivingModels
using AutomotivePOMDPs
using POMDPSimulators
using BeliefUpdaters
using POMDPPolicies
using RLInterface
using POMDPModelTools

# LOGGER = Logger("log1", overwrite=true)
includet("environment.jl")
includet("generative_mdp.jl")
# includet("masking.jl")
includet("cooperative_IDM.jl")
includet("belief_updater.jl")
includet("belief_mdp.jl")
includet("overlays.jl")
includet("make_gif.jl")

rng = MersenneTwister(1)

mdp = GenerativeMergingMDP(random_n_cars=true, min_cars=12,max_cars=14, traffic_speed = :mixed, driver_type=:random,
                           observe_cooperation=true)

s0 = initialstate(mdp, rng)

scene = Scene()
vehstate = vehicle_state(38.0, merge_lane(mdp.env), 10.0, mdp.env.roadway)
push!(scene, Vehicle(vehstate, VehicleDef(), EGO_ID))
vehstate = vehicle_state(80.0, main_lane(mdp.env), 10.0, mdp.env.roadway)
push!(scene, Vehicle(vehstate, VehicleDef(), 2))
vehstate = vehicle_state(107.0, main_lane(mdp.env), 10.0, mdp.env.roadway)
push!(scene, Vehicle(vehstate, VehicleDef(), 3))

c = AutoViz.render(scene, mdp.env.roadway,
            SceneOverlay[
                #  DistToMergeOverlay(target_id=EGO_ID, env=mdp.env, line_width=0.3),
                #  MergingNeighborsOverlay(target_id=EGO_ID, env=mdp.env)
            ],
            cam = StaticCamera(VecE2(-5.0, -5.0), 25.0),
            car_colors=get_car_type_colors(s0.scene, mdp.driver_models),
            canvas_width=1000,
            canvas_height=400)

using Cairo

write_to_png(c, "roadway.png")

c = AutoViz.render([ArrowCar([0.,0.], color=RGBA(0.8, 0.8, 0.8, 1))])

write_to_png(c, "ego_shade2.png")


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/GradientDescent.jl")
include("../data/data.jl")

X, y = getFullData()

lambda = 0.1
iter = 250

x0 = zeros(size(X,2),1)

f = f_logreg(X,y,x0,lambda)
g = g_logreg(X,y,x0,lambda)

@time (x, status,vals1,stops1 ) = qn(x0,f_logreg(X,y,x0,lambda),g_logreg(X,y,x0,lambda),1e-12,iter,20,20,1e-4,.9,"wolfe","bfgs")

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
# Iex
__precompile__()
module Iex

# Connections
include("connection/Connect.jl")
export Connection, status, set_connection_url!, set_connection_version!,
    set_connection_token!

# Accounts
include("account/Account.jl")
export usage, metadata

# Utils
include("stock/Utils.jl")
export get_stocks

# Stocks
include("stock/Prices.jl")
export history, batch, book, symbol_quote, price,
    ohlc, largest_trades, batched_closes

# Profiles
include("stock/Profiles.jl")
export company, insider_roster, insider_summary, insider_transactions,
    logo, peers

    # Fundamentals
include("stock/Fundamentals.jl")
export balance_sheet, dividends, earnings, income, split

# Research
include("stock/Research.jl")
export stats

# News
include("stock/News.jl")
export news

# Reference
include("reference/Reference.jl")
export symbols

# Options
include("stock/Options.jl")
export bsopm_call, bsopm_impvol, calc_term, calc_hr, fortify_mv_df!
export compute_mv_delta, lme_mv_delta

const fields_to_merge = [:id, :ipodate, :enddate, :dadate, :mergerdate, :closingdate]
export fields_to_merge

end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
@testset "Pipeline" begin
    @test eltype(Pipeline) <: ImageTransformation
end

@testset "LinearPipeline" begin
    @testset "concrete pipeline and image" begin
        pl = [FlipX(), FlipY(), CropRatio(1.5), Resize(64,64)]
        @imagetest "LinearPipeline" transform(pl, testimg)

        testimgs = Array{Image,1}()
        push!(testimgs, testimg)
        push!(testimgs, testimg)
        imgs = transform(pl, testimgs)
        @imagetest "LinearPipeline" imgs[1]
        @imagetest "LinearPipeline" imgs[2]
    end
end


text:
# ODE SIU/SICU models
# variables are concentrations, X/N
# Z = 1-S
# U = are noninfectious include both recovered and dead

function sicuprep(data)
    df = data
    df[:,:uninfected] = df[:,:cum_recover] + df[:,:cum_death]
    return df
end

# SICUq ODE
# Cases are quarantined with effectiveness q
# 3+1 parameters
function sicu!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma = p
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I)*(1-Z) - sigmac*I - sigma*I
	du[4] = dZ = beta*(I)*(1-Z)
end

# SICUq ODE
# Cases are quarantined with effectiveness q
# 4+1 parameters
function sicuq!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,q = p
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I+q*C)*(1-Z) - sigmac*I - sigma*I
	du[4] = dZ = beta*(I+q*C)*(1-Z)
end

# SICUf ODE
# Cases are quarantined with factor q
# Cases recover or die with factor f
# 4+1 parameters
function sicuf!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,f = p
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I)*(1-Z) - sigmac*I - f*sigma*I
	du[4] = dZ = beta*(I)*(1-Z)
end

# SICUqf ODE
# Cases are quarantined with factor q
# Cases recover or die with factor f
# 5+1 parameters
function sicuqf!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,q,f = p
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I+q*C)*(1-Z) - sigmac*I - f*sigma*I
	du[4] = dZ = beta*(I+q*C)*(1-Z)
end

# Mitigation applied
# 5 + 1 parameters
function sicum!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,mc,ml = p
    beta *= mc + (1-mc)/(1 + exp(t - 35. - ml))
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I)*(1-Z) - sigmac*I - sigma*I
	du[4] = dZ = beta*(I)*(1-Z)
end

# Mitigation applied after C > 500
# 6+1
function sicufm!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,f,mc,ml = p
    beta *= mc + (1-mc)/(1 + exp(t - 35. - ml))
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I)*(1-Z) - sigmac*I - f*sigma*I
	du[4] = dZ = beta*(I)*(1-Z)
end

# Mitigation applied after C > 500
# 6+1
function sicuqm!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,q,mc,ml = p
    beta *= mc + (1-mc)/(1 + exp(t - 35. - ml))
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I+q*C)*(1-Z) - sigmac*I - sigma*I
	du[4] = dZ = beta*(I+q*C)*(1-Z)
end

# Mitigation applied after C > 500
# 7+1
function sicuqfm!(du,u,p,t)
	C,U,I,Z = u
	beta,sigmac,sigma,q,f,mc,ml = p
    beta *= mc + (1-mc)/(1 + exp(t - 35. - ml))
	du[1] = dC = sigmac*I - sigma*C
	du[2] = dU = sigma*C
	du[3] = dI = beta*(I+q*C)*(1-Z) - sigmac*I - f*sigma*I
	du[4] = dZ = beta*(I+q*C)*(1-Z)
end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
T=6
tsname="joLinearOperator - Product"
@testset "$tsname" begin
for t=1:T # start test loop

if t<=2
    (m,n)=(5,5)
else
    (m,n)=(rand(3:7),rand(3:7))
end

if t%2==1
    tname="loop $t for real($m,$n)"
    a=rand(m,n)
    b=rand(m,n)
    vn=rand(n)
    mvn=rand(n,2)
    vm=rand(m)
    mvm=rand(m,2)
    A=joMatrix(a)
    B=joLinearFunctionAll(m,n,v->b*v,v->transpose(b)*v,v->adjoint(b)*v,v->conj(b)*v,v->b\v,v->transpose(b)\v,v->adjoint(b)\v,v->conj(b)\v,eltype(vn),eltype(vm);fMVok=true,iMVok=true)
else
    tname="loop $t for complex($m,$n)"
    a=rand(ComplexF64,m,n)
    b=rand(ComplexF64,m,n)
    vn=rand(ComplexF64,n)
    mvn=rand(ComplexF64,n,2)
    vm=rand(ComplexF64,m)
    mvm=rand(ComplexF64,m,2)
    A=joLinearFunctionAll(m,n,v->a*v,v->transpose(a)*v,v->adjoint(a)*v,v->conj(a)*v,v->a\v,v->transpose(a)\v,v->adjoint(a)\v,v->conj(a)\v,eltype(vn),eltype(vm);fMVok=true,iMVok=true)
    B=joMatrix(b)
end
c=a*adjoint(b)
C=A*adjoint(B)
afac=rand(eltype(c))
Afac=joNumber(afac,C)
mfac=rand(eltype(c))
Mfac=joNumber(mfac,C)

verbose && println("$tsname $tname")
    @testset "$tname A*B" begin
        @test size(C)==size(c)
        @test length(C)==length(c)
        @test norm(jo_full(C)-c)<joTol
        @test norm(elements(C)-c)<joTol
        @test norm(C)-norm(c)<joTol
        for i=1:2
            @test norm(C,i)-norm(c,i)<joTol
        end
        @test norm(C,Inf)-norm(c,Inf)<joTol
        @test norm(C)-norm(c)<joTol
        for i=1:2
            @test norm(C,i)-norm(c,i)<joTol
        end
        @test norm(C,Inf)-norm(c,Inf)<joTol
        @test norm(elements(adjoint(C))-adjoint(c))<joTol
        @test norm(elements(transpose(C))-transpose(c))<joTol
        @test norm(elements(+C)-(+c))<joTol
        @test norm(elements(-C)-(-c))<joTol
        @test norm(C*vm-c*vm)<joTol
        @test norm(C*mvm-c*mvm)<joTol
        @test norm(adjoint(C)*vm-adjoint(c)*vm)<joTol
        @test norm(adjoint(C)*mvm-adjoint(c)*mvm)<joTol
        @test norm(transpose(C)*vm-transpose(c)*vm)<joTol
        @test norm(transpose(C)*mvm-transpose(c)*mvm)<joTol
        @test norm(adjoint(transpose(C))*vm-adjoint(transpose(c))*vm)<joTol
        @test norm(adjoint(transpose(C))*mvm-adjoint(transpose(c))*mvm)<joTol
        @test norm(transpose(adjoint(C))*vm-transpose(adjoint(c))*vm)<joTol
        @test norm(transpose(adjoint(C))*mvm-transpose(adjoint(c))*mvm)<joTol
        @test norm(elements(transpose(C)*C)-(transpose(c)*c))<joTol
        @test norm(elements(C*transpose(C))-(c*transpose(c)))<joTol
        @test norm(elements(adjoint(C)*C)-(adjoint(c)*c))<joTol
        @test norm(elements(C*adjoint(C))-(c*adjoint(c)))<joTol
        @test norm(elements(mfac*C)-(mfac*c))<joTol
        @test norm(elements(mfac*transpose(C))-(mfac*transpose(c)))<joTol
        @test norm(elements(mfac*adjoint(C))-(mfac*adjoint(c)))<joTol
        @test norm(elements(mfac*conj(C))-(mfac*conj(c)))<joTol
        @test norm(elements(C*mfac)-(c*mfac))<joTol
        @test norm(elements(transpose(C)*mfac)-(transpose(c)*mfac))<joTol
        @test norm(elements(adjoint(C)*mfac)-(adjoint(c)*mfac))<joTol
        @test norm(elements(conj(C)*mfac)-(conj(c)*mfac))<joTol
        @test norm(elements(Mfac*C)-(mfac*c))<joTol
        @test norm(elements(Mfac*transpose(C))-(mfac*transpose(c)))<joTol
        @test norm(elements(Mfac*adjoint(C))-(mfac*adjoint(c)))<joTol
        @test norm(elements(Mfac*conj(C))-(mfac*conj(c)))<joTol
        @test norm(elements(C*Mfac)-(c*mfac))<joTol
        @test norm(elements(transpose(C)*Mfac)-(transpose(c)*mfac))<joTol
        @test norm(elements(adjoint(C)*Mfac)-(adjoint(c)*mfac))<joTol
        @test norm(elements(conj(C)*Mfac)-(conj(c)*mfac))<joTol
        @test norm(elements(C+afac)-(c.+afac))<joTol
        @test norm(elements(transpose(C)+afac)-(transpose(c).+afac))<joTol
        @test norm(elements(adjoint(C)+afac)-(adjoint(c).+afac))<joTol
        @test norm(elements(conj(C)+afac)-(conj(c).+afac))<joTol
        @test norm(elements(afac+C)-(afac.+c))<joTol
        @test norm(elements(afac+transpose(C))-(afac.+transpose(c)))<joTol
        @test norm(elements(afac+adjoint(C))-(afac.+adjoint(c)))<joTol
        @test norm(elements(afac+conj(C))-(afac.+conj(c)))<joTol
        @test norm(elements(C+Afac)-(c.+afac))<joTol
        @test norm(elements(transpose(C)+Afac)-(transpose(c).+afac))<joTol
        @test norm(elements(adjoint(C)+Afac)-(adjoint(c).+afac))<joTol
        @test norm(elements(conj(C)+Afac)-(conj(c).+afac))<joTol
        @test norm(elements(Afac+C)-(afac.+c))<joTol
        @test norm(elements(Afac+transpose(C))-(afac.+transpose(c)))<joTol
        @test norm(elements(Afac+adjoint(C))-(afac.+adjoint(c)))<joTol
        @test norm(elements(Afac+conj(C))-(afac.+conj(c)))<joTol
        if (issquare(C) && m==n)
            @test norm(C\vm-c\vm)<joTol^(2/3)
            @test norm(C\mvm-c\mvm)<joTol^(2/3)
            @test norm(adjoint(C)\vn-adjoint(c)\vn)<joTol^(2/3)
            @test norm(adjoint(C)\mvn-adjoint(c)\mvn)<joTol^(2/3)
            @test norm(transpose(C)\vn-transpose(c)\vn)<joTol^(2/3)
            @test norm(transpose(C)\mvn-transpose(c)\mvn)<joTol^(2/3)
        end
    end

end # end test loop
end


text:
const initialObjValue = -787.0936089997643
const objGap = -0.005*initialObjValue

const riskLowerBound = 5E-04

user_initial_guess = Float64[9.885150021498452 12.345440853151693 8.09804788503498 8.222368074757982 6.3193792865228815 6.718181558496166 7.002166961525422 9.830082624476422 8.539306613084479 8.484855963002648 9.426512704571412 9.516220380584562 7.744186023270856 7.391250853029938 10.394023349415885 10.337654518499786 10.97371413162148 10.878109402087683 8.866970927667884 7.408718801021087 10.300793554219874 12.851671633620022 6.690563099618595 7.878893519026003 5.120913794479996 5.410516080037175 15.778854580767362 11.098267189826315 8.835172095205632 9.553025219588923 8.812089746955834 6.714709214557973 9.534719877178526 7.496042408525231 7.147518921513608 13.648614981894145 7.260634989597417 11.456064390055559 6.052902513583182 10.38812763888143 7.974679881251275 12.609922725925337 7.332591436859364 9.282090116430332 10.423554048557984 8.6365161818614 6.6563134473492065 7.9620545754711785 7.760637514077624 7.340893985975288 6.1273745906303 8.777373054803983 9.756198983370405 5.252333354077139 9.362159429451127 7.9283726039851725 8.929033522853386 5.605079715029882 10.164748098627825 8.273980272880728 5.541245266355158 10.46737031629924 8.373117753289959 8.141746327851461 5.159213859166618 6.569206783607733 7.070161983412617 7.36773096283932 12.203568372581522 5.447947127042287 5.798216811196845 10.632394919667618 5.945806275939606 4.8214313465965954 9.1198821733236 6.195131038873298 3.9639323549948844 8.309316080188168 7.800418371991007 6.410109660535694 5.8270294087652434 3.804621844846956 4.862054184991828 8.072274207443543 6.3048233792870505 5.889459966371539 5.322580485871152 5.157935354858375 6.82591176323608 3.622878989501982 4.595768695801355 8.550473558613746 3.878678798810665 4.599766791623326 7.305491078847105 3.4502452801710857 8.87967213705522 5.392859227818847 7.115728686842548 5.729063691432691]
user_initial_guess = user_initial_guess'

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
@testset "data.jl" begin
    @testset "dataset(name)" begin
        # nonexistent file throws
        @test_throws ErrorException RatingCurves.dataset("nonexistant")

        # 50408 loading
        df = RatingCurves.dataset("50408")
        @test size(df, 1) == 194
        @test size(df, 2) == 3
    end

end


text:
################################################################################
# Regularizer
################################################################################

@with_kw mutable struct Regularizer{T}

	# Primal regularization
	"Jacobian regularization for primal variables."
	x::T=1e-3
	u::T=1e-3

	# Dual regularization
	"Jacobian regularization for dual variables."
	λ::T=1e-3
end

import Base.map!
function map!(reg::Regularizer, f)
	for name in fieldnames(Regularizer)
		setfield!(reg, name, f(getfield(reg, name)))
	end
	return nothing
end

function set!(reg::Regularizer, val)
	f(x) = val
	map!(reg, f)
	return nothing
end

function mult!(reg::Regularizer, val)
	f(x) = x*val
	map!(reg, f)
	return nothing
end


text:
export Categorical

"""
Description:

    Categorical factor node

    The categorical node defines a one-dimensional probability
    distribution over the normal basis vectors of dimension d

    out ∈ {0, 1}^d where Σ_k out_k = 1
    p ∈ [0, 1]^d, where Σ_k p_k = 1

    f(out, p) = Cat(out | p)
              = Π_i p_i^{out_i}

Interfaces:

    1. out
    2. p

Construction:

    Categorical(id=:some_id)
"""
mutable struct Categorical <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Categorical(out, p; id=generateId(Categorical))
        @ensureVariables(out, p)
        self = new(id, Array{Interface}(undef, 2), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:p] = self.interfaces[2] = associate!(Interface(self), p)

        return self
    end
end

slug(::Type{Categorical}) = "Cat"

format(dist::ProbabilityDistribution{Univariate, Categorical}) = "$(slug(Categorical))(p=$(format(dist.params[:p])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{Categorical}; p=[1/3, 1/3, 1/3]) = ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>p))
ProbabilityDistribution(::Type{Categorical}; p=[1/3, 1/3, 1/3]) = ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>p))

dims(dist::ProbabilityDistribution{Univariate, Categorical}) = ()

vague(::Type{Categorical}, n_factors::Int64=3) = ProbabilityDistribution(Univariate, Categorical, p=(1/n_factors)*ones(n_factors))
vague(::Type{Categorical}, n_factors::Tuple) = vague(Categorical, n_factors[1])

isProper(dist::ProbabilityDistribution{Univariate, Categorical}) = (abs(sum(dist.params[:p])-1.) < 1e-6)

unsafeMean(dist::ProbabilityDistribution{Univariate, Categorical}) = deepcopy(dist.params[:p])

unsafeMeanVector(dist::ProbabilityDistribution{Univariate, Categorical}) = deepcopy(dist.params[:p])

function unsafeMode(dist::ProbabilityDistribution{Univariate, Categorical})
    i = findfirst(dist.params[:p] .== maximum(dist.params[:p])) # Index of first maximum
    m = zeros(length(dist.params[:p]))
    m[i] = 1.0

    return m
end

logPdf(dist::ProbabilityDistribution{Univariate, Categorical}, x) = sum(x .* log.(dist.params[:p]))

function sample(dist::ProbabilityDistribution{Univariate, Categorical})
    isProper(dist) || error("Cannot sample from improper distribution $dist")

    p = dist.params[:p]
    p_sum = cumsum(p);
    z = rand()

    d = length(p)
    idx = 0
    for i = 1:d
        if z < p_sum[i]
            idx = i
            break
        end
        idx = d
    end

    x = spzeros(Float64, d)
    x[idx] = 1.0

    return x
end

naturalParams(dist::ProbabilityDistribution{Univariate, Categorical}) = vcat(log.(dist.params[:p][1:end-1]./dist.params[:p][end]), 0.0) # Variant 3 of https://en.wikipedia.org/wiki/Exponential_family

standardDistribution(V::Type{Univariate}, F::Type{Categorical}; η::Vector) = ProbabilityDistribution(V, F, p=exp.(η)./sum(exp.(η)))

logNormalizer(::Type{Univariate}, ::Type{Categorical}; η::Vector) = log(sum(exp.(η)))

logPdf(V::Type{Univariate}, F::Type{Categorical}, x::Vector; η::Vector) = x'*η - logNormalizer(V, F; η=η)

function prod!( x::ProbabilityDistribution{Univariate, Categorical},
                y::ProbabilityDistribution{Univariate, Categorical},
                z::ProbabilityDistribution{Univariate, Categorical}=ProbabilityDistribution(Univariate, Categorical, p=ones(size(x.params[:p]))./length(x.params[:p])))

    # Multiplication of 2 categorical PMFs: p(z) = p(x) * p(y)
    z.params[:p][:] = clamp.(x.params[:p] .* y.params[:p], tiny, Inf) # Soften vanishing probabilities to enforce the convention [1, 0]*[0, 1] ∝ [1/2, 1/2]
    norm = sum(z.params[:p])
    z.params[:p] = z.params[:p]./norm

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{Univariate, Categorical},
                            y::ProbabilityDistribution{Multivariate, PointMass},
                            z::ProbabilityDistribution{Multivariate, PointMass}=ProbabilityDistribution(Multivariate, PointMass, m=[0.0]))

    z.params[:m] = deepcopy(y.params[:m])

    return z
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, Categorical})
    -sum(dist.params[:p].*log.(clamp.(dist.params[:p], tiny, Inf))) # Soften vanishing propbabilities to enforce the convention 0 log 0 = 0
end

# Average energy functional
function averageEnergy(::Type{Categorical}, marg_out::ProbabilityDistribution, marg_p::ProbabilityDistribution{Multivariate})
    -sum(unsafeMean(marg_out).*unsafeLogMean(marg_p))
end

function averageEnergy(::Type{Categorical}, marg_out::ProbabilityDistribution{Multivariate,PointMass}, marg_p::ProbabilityDistribution{Multivariate})
    -sum(marg_out.params[:m].*unsafeLogMean(marg_p))
end


text:
using Restful
import Restful.json

const _web = Restful.app()
const extra_header = Ref{String}("")

_web.get("/") do req, res, route
    theme = get(ENV, "WEB_DISPLAY_THEME", "light")
    res.html(extra_header[] * read(joinpath(@__DIR__, "..", "assets", "$theme.html"), String))
end

_web.get("/hist", json) do req, res, route
    from = try parse(Int, req.params["from"]) catch; 1 end
    from > length(_display.hist) && wait(_display.cond)

    res.json(map(first, _display.hist[from:end]))
end

_web.get("/hist/:id") do req, res, route
    id = parse(Int, route.id)
    id > length(_display.hist) && return res.code(404)

    res.status = 200
    res.content_type, res.body = _display.hist[id]
end

_web.delete("/hist/:id") do req, res, route
    id = parse(Int, route.id)
    id > length(_display.hist) && return res.code(404)

    _display.hist[id] = ("", UInt8[])
    res.code(200)
end


text:
module Match
using MLStyle.Err
using MLStyle.Private

export @match, Pattern, Case, Failed, failed, PatternDef, pattern_match, app_pattern_match, register_app_pattern, register_meta_pattern

struct Failed end
failed = Failed()

struct Pattern
    shape
    guard
end

struct Case
    patterns :: Vector{Pattern}
    body
end

# """
# e.g:
# register_meta_pattern((expr :: Expr) -> <conditional>) do  expr, guard, tag, mod
#   <return quote>
# end
# """
global meta_dispatchers = Vector{Tuple{Function, Function}}()

# """
# e.g
#
# register_app_pattern(ty) do args, guard, tag, mod
#     <return quote>
# end
# where `ty` is constructor or any other custom pattern(see @pattern)
# """
global app_dispatchers = Dict{Any, Function}()

global _count = 0
function mangling(apply, repr_ast)
    global _count = _count + 1
    let symbol = Symbol("<", repr(repr_ast), ".", _count, ">")
        ret = apply(symbol)
        global _count = _count - 1
        ret
    end
end

@inline function collect_guard(arg)
    if isa(arg, Expr) && arg.head == :curly
    let args = arg.args
        Pattern(args[1],
                Expr(:block,
                     reduce( (a, b) -> Expr(:&&, a, b), args[2:end])
                     ))
    end
    else
        Pattern(arg, nothing)
    end
end

function collect_fallthrough(c :: Channel, arg)

    if isa(arg, Expr) && arg.head == :call && arg.args[1] == :|

        let nested = arg.args[2:end-1], last = collect_guard(arg.args[end])

            if isempty(nested)
                put!(c, last)
            else
                foreach(nested) do n
                    collect_fallthrough(c, n)
                end
                put!(c, last)
            end

        end

    else
        put!(c, collect_guard(arg))
    end
end

function collect_fallthrough(arg)
    Channel(c -> collect_fallthrough(c, arg))
end

_warned = false
function merge_cases(cases :: Vector{Case})

    global _warned
    if !_warned
        @warn "Pattern optimization not implemented yet."
        _warned = true
    end

    cases
end

pattern_match_maker(tag :: Symbol, case :: Case, mod :: Module) =
    if isempty(case.patterns)
        InternalException("Internal Error 0") |> throw
    else let body = case.body,
             cond = map(case.patterns) do pattern
                    pattern_match(pattern.shape, pattern.guard, tag, mod)
                    end |>
                    function (last)
                        reduce((a, b) -> Expr(:||, a, b), last)
                    end
        quote
            if $cond
                $(case.body)
            else
                $failed
            end
        end
        end
    end

function make_dispatch(tag, cases, mod :: Module)
    let cases =
        map(cases) do case
                if !(case.head == :call && case.args[1] == :(=>))
                    SyntaxError("A case should be something like <pattern> => <body>") |> throw
                end

                let args = case.args, pattern = args[2], body = args[3]
                    Case(append!([], collect_fallthrough(pattern)), body)
                end
        end

    map(it -> pattern_match_maker(tag, it, mod), cases)
    end
end

macro match(target, pattern_def)

    if !isa(pattern_def, Expr) || !(pattern_def.head in (:braces, :bracescat, :block))
        SyntaxError("Invalid leading marker of match body.") |> throw
    end

    mangling(target) do tag_sym
        final =
            quote
                throw(($InternalException)("Non-exhaustive pattern found!"))
            end

        mangled = Symbol("case", ".", "test")

        args =
            if pattern_def.head === :block
                filter(it -> !isa(it, LineNumberNode), pattern_def.args)
            else
                pattern_def.args
            end
        for dispatched in make_dispatch(tag_sym, args[end:-1:1], __module__)
            final =
                quote
                    let $mangled =
                        let
                            $dispatched
                        end
                    if $mangled === $failed
                        $final
                    else
                        $mangled
                    end
                    end
                end
        end
        quote
              let $tag_sym = $target
                $final
              end
        end |> esc
        end
end

function pattern_match(expr :: Expr, guard, tag, mod :: Module)

    for (dispatcher_test, dispatch) in meta_dispatchers
        if dispatcher_test(expr)
            return dispatch(expr, guard, tag, mod)
        end
    end
    PatternUnsolvedException(expr) |> throw

end

function app_pattern_match(destructor, args, guard, tag, mod :: Module)
    global app_dispatchers
    fn =
        let destructor = get_most_union_all(destructor, mod)
                get(app_dispatchers, destructor, nothing)
        end

    if fn === nothing
        PatternUnsolvedException(:($destructor($(args...)))) |> throw
    else
        fn(args, guard, tag, mod)
    end

end

function register_meta_pattern(dispatch :: Function, dispatcher_test :: Function)
    global meta_dispatchers
    push!(meta_dispatchers, (dispatcher_test, dispatch));
end

function register_app_pattern(dispatch :: Function, destructor :: Any)
    global app_dispatchers
    push!(app_dispatchers, (destructor => dispatch))
end

# """
# builtin meta pattern(for application pattern)
# """
register_meta_pattern(expr :: Expr -> expr.head === :call) do expr, guard, tag, mod
    destructor = expr.args[1]
    app_pattern_match(destructor, expr.args[2:end], guard, tag, mod)
end

# """
# builtin application pattern(for dictionaries)
#
# x = "9"
# dict = Dict(1 => Dict(x => 3), 2 => 3)
#
# @match dict begin
#     Dict(1 => Dict(&x => a), 2 =>b) => a == b
# end
#
# # => true
#
# """
register_app_pattern(Dict) do args, guard, tag, mod

    matching =
        map(args) do kv
            if !(isa(kv, Expr) && kv.head === :call && (@eval mod $(kv.args[1])) === Pair)
                SyntaxError("Dictionary destruct must take patterns like Dict(<expr> => <pattern>)") |> throw
            end
            let (k, v) = kv.args[2:end]
                mangling(tag) do tag!
                    let tag! = Symbol(tag!, "[", k, "]"),
                        action = pattern_match(v, nothing, tag!, mod)

                        quote
                            $tag! = $get($tag, $k) do
                                    $failed
                            end
                            if $failed !== $tag!
                                $action
                            else
                                false
                            end
                        end
                    end
                end
            end
        end |>
        function (last)
            reduce((a, b) -> Expr(:&&, a, b), last, init=:($isa($tag, $Dict)))
        end

    check_ty = :($isa($tag, $Dict))

    if guard === nothing
        quote
            $check_ty && $matching
        end
    else
        quote
            $check_ty && $matching && $guard
        end
    end
end

register_app_pattern(in) do args, guard, tag, mod

    pattern, name = args

    let pat1 = pattern_match(pattern, nothing, tag, mod),
        pat2 = quote
                  $name = $tag
                  true
              end

        :($pat1 && $pat2) |>
        function (last)
            if guard === nothing
                last
            else
                :($last && $guard)
            end
        end
    end
end

module PatternDef
    import MLStyle.Match: register_app_pattern, register_meta_pattern

    App  = register_app_pattern
    Meta = register_meta_pattern
end

end # module


text:
using KitBase, FluxReconstruction, LinearAlgebra, OrdinaryDiffEq
using ProgressMeter: @showprogress
using Plots

pyplot()
cd(@__DIR__)

begin
    set = Setup(
        "gas",
        "cylinder",
        "2d0f",
        "hll",
        "nothing",
        1, # species
        3, # order of accuracy
        "positivity", # limiter
        "euler",
        0.1, # cfl
        1.0, # time
    )

    ps0 = KitBase.CSpace2D(1.0, 6.0, 60, 0.0, π, 50)
    deg = set.interpOrder - 1
    ps = FRPSpace2D(ps0, deg)

    vs = nothing
    gas = Gas(
        1e-6,
        2.0, # Mach
        1.0,
        1.0, # K
        5 / 3,
        0.81,
        1.0,
        0.5,
    )
    ib = nothing

    ks = SolverSet(set, ps0, vs, gas, ib)
end

n1 = [[0.0, 0.0] for i = 1:ps.nr+1, j = 1:ps.nθ]
for i = 1:ps.nr+1, j = 1:ps.nθ
    angle = sum(ps.dθ[1, 1:j-1]) + 0.5 * ps.dθ[1, j]
    n1[i, j] .= [cos(angle), sin(angle)]
end

n2 = [[0.0, 0.0] for i = 1:ps.nr, j = 1:ps.nθ+1]
for i = 1:ps.nr, j = 1:ps.nθ+1
    angle = π / 2 + sum(ps.dθ[1, 1:j-1])
    n2[i, j] .= [cos(angle), sin(angle)]
end

u0 = zeros(ps.nr, ps.nθ, deg + 1, deg + 1)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    u0[i, j, k, l] =
        max(exp(-3 * ((ps.xpg[i, j, k, l, 1] - 2)^2 + (ps.xpg[i, j, k, l, 2] - 2)^2)), 1e-2)
end

function dudt!(du, u, p, t)
    J, ll, lr, dhl, dhr, lpdm, ax, ay = p

    nx = size(u, 1)
    ny = size(u, 2)
    nsp = size(u, 3)

    f = zeros(nx, ny, nsp, nsp, 2)
    for i in axes(f, 1), j in axes(f, 2), k = 1:nsp, l = 1:nsp
        fg, gg = ax * u[i, j, k, l], ay * u[i, j, k, l]
        f[i, j, k, l, :] .= inv(J[i, j][k, l]) * [fg, gg]
    end

    u_face = zeros(nx, ny, 4, nsp)
    f_face = zeros(nx, ny, 4, nsp, 2)
    for i in axes(u_face, 1), j in axes(u_face, 2), l = 1:nsp
        u_face[i, j, 1, l] = dot(u[i, j, l, :], ll)
        u_face[i, j, 2, l] = dot(u[i, j, :, l], lr)
        u_face[i, j, 3, l] = dot(u[i, j, l, :], lr)
        u_face[i, j, 4, l] = dot(u[i, j, :, l], ll)

        for n = 1:2
            f_face[i, j, 1, l, n] = dot(f[i, j, l, :, n], ll)
            f_face[i, j, 2, l, n] = dot(f[i, j, :, l, n], lr)
            f_face[i, j, 3, l, n] = dot(f[i, j, l, :, n], lr)
            f_face[i, j, 4, l, n] = dot(f[i, j, :, l, n], ll)
        end
    end

    fx_interaction = zeros(nx + 1, ny, nsp)
    for i = 2:nx, j = 1:ny, k = 1:nsp
        au =
            (f_face[i, j, 4, k, 1] - f_face[i-1, j, 2, k, 1]) /
            (u_face[i, j, 4, k] - u_face[i-1, j, 2, k] + 1e-6)
        fx_interaction[i, j, k] = (
            0.5 * (f_face[i, j, 4, k, 1] + f_face[i-1, j, 2, k, 1]) -
            0.5 * abs(au) * (u_face[i, j, 4, k] - u_face[i-1, j, 2, k])
        )
    end
    fy_interaction = zeros(nx, ny + 1, nsp)
    for i = 1:nx, j = 2:ny, k = 1:nsp
        au =
            (f_face[i, j, 1, k, 2] - f_face[i, j-1, 3, k, 2]) /
            (u_face[i, j, 1, k] - u_face[i, j-1, 3, k] + 1e-6)
        fy_interaction[i, j, k] = (
            0.5 * (f_face[i, j, 1, k, 2] + f_face[i, j-1, 3, k, 2]) -
            0.5 * abs(au) * (u_face[i, j, 1, k] - u_face[i, j-1, 3, k])
        )
    end

    rhs1 = zeros(nx, ny, nsp, nsp)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp
        rhs1[i, j, k, l] = dot(f[i, j, :, l, 1], lpdm[k, :])
    end
    rhs2 = zeros(nx, ny, nsp, nsp)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp
        rhs2[i, j, k, l] = dot(f[i, j, k, :, 2], lpdm[l, :])
    end

    for i = 2:nx-1, j = 2:ny-1, k = 1:nsp, l = 1:nsp
        du[i, j, k, l] = -(
            rhs1[i, j, k, l] +
            rhs2[i, j, k, l] +
            (fx_interaction[i, j, l] - f_face[i, j, 4, l, 1]) * dhl[k] +
            (fx_interaction[i+1, j, l] - f_face[i, j, 2, l, 1]) * dhr[k] +
            (fy_interaction[i, j, k] - f_face[i, j, 1, k, 2]) * dhl[l] +
            (fy_interaction[i, j+1, k] - f_face[i, j, 3, k, 2]) * dhr[l]
        )
    end

    return nothing
end

tspan = (0.0, 1.0)
a = [1.0, 1.0]
p = (ps.J, ps.ll, ps.lr, ps.dhl, ps.dhr, ps.dl, a[1], a[2])
prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.01
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Tsit5(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:10#nt
    step!(itg)
end

contourf(ps.x, ps.y, itg.u[:, :, 2, 2], aspect_ratio = 1, legend = true)


text:
module DocumentationGenerator

using Pkg
using Pkg.TOML
using Pkg.Types
using Documenter

"""
    default_docs(package, root, pkgroot)

Generates a default documentation for a package without Documenter.jl docs.
"""
function default_docs(package, root, pkgroot)
    @info("Generating default fallback docs (readme + exported methods).")
    doc_source = joinpath(root, "src")
    mkpath(doc_source)
    pages = ["Docstrings" => "autodocs.md"]
    handle_readme(pages, pkgroot, doc_source)
    pkg_sym = Symbol(package)
    @eval Module() begin
        using Pkg
        Pkg.add("Documenter")
        using Documenter
        using $pkg_sym
        open(joinpath($doc_source, "autodocs.md"), "w") do io
            println(io, """
            ```@autodocs
            Modules = [$($pkg_sym)]
            ```
            """)
        end
        makedocs(
            format = Documenter.HTML(),
            sitename = "$($package).jl",
            modules = [$pkg_sym],
            root = $root,
            pages = $(reverse(pages))
        )
    end
end

function handle_readme(pages, pkgroot, doc_source)
    for file in readdir(pkgroot)
        if occursin("readme", lowercase(file))
            readme = joinpath(pkgroot, file)
            if isfile(readme)
                newreadmepath = joinpath(doc_source, "index.md")
                rendergfm(readme, newreadmepath)
                copylocallinks(readme, newreadmepath)
                push!(pages, "Readme" => "index.md")
                break
            end
        end
    end
end

function rendergfm(file, fileout)
    try
        commonmarker = find_ruby_gem("commonmarker")
        cmd = `$(commonmarker) $(file)`
        rendered = read(cmd, String)
        open(fileout, "w") do io
            # lots of backticks so the @raw block doesn't end prematurely
            println(io, "````````````@raw html\n", rendered, "\n````````````")
        end
    catch err
        cp(file, fileout)
        @error("Rendering GFM failed. Falling back to Julia implementation.", error = err)
    end
end

"""
    readme_docs(package, root, pkgroot)

Generates README based fallback docs when the package installs but can't be loaded.
"""
function readme_docs(package, root, pkgroot)
    @info("Generating readme-only fallback docs.")
    doc_source = joinpath(root, "src")
    mkpath(doc_source)
    pages = []
    handle_readme(pages, pkgroot, doc_source)
    pkg_sym = Symbol(package)
    @eval Module() begin
        using Pkg
        Pkg.add("Documenter")
        using Documenter
        makedocs(
            format = Documenter.HTML(),
            sitename = "$($package).jl",
            modules = [Module()],
            root = $root,
            pages = $(reverse(pages))
        )
    end
end

function find_ruby_gem(gem)
    which = Sys.iswindows() ? "where" : "which"
    so = IOBuffer()
    p = run(pipeline(`$which $gem`, stdout = so), wait = false)
    # gem is on path
    success(p) && return chomp(String(take!(so)))
    so = IOBuffer()
    p = run(pipeline(`$which gem`, stdout = so), wait = false)
    if !success(p)
        @error("`gem` not found on PATH")
        return ""
    end
    gempath = chomp(String(take!(so)))

    p = run(pipeline(`$gempath environment gempath`, stdout = so), wait = false)
    if success(p)
        paths = chomp(String(take!(so)))
        paths = split(paths, Sys.iswindows() ? ';' : ':')
        for path in paths
            gempath = joinpath(path, "bin", gem)
            isfile(gempath) && return gempath
        end
    end

    if Sys.isunix()
        fallback = joinpath("/usr/local/bin", gem)
        isfile(fallback) && return fallback
    end

    return ""
end

using Markdown
"""
    copylocallinks(originalreadme, readmepath)

Copy relative link targets in the `originalreadme` such that they are reachable from `readmepath` as well.
"""
function copylocallinks(originalreadme, readmepath)
    basepath = normpath(joinpath(originalreadme, ".."))
    newbasepath = normpath(joinpath(readmepath, ".."))
    contents = String(read(readmepath))
    md = Markdown.parse(contents)
    links = []
    recurseMDcontents(md, links)
    for link in links
        occursin(r"^http\:\/\/.*", link) && continue
        isabspath(link) && continue
        isfile(joinpath(basepath, link)) || continue

        ispath(dirname(joinpath(newbasepath, link))) || mkpath(dirname(joinpath(newbasepath, link)))
        try
            cp(joinpath(basepath, link), joinpath(newbasepath, link))
        catch err
            @error("copying asset `$(link)` failed", exception=err)
        end
    end
end
function recurseMDcontents(md, links)
    isdefined(md, :content) && foreach(c -> recurseMDcontents(c, links), md.content)
end
recurseMDcontents(md::Markdown.Link, links) = push!(links, md.url)
recurseMDcontents(md::Markdown.Image, links) = push!(links, md.url)

function get_method_from_registry(pspec, registry, root)
    if isfile(registry)
        uuid = string(pspec.uuid)
        toml = Pkg.TOML.parsefile(registry)
        if haskey(toml, uuid)
            pkg = toml[uuid]
            if haskey(pkg, "method") && haskey(pkg, "location")
                return (pkg["method"], pkg["location"])
            else
                @warn("Invalid registry entry for $(pspec.name).")
            end
        else
            @info("$(pspec.name) ($(uuid)) not found in documentation registry.")
        end
    else
        @warn("No registry found. Falling back to `vendored` docs.")
    end
    return ("vendored", joinpath(root, "docs"))
end

function parseall(str)
    pos = firstindex(str)
    exs = []
    while pos <= lastindex(str)
        ex, pos = Meta.parse(str, pos)
        push!(exs, ex)
    end
    if length(exs) == 0
        throw(ParseError("end of input"))
    elseif length(exs) == 1
        return exs[1]
    else
        return Expr(:block, exs...)
    end
end

"""
Rewrites the docs/make.jl for a given package to not deploy and respect
our build folders. This respects any custom scripts needed to build the docs
while disabling deploy + putting the right args into makedocs
"""
function rewrite_makefile(makefile)
    ast = parseall(read(makefile, String))
    make_expr = Expr(:block)
    push!(make_expr.args, :(using Pkg))
    buildpath = joinpath(dirname(makefile), "build")
    for elem in ast.args
        # skip deply(...) - we don't want to deploy
        Meta.isexpr(elem, :call) && elem.args[1] == :deploydocs && continue
        if Meta.isexpr(elem, :call) && elem.args[1] == :makedocs
            # rewrite makedoc call to respect our requirements
            new_args = []
            for argument in elem.args
                if Meta.isexpr(argument, :kw)
                    name, arg = argument.args
                    # assure that we generate HTML
                    if name == :format
                        html = :(Documenter.HTML())
                        if Meta.isexpr(arg, :call) && arg.args[1] == :(Documenter.HTML)
                            append!(html.args, arg.args[2:end])
                        end
                        argument.args[2] = html
                    end
                    # filter out root + build dir
                    if name == :build
                        # if there is a custom build folder, record it!
                        cd(dirname(makefile)) do
                            buildpath = abspath(arg)
                        end
                    end
                end
                push!(new_args, argument)
            end
            push!(new_args,
                Expr(:kw, :root, dirname(makefile)),
            )
            elem = Expr(:call, new_args...)
        end
        # Pkg.add packages required by package build script
        if Meta.isexpr(elem, :using)
            pkgs = Symbol[]
            for arg in elem.args
                if Meta.isexpr(arg, :.) && length(arg.args) == 1
                    push!(pkgs, arg.args[1])
                end
            end
            push!(make_expr.args, Expr(:call, :(Pkg.add), [string.(pkgs)...]))
        end

        push!(make_expr.args, elem)
    end
    return make_expr, buildpath
end

struct PkgNoWork <: Exception
    name::String
end

function install_and_use(pspec)
    # add dependencies as well? will be trivial with a package containing Project.toml
    try
        Pkg.add(pspec)
    catch e
        @warn("Pkg build failed for ", pspec, " due to ", e)
    end
    pkg_sym = Symbol(pspec.name)

    # needs to be main, because documenter wants to have the symbol in main -.-
    pkg_module = try
        @eval(Main, (using $pkg_sym; $pkg_sym))
    catch e
        @warn("Could not load `$pkg_sym`.", exception = e)
        nothing
    end
    pkgdir = Base.find_package(pspec.name)
    if pkgdir == nothing
       throw(PkgNoWork("Package not added"))
    end
    return pkg_module, normpath(joinpath(dirname(Base.find_package(pspec.name)), ".."))
end

function readstr_buffer(x::IOStream)
    return read(x, String)
end

function readstr_buffer(x::Base.GenericIOBuffer{Array{UInt8,1}})
    return String(take!(x))
end

"""
    run_with_timeout(
        command; log=stdout, timeout = 5*60, name = "",
        wait_time = 1, verbose = true
    )

Runs `command` and pipes all output to `log`. The process will be terminated after
`timeout` seconds without any output. `name` describes the process in log messages,
and `verbose` determines whether meta-logs ("process started" etc.) will be printed.
"""
function run_with_timeout(
        command; log=stdout, timeout = 5*60, name = "",
        wait_time = 1, verbose = true, kill_timeout = 60
)

    out_io = IOBuffer()
    err_io = IOBuffer()
    out_file, err_file = "", ""
    if VERSION < v"1.1"
        out_file, out_io = mktemp()
        err_file, err_io = mktemp()
    end
    pipe = pipeline(command, stdout = out_io, stderr = err_io)
    process = run(pipe, wait = false)
    if VERSION < v"1.1"
        out_io = open(out_file)
        err_io = open(err_file)
    end
    timeout_start = time()
    task = @async begin
        logfallback = false
        io = try
            log isa String ? open(log, "w") : log
        catch err
            @error "Error opening logfile, falling back to stdout" error=err
            logfallback = true
            stdout
        end
        try
            tstart = time()
            verbose && @info("starting $name")
            while process_running(process)
                elapsed = (time() - timeout_start)
                if elapsed > timeout
                    verbose && @info("Terminating $name")
                    kill(process)
                    # Handle scenarios where SIGTERM is blocked/ignored/handled by the process
                    start_time = time()
                    while process_running(process)
                        if time() - start_time > kill_timeout
                            verbose && @info("Killing $name")
                            kill(process, signum = SIGKILL)
                        end
                        sleep(5)
                    end
                    break
                end
                errstr, outstr = readstr_buffer.((out_io, err_io))
                is_silent = length(errstr) == 0 && length(outstr) == 0
                isempty(outstr) || println(io, outstr)
                isempty(errstr) || println(io, errstr)

                # if something printed reset timeout
                if !is_silent
                    timeout_start = time()
                end

                sleep(wait_time)
            end

            verbose && @info("$name completed in $(round(time() - tstart, digits=1)) seconds")
        catch err
            @error "Error while running $(name) with timeout." error=err
        finally
            errstr, outstr = readstr_buffer.((out_io, err_io))
            isempty(outstr) || println(io, outstr)
            isempty(errstr) || println(io, errstr)

            flush(io)
            if log isa String && !logfallback
                close(io)
            end
        end
    end
    return process, task
end

"""
    installable_on_version(version = VERSION; registry=joinpath(homedir(), ".julia/registries/General"))

Returns a vector of named tuples `(name, url, versions)` of packages in `registry`
compatible with Julia version `version`.
"""
function installable_on_version(version = VERSION; registry=joinpath(homedir(), ".julia/registries/General"))
    allpkgs = []
    for initial in filter!(isdir, joinpath.(registry, readdir(registry)))
        for pkg in filter!(isdir, joinpath.(registry, initial, readdir(initial)))
            "Compat.toml" in readdir(pkg) || continue
            pkgtoml = Pkg.TOML.parsefile(joinpath(pkg, "Package.toml"))
            versions = Pkg.TOML.parsefile(joinpath(pkg, "Versions.toml"))
            compat = Pkg.TOML.parsefile(joinpath(pkg, "Compat.toml"))
            for pkgver in keys(compat)
                try
                    if haskey(compat[pkgver], "julia")
                        if any(in.(version, Pkg.Types.VersionRange.(compat[pkgver]["julia"])))
                            push!(allpkgs,
                                (
                                    name = pkgtoml["name"],
                                    url = pkgtoml["repo"],
                                    versions = [
                                                VersionNumber(v) for v in keys(versions) if
                                                    VersionNumber(v) in VersionRange(pkgver)
                                               ]
                                )
                            )
                        end
                    end
                catch err
                    @error err
                end
            end
        end
    end
    allpkgs
end

const DOCS_REGISTRY = "https://github.com/JuliaDocs/DocumentationGeneratorRegistry.git"
function download_registry(basepath; sync = true)
    tomlpath = joinpath(basepath, "DocumentationGeneratorRegistry", "Registry.toml")
    if sync
        try
            rm(joinpath(basepath, "DocumentationGeneratorRegistry"), force = true, recursive = true)
            cd(basepath)
            run(`git clone --depth=1 $DOCS_REGISTRY DocumentationGeneratorRegistry`)
            @assert isfile(tomlpath)
            return tomlpath
        catch err
            @warn("Couldn't download docs registry.", exception = err)
        end
    else
        if isfile(tomlpath)
            return tomlpath
        else
            @warn("No registry found at `$(tomlpath)`. Cloning again.")
            return download_registry(basepath; sync = true)
        end
    end
    return ""
end

get_docs_dir(name, uuid) = get_docs_dir(name, UUID(uuid))
get_docs_dir(name, uuid::UUID) = joinpath(name, Base.package_slug(uuid, 5))

"""
    build_documentation(name, url, version; basepath=joinpath(@__DIR__, ".."))

Asynchronously build documentation for package `name`@`version` at `url` and
save the HTML docs to `\$basepath/build` with logs in `\$basepath/logs`.

Note that this will overwrite previous builds/logs.
"""
function build_documentation(uuid, name, url, version;
                             basepath = joinpath(@__DIR__, ".."),
                             envpath = normpath(joinpath(@__DIR__, "..")),
                             juliacmd = first(Base.julia_cmd()),
                             registry_path = "",
                             deployment_url = "")
    workerfile = joinpath(@__DIR__, "worker_work.jl")
    buildpath = joinpath(basepath, "build")
    logpath = joinpath(basepath, "logs")

    isdir(buildpath) || mkpath(buildpath)
    isdir(logpath) || mkpath(logpath)

    builddir = joinpath(buildpath, get_docs_dir(name, uuid), string(version))
    isdir(builddir) || mkpath(builddir)
    logfile = joinpath(logpath, "$(name)-$(uuid) $version.log")
    cmd = `$(juliacmd) --project=$(envpath) --color=no --compiled-modules=no -O0 $workerfile $uuid $name $url $version $builddir $registry_path $deployment_url`
    process, task = run_with_timeout(cmd, log=logfile, name = string("docs build for package ", name, " (", uuid, ")"))
    return process
end

get_uuids(pkg_name::String) = Pkg.Types.registered_uuids(Pkg.Types.Context().env, pkg_name)

function build_documentations(
        packages;
        processes::Int = 8, sleeptime = 0.5,
        juliacmd = first(Base.julia_cmd()),
        basepath = joinpath(@__DIR__, ".."),
        envpath = normpath(joinpath(@__DIR__, "..")),
        filter_versions = last,
        sync_registry = true,
        deployment_url = "pkg.julialang.org/docs"
    )
    regpath = download_registry(basepath, sync = sync_registry)

    process_queue = []
    for package in packages
        uuid = get(package, :uuid, nothing)
        if uuid == nothing
            uuids = get_uuids(package.name)
            length(uuids) > 1 && (@error "Package $(name) has multiple uuids", uuids)
            uuid = uuids[1]
        end
        #those somehow get stuck - might be random
        while length(process_queue) >= processes
            filter!(process_running, process_queue)
            sleep(sleeptime)
        end
        process = nothing
        if !haskey(package, :latest_docs_version)
            for version in vcat(filter_versions(sort(package.versions, rev=true)))
                process = build_documentation(uuid, package.name, package.url, version, envpath = envpath,
                                                    basepath = basepath, juliacmd = juliacmd,
                                                    registry_path = regpath, deployment_url = deployment_url)
                push!(process_queue, process)
            end
        else
            process = build_documentation(uuid, package.name, package.repo, package.latest_docs_version,
                                                envpath = envpath, basepath = basepath, juliacmd = juliacmd,
                                                registry_path = regpath, deployment_url = deployment_url)
            push!(process_queue, process)
        end
    end
    for proc in process_queue
        wait(proc)
    end
end

end


text:
#==============================================================================#
#                                  Upsample                                    #
#==============================================================================#
# length        = desired length of WGN noise vector
# power         = desired power level
# units         = "unitless", "watts", "dBm", "dBW". Case insensitive.
# impedence     = resistance of the load
# returnComplex = true/false. Whether or not to return complex noise
#
# Returned vector is in volts. The function uses paramter 'units' and 'impedence' to determine the RMS voltage level
function wgn(length::Integer; power::Real=1.0, units::String = "unitless", impedence::Real=1.0, returnComplex::Bool=false)
    # TODO: add argument valdiation
    units = lowercase( units )
    if units == "unitless"
        np = power^2/impedence
    elseif units == "watts"
        np = power
    elseif units == "dbm"
        np = 10^((power - 30)/10)
    elseif units == "dbw"
        np = 10^(power/10)
    end

    if returnComplex
        noiseVec = sqrt(impedence*np) * ( randn(length) + randn(length)im )
    else
        noiseVec = sqrt(impedence*np) * randn(length)
    end

    return noiseVec
end


text:
function minimum_degree!(g::AbstractGraph)
    mg = buildmeta(g)
    origin_edges = edges(g)
    for i in 1:nv(g)
        mv = argmindeg(mg)
        sedges = getsaturate(mg, mv)
        add_edges!(mg, sedges)
        sedges::Vector{Vector{Int}} = [edge_m2s(mg, e[1], e[2]) for e in sedges]
        add_edges!(g, sedges)
        rem_vertex!(mg, mv)
    end
    data = Dict("added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function minimum_degree(g::AbstractGraph; kwargs...)
    cg = deepcopy(g)
    data = minimum_degree!(cg)
    return cg, data
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using PowerSystems
cost = VariableCost([(1.0, 1.0), (2.0, 1.1), (3.0, 1.2)])
slopes = get_slopes(cost)
res = [1.0, 10.0, 10.0]
for (ix, v) in enumerate(slopes)
    @test isapprox(v, res[ix])
end

bps = get_breakpoint_upperbounds(cost)
res = [1.0, 0.1, 0.1]
for (ix, v) in enumerate(bps)
    @test isapprox(v, res[ix])
end