text:
import KernelDensityEstimatePlotting: plotKDE

KernelDensityEstimatePlotting.plotKDE(mkd::ManifoldKernelDensity,w...;kw...) = plotKDE(mkd.belief, w...;kw...)
KernelDensityEstimatePlotting.plotKDE(arr::AbstractVector{<:ManifoldKernelDensity},w...;kw...) = plotKDE((x->x.belief).(arr), w...;kw...)

"""
    $(SIGNATURES)

A peneric KDE plotting function that allows marginals of higher dimensional beliefs and various keyword options.

Example:
```julia
p = kde!(randn(3,100))

plotKDE(p)
plotKDE(p, dims=[1;2], levels=3)
plotKDE(p, dims=[1])

q = kde!(5*randn(3,100))
plotKDE([p;q])
plotKDE([p;q], dims=[1;2], levels=3)
plotKDE([p;q], dims=[1])
```
"""
function KernelDensityEstimatePlotting.plotKDE(
                  fgl::AbstractDFG,
                  sym::Symbol;
                  solveKey::Symbol=:default,
                  dims=nothing,
                  title="",
                  levels::Int=5,
                  fill::Bool=false,
                  layers::Bool=false,
                  c=nothing,
                  overlay=nothing  )
  #
  p = getBelief(getVariable(fgl,sym), solveKey)
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  # mp = marginal(p,mmarg)
  bel = p isa BallTreeDensity ? p : p.belief
  plotKDE(bel, levels=levels, dims=dims, title=string(sym, "  ", title), fill=fill, layers=layers, c=c, overlay=overlay )
end
function KernelDensityEstimatePlotting.plotKDE(
                  fgl::AbstractDFG,
                  syms::Vector{Symbol};
                  solveKey::Symbol=:default,
                  addt::Union{<:AbstractVector{<:BallTreeDensity},AbstractVector{<:ManifoldKernelDensity}}=BallTreeDensity[],
                  dims=nothing,
                  title=nothing,
                  levels=3,
                  layers::Bool=false,
                  c=getColorsByLength(length(addt)),
                  overlay=nothing  )
  #
  # TODO -- consider automated rotisary of color
  # colors = ["black";"red";"green";"blue";"cyan";"deepskyblue"; "yellow"]
  # COLORS = repmat(colors, 10)
  # COLORS = getColorsByLength(length(syms))
  MP = BallTreeDensity[]
  LEG = String[]
  # mmarg = Int[]
  for sym in syms
    p = getBelief(getVariable(fgl,sym), solveKey)
    # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
    # mp = marginal(p,mmarg)
    if p isa BallTreeDensity
      push!(MP, p)
    else
      push!(MP, p.belief)
    end
    push!(LEG, string(sym))
  end
  for p in addt
    # mp = marginal(p,mmarg)
    push!(MP, p)
    push!(LEG, "add")
  end
  plotKDE(MP, c=c, levels=levels, dims=dims, legend=LEG, title=title, layers=layers, overlay=overlay)
end

"""
    $(SIGNATURES)

Draw the upward belief from clique `cllb::Symbol` message on variable `lb::Symbol`.

Example:

```julia
plotUpMsgsAtCliq(tree, :x2, :x1)
```

Related

plotKDE, getUpMsgs
"""
function plotUpMsgsAtCliq(treel::AbstractBayesTree,
                          cllb::Symbol,
                          lb::Symbol;
                          show::Bool=true,
                          w=20cm, h=15cm,
                          levels::Int=1,
                          dims::Union{Vector{Int}, Nothing}=nothing )
  #

  cliq = getClique(treel, cllb)
  cliqoutmsg = getUpMsgs(cliq)
  bel = convert(BallTreeDensity, cliqoutmsg.belief[lb])
  plotKDE(bel, dims=dims)
end

function plotTreeUpwardMsgs(fgl::G,
                            bt::AbstractBayesTree;
                            N=300 ) where G <: AbstractDFG
  #
  len = length(bt.cliques)-1
  vv = Array{Gadfly.Compose.Context,1}(undef, len)
  #r = Array{RemoteRef,1}(len)
  i = 0
  for cliq in bt.cliques
      if cliq[1] == 1 println("No upward msgs from root."); continue; end
      @show cliq[2].attributes["label"]
      i+=1
      vv[i] = drawHorDens(fgl, cliq[2].attributes["debug"].outmsg.p, N)
  end
  vv
end

# for some reason we still need msgPlots of right size in the global for these functions to work.
# precall drawTreeUpwardMsgs or drawFrontalDens to make this work properly TODO
function vstackedDensities(msgPlots)
    #msgPlots = f(fg, bt) # drawTreeUpwardMsgs
    evalstr = ""
    for i in 1:length(msgPlots)
        evalstr = string(evalstr, ",msgPlots[$(i)]")
    end
    eval(parse(string("vstack(",evalstr[2:end],")")))
end

"""
    $(SIGNATURES)

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.

DevNotes
- TODO, standardize around ::MIME="image/svg", see JuliaRobotics/DistributedFactorGraphs.jl#640
"""
function plotLocalProduct(fgl::AbstractDFG,
                          lbl::Symbol;
                          solveKey::Symbol=:default,
                          N::Int=100,
                          dims::Vector{Int}=Int[],
                          levels::Int=1,
                          show::Bool=false,
                          dirpath="/tmp/",
                          mimetype::AbstractString="svg",
                          sidelength=20cm,
                          title::String="Local product ($solveKey), ",
                          xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing  )
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getBelief(getVariable(fgl, lbl), solveKey))
  push!(lbls, "curr")
  pl = nothing
  pp, parr, lb, sinfd = IIF.localProduct(fgl, lbl, N=N, solveKey=solveKey)

  partials = []

  # another sanity check
  xmin !== nothing && xmax !== nothing && xmin == xmax ? error("xmin must be less than xmax") : nothing
  ymin !== nothing && ymax !== nothing && ymin == ymax ? error("ymin must be less than ymax") : nothing

  # helper functions
  function plotDirectProducts()
      if pp != parr[1]
        push!(arr,pp)
        push!(lbls, "prod")
        for a in parr
          push!(arr, a)
        end
        @show lb, lbls
        lbls = union(lbls, string.(lb))
      end
      dims = length(dims) > 0 ? dims : collect(1:Ndim(pp))
      colors = getColorsByLength(length(arr))
      plotKDE(arr, dims=dims, levels=levels, c=colors, title=string(title,lbl), legend=string.(lbls)) #
  end
  function plotPartialProducts()
      # stack 1d plots to accomodate all the partials
      PL = []
      lbls = String["prod";"curr";string.(lb)]
      pdims = sort(collect(keys(partials)))
      for dimn in pdims
        vals = partials[dimn]
        proddim = marginal(pp, [dimn])
        colors = getColorsByLength(length(vals)+2)
        newbel_ = getBelief(getVariable(fgl, lbl),solveKey)
        newbel = newbel_ isa ManifoldKernelDensity ? newbel_.belief : newbel_
        pl = plotKDE( [proddim;newbel;vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
        push!(PL, pl)
      end
      Gadfly.vstack(PL...)
  end

  @show length(parr), length(partials)
  if length(parr) > 0 && length(partials) == 0
    pl = plotDirectProducts()
  elseif length(parr) == 0 && length(partials) > 0
    pl = plotPartialProducts()
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # set coordinates accordingly
  co = Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
  pl.coord = co

  # now let's export:
  backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,0.75*sidelength), pl)
  driver = mimetype in ["pdf"] ? "evince" : "eog"
  show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pl
end

"""
    $(SIGNATURES)

Plot---for the cylindrical manifold only---the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.  The linear and circular belief products are returned as a tuple.
"""
function plotLocalProductCylinder(fgl::G,
                                  lbl::Symbol;
                                  N::Int=100,
                                  levels::Int=1,
                                  show=false,
                                  dirpath="/tmp/",
                                  mimetype::AbstractString="svg",
                                  sidelength=30cm,
                                  scale::Float64=0.2 ) where G <: AbstractDFG
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  carr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getKDE(getVariable(fgl, lbl)))
  push!(lbls, "curr")
  pl = nothing
  plc = nothing
  pp, parr, partials, lb = IncrementalInference.localProduct(fgl, lbl, N=N)
  if length(parr) > 0 && length(partials) == 0
    if pp != parr[1]
      push!(arr, marginal(pp,[1]))
      push!(lbls, "prod")
      push!(carr, marginal(pp, [2]))
      for a in parr
        push!(arr, marginal(a,[1]))
        push!(carr, marginal(a,[2]))
      end
      lbls = union(lbls, lb)
    end
    colors = getColorsByLength(length(arr))
    pll = plotKDE(arr, dims=[1;], levels=levels, c=colors, legend=lbls, title=string("Local product, ",lbl))
    plc = AMP.plotKDECircular(carr, c=colors, scale=scale)
    # (pll, plc)
  elseif length(parr) == 0 && length(partials) > 0
    # stack 1d plots to accomodate all the partials
    PL = []
    PLC = []
    lbls = ["prod";"curr";lb]
    pdims = sort(collect(keys(partials)))
    for dimn in pdims
      vals = partials[dimn]
      proddim = marginal(pp, [dimn])
      colors = getColorsByLength(length(vals)+2)
      if dimn == 1
        pll = plotKDE([proddim;getKDE(getVariable(fgl, lbl));vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
        push!(PL, pl)
      else
        plc = AMP.plotKDECircular([proddim; marginal(getKDE(getVariable(fgl, lbl)),[2]);vals], c=colors, scale=scale)
        push!(PLC, plc)
      end
    end
    pl = Gadfly.vstack(PL..., PLC...)
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # now let's export:
  # backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  # Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,sidelength), pl)
  # driver = mimetype in ["pdf"] ? "evince" : "eog"
  # show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pll, plc
end

"""
    $SIGNATURES

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference. String version is obsolete and will be deprecated.
"""
plotLocalProduct(fgl::AbstractDFG, lbl::AbstractString; N::Int=100, dims::Vector{Int}=Int[] ) = plotLocalProduct(fgl, Symbol(lbl), N=N, dims=dims)

"""
    $SIGNATURES

Project (convolve) to and take product of variable in Bayes/Junction tree.

Notes
- assume cliq and var sym the same, unless both specified.
- `cliqsym` defines a frontal variable of a clique.
"""
function plotTreeProductUp(fgl::G,
                           treel::AbstractBayesTree,
                           cliqsym::Symbol,
                           varsym::Symbol=cliqsym;
                           levels::Int=1,
                           dims::Vector{Int}=Int[]  ) where G <: AbstractDFG
  #
  # build a subgraph copy of clique
  cliq = getClique(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraph(fgl, syms)

  # add upward messages to subgraph
  msgs = fetchMsgsUpChildren(treel,cliq, TreeBelief)
  # @show typeof(msgs)
  addMsgFactors!.(subfg, msgs, IIF.UpwardPass)

  # predictBelief
  # stuff = treeProductUp(fgl, treel, cliqsym, varsym)
  # plotKDE(manikde!(stuff[1], getManifolds(fgl, varsym)))
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Up $(cllbl) | ", levels=levels, dims=dims)
end

function plotTreeProductDown( fgl::AbstractDFG,
                              treel::AbstractBayesTree,
                              cliqsym::Symbol,
                              varsym::Symbol=cliqsym;
                              levels::Int=1  )
  #
  # build a subgraph copy of clique
  cliq = whichCliq(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraphFromLabels!(fgl,syms)

  # add upward messages to subgraph
  msgs = getCliqParentMsgDown(treel,cliq)
  addMsgFactors!(subfg, msgs)

  # predictBelief
  # stuff = treeProductUp(fgl, treel, cliqsym, varsym)
  # plotKDE(manikde!(stuff[1], getManifolds(fgl, varsym)))
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Dwn $(cllbl) | ", levels=levels)
end

"""
    $SIGNATURES

Overlay plot all upward messages from cliques.
"""
function plotCliqUpMsgs(fg::G,
                        tree::AbstractBayesTree,
                        sym::Symbol;
                        show::Bool=true,
                        dims::Vector{Int}=Int[],
                        levels::Int=1,
                        c=nothing,
                        title="up msgs on $(sym)"    ) where G <: AbstractDFG
  #
  # get all msgs
  allmsgs = getTreeCliqUpMsgsAll(tree)
  # stack messages by variable
  sckmsgs = stackCliqUpMsgsByVariable(tree, allmsgs)

  if !haskey(sckmsgs, sym)
    @warn "plotCliqUpMsgs -- tree does not have up messages for $sym."
    return nothing
  end

  # prepend current estimate too
  Xs = getBelief(fg, sym)
  # vectorize beliefs
  beliefs = BallTreeDensity[Xs;]
  lbls = String["curr,-1";]
  for msg in sckmsgs[sym]
    push!(beliefs, msg.belief)
    push!(lbls, "$(msg.cliqId),$(msg.depth)")
  end

  # ignoring legend and color information
  cc = c === nothing ? getColorsByLength(length(beliefs)) : c

  # plot and return
  plotKDE(beliefs, levels=levels, c=cc, title=title, legend=lbls)
end

"""
    $SIGNATURES

Development function to plot the same variable from both factor graphs together.
"""
function plotPairVariables(dfg1::G1,
                           dfg2::G2,
                           sym::Symbol;
                           dims=nothing,
                           levels::Int=3,
                           title::String="") where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  X1 = getBelief(dfg1, sym)
  X2 = getBelief(dfg2, sym)

  plotKDE([X1;X2], c=["black";"red"], legend=["dfg1";"dfg2"], dims=dims, levels=levels, title=title*" $sym")
end

function plotPairPose2(dfg1::G1,
                           dfg2::G2,
                           sym::Symbol;
                           dims=nothing,
                           levels::Int=3,
                           title::String="") where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  X1 = getKDE(dfg1, sym)
  X2 = getKDE(dfg2, sym)

  plotPose(Pose2(), [X1;X2], title*" $sym", c=["black";"red"], levels=levels)
  # , legend=["dfg1";"dfg2"]
end

"""
    $SIGNATURES

Plot new proposal (convolution) for factor x of variable y given all other -- i.e. y|X

Notes
- plot new result on `towards` as first color.
- Plot other variables in factor on 2nd to m colors.
"""
function plotVariableGivenFactor(dfg::G,
                                 fct::Symbol,
                                 towards::Symbol;
                                 levels::Int=2,
                                 dims=nothing  ) where G <: AbstractDFG
  #
  pts = approxConv(dfg,fct,towards)
  mani = getManifold(dfg, towards)
  res = manikde!(mani, pts)

  lie = ls(dfg, fct)
  setdiff!(lie, [towards])

  otr = map(x->getKDE(dfg,x),lie)
  lbls = string.([towards;lie])

  pl = plotKDE([res;otr],dims=dims,levels=levels,legend=lbls)

  return pl
end

"""
    $SIGNATURES

Plot the downward messages currently stored in a clique.
"""
function plotCliqDownMsgs(tree::AbstractBayesTree,
                          frnt::Symbol;
                          show::Bool=true,
                          levels::Int=2,
                          dims=nothing,
                          existing=nothing  )
  #
  cliq = getClique(tree,frnt)
  msgs = IIF.getMessageBuffer(cliq).downRx.belief

  PL = []

  for (key, beldim) in msgs
    npl = plotKDE(manikde!(beldim.val, beldim.softtype), levels=levels, title="dwn msg $key", dims=dims)
    existing === nothing ? nothing : union!(npl.layers, existing.layers)
    push!(PL, npl)
  end

  existing === nothing ? nothing : push!(PL, existing)
  pl = vstack(PL...)

  folderpath = "/tmp/caesar/random/"
  filepath = folderpath*"downmsgs_cliq$(cliq.index).pdf"
  Base.mkpath(folderpath)
  pl |> PDF(filepath, 20cm, length(PL)*12cm)
  @async run(`evince $filepath`)
  pl
end


text:
include("linear_system_solvers.jl")
#using advanced_timer
using Main.advanced_timer
#using SparseArrays

start_advanced_timer()
solver_julia = linear_solver_JULIA(:unsymmetric)
initialize!(solver_julia)
ls_factor!(solver_julia, SparseArrays.speye(10), 10, 0)
ls_solve(solver_julia, rand(10))
#mumps_sym(:symmetric)

solver_mumps = linear_solver_MUMPS(:unsymmetric)
initialize!(solver_mumps)
ls_factor!(solver_mumps, SparseArrays.speye(10), 10, 0)
ls_solve(solver_mumps, rand(10))
pause_advanced_timer()

## attempt to produce minimal example that re-produces seg-faults
##

using MUMPS, MPI
MPI.Init()

cntl = MUMPS.default_cntl64[:]; # copy
cntl[1] = 0.01

icntl = MUMPS.default_icntl[:]; # copy

icntl[4] = 1;
icntl[10] = 2; # 2 iterative refinement steps
icntl[14] = 200.0 #1000.0;

n = 1000

my_kkt_solver = Symmetric_KKT_solver()
my_kkt_solver.ls_solver = linear_solver_MUMPS(:symmetric)

my_kkt_solver.ls_solver._factor = MUMPS.Mumps{Float64}(MUMPS.mumps_symmetric, icntl, cntl)

function solves!(my_kkt_solver::abstract_KKT_system_solver)
  MUMPS.associate_rhs!(my_kkt_solver.ls_solver._factor, rand(n));
  MUMPS.solve!(my_kkt_solver.ls_solver._factor);
  sol = MUMPS.get_solution(my_kkt_solver.ls_solver._factor);
  @show LinearAlgebra.norm(sol)
end

function fact!(my_kkt_solver::abstract_KKT_system_solver)
    MUMPS.associate_matrix!(my_kkt_solver.ls_solver._factor, SparseArrays.sparse(rand(n,n)));
    MUMPS.factorize!(my_kkt_solver.ls_solver._factor);
end

for i = 1:50
    fact!(my_kkt_solver)
    for i = 1:50
      solves!(my_kkt_solver)
    end
end

ext: jl


text:
using OmegaGrav
using MatterPower
using Plots, LaTeXStrings
using Dierckx
using DelimitedFiles, Printf

# %% Cosmological parameters
As, ns, kpivot = 2.097e-9, 0.9652, 0.05
Ωm, ΩΛ, Ωb = 0.315, 0.685, 0.049
Ωk = 1 - Ωm - ΩΛ
h0 = 0.674
ωm, fb = Ωm * h0^2, Ωb / Ωm

# %% Loop over redshifts
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
Ωgrav_pklin = zeros(7)
Ωgrav_pknl = zeros(7)
Ωgrav_halo = zeros(7)
for ired = 1:7
   # %% Read in linear and non-linear power spectra computed by CLASS
   filename = @sprintf("data/matterpower05_z%1d_pk.dat", ired)
   d = readdlm(filename, comments = true)
   pk = Spline1D(d[:, 1], d[:, 2])
   filename = @sprintf("data/matterpower05_z%1d_pk_nl.dat", ired)
   d = readdlm(filename, comments = true)
   pknl = Spline1D(d[:, 1], d[:, 2])

   # %% Compute `Ωgrav = Ωm * W / 2`
   z = redshift[ired]
   println("redshift = ", z)
   # Linear matter power spectrum
   Ωgrav_pklin[ired] = ograv_pk(pk, z, Ωm)
   println("Ωgrav(linear pk) = ", Ωgrav_pklin[ired])
   # Non-linear matter power spectrum
   Ωgrav_pknl[ired] = ograv_pk(pknl, z, Ωm)
   println("Ωgrav(non-linear pk) = ", Ωgrav_pknl[ired])
   # Halos
   Ωgrav_halo[ired] = ograv_halo(pk, z, Ωm)
   println("Ωgrav(halos) = ", Ωgrav_halo[ired])
end

# %% Plot results and save to "omegagrav_class.pdf"
p = plot(
   redshift,
   -Ωgrav_pklin,
   xlab = "Redshift",
   ylab = L"-\Omega_{\rm grav}(z)~from~CLASS~P(k)",
   lab = "Linear P(k)",
   yaxis = :log10,
   legend = :bottomleft,
   ls = :dash,
   m = 2,
)
p = plot!(redshift, -Ωgrav_pknl, m = 2, lab = "Non-linear P(k)")
p = plot!(redshift, -Ωgrav_halo, m = 2, lab = "Halos")
savefig("omegagrav_class.pdf")
display(p)


text:
module Testbyrow

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@transform with cols(AsTable)" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    d = @transform df cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ transform(df, [] => (() -> (a = 1, b = 2)) => AsTable)

    d = @transform df cols(AsTable) =  (a = :i, b = :g)

    @test d.a == df.i

    d = @transform df begin
        a = 1
        cols(AsTable) =  (a1 = :i, a2 = :g)
    end

    @test d.a1 == df.i
end

@testset "@select with cols(AsTable) = " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    d = @select df cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ select(df, [] => (() -> (a = 1, b = 2)) => AsTable)

    d = @select df cols(AsTable) =  (a = :i, b = :g)

    @test d.a == df.i

    d = @select df begin
        a = 1
        cols(AsTable) =  (a1 = :i, a2 = :g)
    end

    @test d.a1 == df.i
end

@testset "@transform! with cols(AsTable) = " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    d = @transform! copy(df) cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ transform(df, [] => (() -> (a = 1, b = 2)) => AsTable)

    d = @transform! copy(df) cols(AsTable) =  (a = :i, b = :g)

    @test d.a == df.i

    d = @transform! copy(df) begin
        a = 1
        cols(AsTable) =  (a1 = :i, a2 = :g)
    end

    @test d.a1 == df.i
end

@testset "@select! with cols(AsTable) = " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    d = @select! copy(df) cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ select(df, [] => (() -> (a = 1, b = 2)) => AsTable)

    d = @select! copy(df) cols(AsTable) =  (a = :i, b = :g)

    @test d.a == df.i

    d = @select! copy(df) begin
        a = 1
        cols(AsTable) =  (a1 = :i, a2 = :g)
    end

    @test d.a1 == df.i
end

@testset "@combine with cols(AsTable) = " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    g = groupby(df, :g)

    d = @combine g cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ combine(g, [] => (() -> (a = 1, b = 2)) => AsTable)

    d = @combine df cols(AsTable) =  (a = first(:i), b = first(:g))

    @test d == DataFrame(a = 1, b = 1)

    d = @combine df begin
        a = 1
        cols(AsTable) =  (a1 = 1, a2 = 2)
    end

    @test d == DataFrame(a = 1, a1 = 1, a2 = 2)
end

@testset "@by with cols(AsTable) = " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    g = groupby(df, :g)

    d = @by df :g cols(AsTable) =  (a = 1, b = 2)

    @test d ≅ combine(groupby(df, :g), [] => (() -> (a = 1, b = 2)) => AsTable)
end

end


text:
export Layer

# For a linked-list type layred structure in the future.
mutable struct Layer
    power::SubArray{Float64, 2, Array{Float64,3}, Tuple{Base.Slice{Base.OneTo{Int}}, Base.Slice{Base.OneTo{Int}}, Int}, true}
    phase::SubArray{Float64, 2, Array{Float64,3}, Tuple{Base.Slice{Base.OneTo{Int}}, Base.Slice{Base.OneTo{Int}}, Int}, true}
    core::Array{Float64,3}
    prev::Layer
    next::Layer

    function Layer(sz::Int)
        # t_af, t_fa, r_af, r_fa, pr_f, pr_a
        core = Array{Float64,3}(undef, sz, 6, 2)
        node = new(view(core, :, :, 1), view(core, :, :, 2), core)
        node.prev = node
        node.next = node
        return node
    end
end

function setup!(o::Layer, n::VecI, κ::VecI, fs::VecI, d::Real)
    power = o.power
    phase = o.phase
    #### Forward/Backward Transmission
    transmit!(view(power, :, 1), view(phase, :, 1), n, κ; forward=true)
    transmit!(view(power, :, 2), view(phase, :, 2), n, κ; forward=false)
    #### Forward/Backward Reflection
    reflect!(view(power, :, 3), view(phase, :, 3), n, κ; forward=true)
    reflect!(view(power, :, 4), view(phase, :, 4), n, κ; forward=false)
    #### Propagation in layer/air
    propagate!(view(power, :, 5), view(phase, :, 5), n, κ, fs, d)
    propagate!(view(power, :, 6), view(phase, :, 6), fs, d)
end

function apply!(rs::VecIO, θs::VecIO, o::Layer)
    power = o.power
    phase = o.phase
    one2n = eachindex(rs)
    @simd for i in one2n
        @inbounds rs[i] = power[i,1] * power[i,2] * power[i,5] / power[i,6]
    end
    @simd for i in one2n
        @inbounds θs[i] = phase[i,1] + phase[i,2] + phase[i,5] - phase[i,6]
    end
    @inbounds for i in one2n
        tmpr = power[i,4] * power[i,5]
        tmpθ = phase[i,4] + phase[i,5]
        tmpr, tmpθ = inf_geometric_sum(tmpr * tmpr, 2.0 * tmpθ)
        rs[i] *= tmpr
        θs[i] += tmpθ
    end
end


text:
# When running this file, make sure to set the `--code-coverage=all` command-line flag.

# Important note: even if one or more tests fail, we will still exit with status code 0.

# The reason for this is that we always want to upload code coverage, even if some of the
# tests fail. Therefore, even if the `coverage-linux64` pipeline passes, you should not
# assume that all of the tests passed. If you want to know if all of the tests are passing,
# please look at the status of the `tester_linux64` pipeline.

const include_tests = String[]

const exclude_tests = String[]

empty!(Base.DEPOT_PATH)
push!(Base.DEPOT_PATH, mktempdir(; cleanup = true))

module ChooseTests
    include(joinpath(dirname(dirname(@__DIR__)), "test", "choosetests.jl"))
end

const tests = ChooseTests.choosetests() |>
              first |>
              x -> setdiff(x, exclude_tests) |>
              x -> vcat(x, include_tests) |>
              unique |>
              sort

const ncores = min(Sys.CPU_THREADS, Threads.nthreads())

@info "" ncores Sys.CPU_THREADS Threads.nthreads()

try
    Base.runtests(tests; ncores)
catch ex
    @error "" exception=(ex, catch_backtrace())
end

max_stars_repo_name: chengchingwen/julia

max_issues_repo_name: chengchingwen/julia

max_forks_repo_name: GregPlowman/julia


text:
using ModelingToolkit, OrdinaryDiffEq, Base.Threads

# yeah
include("utils.jl")
param=(nAge=nAge, nHS=nHS, nIso=nIso, nTest=nTest)

#include("diseaseProg.jl")
include("diseaseProg_basic.jl")

# todo replace with inital state tensor
stateTensor= 100*ones((nAge, nHS, nIso, nTest))
dstateTensor=100*ones((nAge, nHS, nIso, nTest))

ssize=(nAge, nHS, nIso, nTest)

nIso=4
nTest=4
nHS=8
nAge=9

ssize=(nAge, nHS, nIso, nTest)
param=(nAge=nAge,nIso=nIso,nTest=nTest)
function f!(dstateTensor,stateTensor_vec,p,t)
  trTensor_complete = zeros(Operation,(nAge, nHS, nIso, nTest, nHS, nIso, nTest))
  trTensor_diseaseProgression=trFunc_diseaseProgression(param...)
  for tst=1:nTest, iso=1:nIso, nhs=1:nHS, ntst=1:nTest, niso=1:nIso, hs=1:nHS, ag=1:nAge
    trTensor_complete[ag,hs,iso,tst,nhs,iso,tst]+=trTensor_diseaseProgression[ag,hs,iso,nhs]
  end
  dstateTensor .= false
  for tst=1:nTest, iso=1:nIso, nhs=1:nHS, ntst=1:nTest, niso=1:nIso, hs=1:nHS, ag=1:nAge
    dstateTensor[ag,hs,iso,tst]+= stateTensor[1,nhs,niso,ntst] * trTensor_complete[ag,nhs,niso,ntst,hs,iso,tst]
  end
end
trTensor_complete = zeros(Operation,(nAge, nHS, nIso, nTest, nHS, nIso, nTest))
trTensor_diseaseProgression=trFunc_diseaseProgression(param...)

f!(dstateTensor,stateTensor,nothing,0.0)
u0= 50*ones(nAge, nHS, nIso, nTest)
prob = ODEProblem(f!,u0,(0.0,80.0))
sys = modelingtoolkitize(prob)
simplified_exs = simplify.(sys.eqs) # Can take a bit, but gives more optimized code
simplifiedsys = ODESystem(simplified_exs,sys.iv,sys.states,sys.ps)
generate_function(simplifiedsys,multithread=true) # Generate the Julia code to save

# sum to get all infacted sum(stateTensor[:,:,8,:],[1,2,3])

#using Pkg
#Pkg.add(PackageSpec(url="https://github.com/JuliaDiffEq/SciPyDiffEq.jl", rev="master"))
# using SciPyDiffEq
# SciPyDiffEq.RK23

#
# fun = lambda t,y: dydt_Complete(t,y, **kwargs),
# t_span=(0.,total_days),
# y0 = cur_stateTensor,
# method='RK23',
# t_eval=range(total_days),
# rtol = 1e-3, #default 1e-3
# atol = 1e-3, # default 1e-6
# SciPyDiffEq code is compiling forever
# sol = solve(prob,SciPyDiffEq.RK23(),rtol=1.0e-3,atol=1.0e-3,saveat=1.0*[0:80;])


text:
"""
    TensorProductMesh(domain::RectangleDomain, n1::Integer, n2::Integer, coordinates::Vector{<:Abstract2DPoint})

A tensor-product `AbstractMesh` of `domain` with `n1` elements in dimension 1, and `n2`
in dimension 2. `coordinates` is an optional argument that can be used to specify the coordinates vector. If `coordinates = nothing`, mesh vertices are defaulted to be equispaced.
"""
struct TensorProductMesh{
    FT,
    CT <: Geometry.Abstract2DPoint{FT},
    RD <: RectangleDomain{CT},
} <: AbstractMesh{FT}
    domain::RD
    n1::Int64 # number of elements in x1 direction
    n2::Int64 # number of elements in x2 direction
    faces::Vector{Tuple{Int64, Int64, Int64, Int64, Bool}}
    coordinates::Vector{CT}
end

function TensorProductMesh(
    domain::RectangleDomain,
    n1,
    n2,
    coordinates = nothing,
)

    nelem = n1 * n2
    x1periodic = isnothing(domain.x1boundary)
    x2periodic = isnothing(domain.x2boundary)
    faces = Vector{Tuple{Int64, Int64, Int64, Int64, Bool}}(undef, nelem * 4)

    # Store all mesh faces as (elem, face, opelem, opface, reversed)
    # so that we can determine face pairs via the map (elem, face) to neighbouring (opelem, opface, reversed)
    for e in 1:nelem
        z2s, z1s = fldmod(e - 1, n1)

        # Face 1
        z1 = z1s - 1
        z2 = z2s
        if z1 < 0 && !x1periodic
            faces[(e - 1) * 4 + 1] = (e, 1, 0, 1, false)
        else
            if z1 < 0
                z1 += n1
            end
            opface = 2
            opelem = z2 * n1 + z1 + 1
            faces[(e - 1) * 4 + 1] = (e, 1, opelem, opface, false)
        end

        # Face 2
        z1 = z1s + 1
        if z1 == n1 && !x1periodic
            faces[(e - 1) * 4 + 2] = (e, 2, 0, 2, false)
        else
            if z1 == n1
                z1 -= n1
            end
            opface = 1
            opelem = z2 * n1 + z1 + 1
            faces[(e - 1) * 4 + 2] = (e, 2, opelem, opface, false)
        end

        # Face 3
        z1 = z1s
        z2 = z2s - 1
        if z2 < 0 && !x2periodic
            faces[(e - 1) * 4 + 3] = (e, 3, 0, 3, false)
        else
            if z2 < 0
                z2 += n2
            end
            opface = 4
            opelem = z2 * n1 + z1 + 1
            faces[(e - 1) * 4 + 3] = (e, 3, opelem, opface, false)
        end

        # Face 4
        z2 = z2s + 1
        if z2 == n2 && !x2periodic
            faces[(e - 1) * 4 + 4] = (e, 4, 0, 4, false)
        else
            if z2 == n2
                z2 -= n2
            end
            opface = 3
            opelem = z2 * n1 + z1 + 1
            faces[(e - 1) * 4 + 4] = (e, 4, opelem, opface, false)
        end
    end

    if isnothing(coordinates)
        CT = Domains.coordinate_type(domain)
        coordinates = Vector{CT}(undef, (n1 + 1) * (n2 + 1))
        x1min = Geometry.component(domain.x1x2min, 1)
        x2min = Geometry.component(domain.x1x2min, 2)
        x1max = Geometry.component(domain.x1x2max, 1)
        x2max = Geometry.component(domain.x1x2max, 2)
        # Default equispaced vertex coordinates, if the user has not specified their locations
        range1 = range(x1min, x1max; length = n1 + 1)
        range2 = range(x2min, x2max; length = n2 + 1)
        # Coordinates array, row-major storage
        for i in 1:(n1 + 1)
            for j in 1:(n2 + 1)
                coordinates[(i - 1) * (n2 + 1) + j] = CT(range1[i], range2[j])
            end
        end
    end
    TensorProductMesh(domain, n1, n2, faces, coordinates)
end

function Base.show(io::IO, mesh::TensorProductMesh)
    print(io, mesh.n1, "×", mesh.n2, " TensorProductMesh of ")
    print(io, mesh.domain)
end


text:
module Metatheory

using RuntimeGeneratedFunctions
using Base.Meta

include("docstrings.jl")

RuntimeGeneratedFunctions.init(@__MODULE__)

# TODO document options
# Options
options = Dict{Symbol, Any}(
    :verbose => false,
    :printiter => false,
)

macro log(args...)
    quote options[:verbose] && @info($(args...)) end |> esc
end

export options

# TODO document this interface
include("expr_abstraction.jl")
export iscall
export istree
export getfunsym
export getfunargs
export setfunsym!
export setfunargs!

include("Util/Util.jl")
using .Util
export Util

include("rgf.jl")
include("rule.jl")
export Rule

include("theory.jl")
include("matchcore_compiler.jl")
include("rewrite.jl")
include("match.jl")

include("EGraphs/EGraphs.jl")
export EGraphs

export @metatheory_init

include("Library/Library.jl")
export Library

export @rule
export @theory

export Theory

export rewrite
export @rewrite
export @esc_rewrite
export @compile_theory
export @matcher
export @rewriter

end # module


text:
using StatsBase: mean
using NQCDynamics: TimeCorrelationFunctions
using .TimeCorrelationFunctions: TimeCorrelationFunction

abstract type AbstractReduction end

"""
Sum the outputs from each trajectory.
"""
struct SumReduction <: AbstractReduction end
(reduction::SumReduction)(u, batch, I) = u + sum(batch), false

"""
Average the outputs over all trajectories.
"""
mutable struct MeanReduction <: AbstractReduction
    counter::Int
end

function (reduction::MeanReduction)(u, batch, I)
    failed_trajectories = length.(batch) .!= length(u)
    deleteat!(batch, failed_trajectories)

    previous = u .* reduction.counter
    reduction.counter += length(batch)
    current = previous .+ sum(batch)

    (current ./ reduction.counter, false)
end

function select_reduction(reduction::Symbol)
    if reduction === :mean
        return MeanReduction(0)
    elseif reduction === :sum
        return SumReduction()
    elseif reduction === :append
        return (u,data,I)->(append!(u,data),false)
    else
        throw(ArgumentError("`reduction` $reduction not recognised."))
    end
end

function get_u_init(::Union{MeanReduction, SumReduction}, saveat, stripped_kwargs, tspan, u0, output::TimeCorrelationFunction)
    if saveat != []
        if saveat isa Number
            savepoints = tspan[1]:saveat:tspan[2]
        else
            savepoints = saveat
        end
    elseif haskey(stripped_kwargs, :dt)
        dt = stripped_kwargs[:dt]
        savepoints = tspan[1]:dt:tspan[2]
    else
        throw(ArgumentError("Make sure to pass either `saveat` or `dt` as keyword arguments
            to compute average/sum over many trajectories."))
    end
    u_init = [TimeCorrelationFunctions.correlation_template(output) for _ ∈ savepoints]
    return u_init
end

function get_u_init(::Union{MeanReduction, SumReduction}, saveat, stripped_kwargs, tspan, u0, output::AbstractOutput)
    return output_template(output, u0)
end

get_u_init(reduction, saveat, stripped_kwargs, tspan, u0, output) = []


text:
#############################################################################
# setdiff.jl
# Handles set differences of combinatorial sets.
#############################################################################

import Base: setdiff

export setdiff, get_elements

mutable struct SetDiffAtom <: CombiSet
  head::Symbol
  id_hash::UInt64
  children::Tuple
  elements::ValOrNothing
  baseset::AbstractArray
  cardinality::Int

  function SetDiffAtom(elements::ValOrNothing, children::Tuple, baseset::AbstractArray)
    this = new(:setdiff, 0, children, elements, baseset, length(baseset))
    this.id_hash = object_id(this)
    return this
  end
end

function setdiff(set1, set2::AbstractArray)
  set01 = Set(set1)
  set02 = Set(set2)
  return collect(setdiff(set01, set02))
end

function setdiff(set1::CombiSet, set2::AbstractArray)
  newset = SetDiffAtom(nothing, (set1, Constant(set2)), set1.baseset)
  return newset
end

function setdiff(set1::AbstractArray, set2::CombiSet)
  newset = SetDiffAtom(nothing, (Constant(set1), set2), set2.baseset)
  return newset
end

function setdiff(set1::CombiSet, set2::CombiSet)
  baseset = collect(union(Set(set1.baseset), Set(set2.baseset)))
  newset = SetDiffAtom(nothing, (set1, set2), baseset)
  return newset
end

function get_elements(x::SetDiffAtom)
  if isa(x.children[1], Constant)
    set1 = Set(x.children[1].value)
    set2 = x.children[2].elements
  elseif isa(x.children[2], Constant)
    set1 = x.children[1].elements
    set2 = Set(x.children[2].value)
  else
    set1 = x.children[1].elements
    set2 = x.children[2].elements
  end
  x.elements = setdiff(set1, set2)
  return x.elements
end


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Typography"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
   showdocs([fontsize, fontweight, fontcolor, fontstyle, fontfamily, fonttype,
             fontcase, textalign, lineheight, letterspacing, title, heading,
             blockquote, caption, emph, code])
)

$(vskip(4em))

"""

function main(window)
    docpage(intro)
end


text:
include("../src/dependencies.jl")

function plot_data_callbacks(time::Float64, state::NetworkState, event::TimedEvent)
    return nothing
end

function proportion_test(params::NetworkParameters; max_time::Float64, lambda_range::Vector{Float64})
    print("Begin proportion test [")
    mean_system_job_totals = zeros(Float64, length(lambda_range))
    proportion_orbiting_jobs = zeros(Float64, length(lambda_range))

    for i in 1:length(lambda_range)
        print("-")
        params.λ = lambda_range[i]
        Random.seed!(1)
        sim_states = simulate(NetworkState(0, zeros(Int8, params.L), 0), TimedEvent(ArrivalEvent(),0.0, 0, 0.0), params, max_time = max_time, callback = plot_data_callbacks)

        # Test Criteria
        proportion_data = sim_states.orbiting_jobs/sim_states.number_in_system
        !isnan(proportion_data) ? proportion_orbiting_jobs[i] = proportion_data : proportion_orbiting_jobs[i] = 0.0
        @assert proportion_orbiting_jobs[i] >= 0.0 && proportion_orbiting_jobs[i] <= 1.0

    end
    print("]")
    println(proportion_orbiting_jobs)
    println("proportion_test - PASSED!")
end

"""
Call simulate function below for lambdas in lambda_range
"""

test_params = NetworkParameters(  L=3,
                                gamma_shape = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(3),
                                P = [0 1.0 0;
                                    0 0 1.0;
                                    0.5 0 0],
                                Q = [0 0.5 0;
                                     0 0 0.5;
                                     0.5 0 0],
                                p_e = AnalyticWeights([1.0, 0, 0]),
                                K = fill(5,3))

lambdas = [i for i in 0.2:0.2:1.2]

proportion_test(test_params, max_time = 1000.0, lambda_range = lambdas)


text:
# Copyright: Wen-Chih Chen
# The example detects outliers based on Chen and Johnson (2010)
# Chen, W. C., & Johnson, A. L. (2010). A unified model for detecting efficient and inefficient outliers
# in data envelopment analysis. Computers & Operations Research, 37(2), 417-425.
# it will write analysis result to output.csv

using JuMP
using Gurobi # Gurobi is used as the LP solver
# uncomment the following for a large-scale problem or with size-limited solvers
# using FrontierEfficiencyAnalysis

function outlierDetection()
    # ---------------------------------
    # data section
    # ---------------------------------
    # read data from a csv file
    println(">> read data")
    data = readcsv("outlierExample.csv") # input User's (.csv) data path
    # each row is for a DMU; each column is an input or output

    # you need to set the inputNo correspondingly
    inputNo = 1 # the first two columns are the inputs while the others are the outputs
    # determine the size of the data
    scale, dimension = size(data) # scale is the number of DMU, dimension is the total number of inputs and outputs
    inputs = data[:, 1:inputNo] # input data
    outputs = data[:, inputNo+1:dimension] # output data
    # end of data section
    # we have input data in inputs; output data in outputs
    # ---------------------------------

    println(">> solving ...")
    # output-oriented
    # outerImpact: vector for the total delta_o for removing the corresponding DMU
    # innerImpact: vector for the total delta_i for removing the corresponding DMU
    # outerImpactNo: vector for the no. of DMU impacted due to removing the corresponding DMU
    # innerImpactNo: vector for the no. of DMU impacted due to removing the corresponding DMU
    outerImpact, innerImpact, outerImpactNo, innerImpactNo = getDeltaOnY(inputs, outputs)
    totalImpact = outerImpact + innerImpact # eq. (8)

    # # for input-oriented models
    # outerImpact, innerImpact, outerImpactNo, innerImpactNo = getDeltaOnX(inputs, outputs)
    # totalImpact = outerImpact + innerImpact # eq. (15)

    # output section: export reslts to output.csv
    writeResult("output.csv", [collect(1:scale)';outerImpact'; outerImpactNo'; innerImpact'; innerImpactNo']')

    println("done!")

end

# detect outlier measures in output-oriented models
# return (6) and (7) in Chen and Johnson (2010) in absolute value
function getDeltaOnY(inputs, outputs)
    scale = size(inputs)[1]
    outerDistWR = zeros(scale)
    innerDistWR = zeros(scale)
    for k = 1:scale
        outerDistWR[k] = outerDistOnY(inputs, outputs, inputs[k,:], outputs[k,:])
        innerDistWR[k] = innerDistOnY(inputs, outputs, inputs[k,:], outputs[k,:])
    end

    outerDistDelta = zeros(scale)
    innerDistDelta = zeros(scale)
    outerImpactNo = zeros(scale)
    innerImpactNo = zeros(scale)
    for delPt = 1:scale # DMU to be reomoved (R)
        set = collect(1:scale)
        deleteat!(set, delPt)

        # determine the total impact of removing delPt for each DMU except of delPt itself
        for k = 1:scale
            if k != delPt
                outerDiff = abs(outerDistWR[k]-outerDistOnY(inputs[set, :], outputs[set, :], inputs[k, :], outputs[k, :]))
                outerDistDelta[delPt] = outerDistDelta[delPt] + outerDiff
                if outerDiff > 0
                    outerImpactNo[delPt] = outerImpactNo[delPt] + 1
                end
                innerDiff = abs(innerDistWR[k]-innerDistOnY(inputs[set, :], outputs[set, :], inputs[k, :], outputs[k, :]))
                innerDistDelta[delPt] = innerDistDelta[delPt] + innerDiff
                if innerDiff > 0
                    innerImpactNo[delPt] = outerImpactNo[delPt] + 1
                end
            end
        end
    end
    return outerDistDelta, innerDistDelta, outerImpactNo, innerImpactNo
end

# detect outlier measures in input-oriented models
# return (13) and (14) in Chen and Johnson (2010) in absolute value
function getDeltaOnX(inputs, outputs)
    scale = size(inputs)[1]
    outerDistWR = zeros(scale)
    innerDistWR = zeros(scale)
    for k = 1:scale
        outerDistWR[k] = outerDistOnX(inputs, outputs, inputs[k,:], outputs[k,:])
        innerDistWR[k] = innerDistOnX(inputs, outputs, inputs[k,:], outputs[k,:])
    end

    outerDistDelta = zeros(scale)
    innerDistDelta = zeros(scale)
    outerImpactNo = zeros(scale)
    innerImpactNo = zeros(scale)
    for delPt = 1:scale # DMU to be reomoved (R)
        set = collect(1:scale)
        deleteat!(set, delPt)

        # determine the total impact of removing delPt for each DMU except of delPt itself
        for k = 1:scale
            if k != delPt
                outerDiff = abs(outerDistWR[k]-outerDistOnX(inputs[set, :], outputs[set, :], inputs[k, :], outputs[k, :]))
                outerDistDelta[delPt] = outerDistDelta[delPt] + outerDiff
                if outerDiff > 0
                    outerImpactNo[delPt] = outerImpactNo[delPt] + 1
                end
                innerDiff = abs(innerDistWR[k]-innerDistOnX(inputs[set, :], outputs[set, :], inputs[k, :], outputs[k, :]))
                innerDistDelta[delPt] = innerDistDelta[delPt] + innerDiff
                if innerDiff > 0
                    innerImpactNo[delPt] = outerImpactNo[delPt] + 1
                end
            end
        end
    end
    return outerDistDelta, innerDistDelta, outerImpactNo, innerImpactNo
end

# return the distance to the output boundary wrt to DMU t in an output-oriented analysis
# eq. (1) in Chen and Johnson (2010)
# inputs: input data; outputs: output data; t: DMU t under evaluation
function outerDistOnY(inputs, outputs, inputRHS, outputRHS)
    scale, inputNo = size(inputs)
    outputNo = size(outputs)[2]

    # modeling begins
    model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
    @variable(model, Lambda[1:scale] >= 0)
    @variable(model, Theta)
    @objective(model, Max, Theta)

    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) == inputRHS[i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) == Theta*outputRHS[j])
    @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)

    # use solveDEA(model) instead of solve(model) in the following for a large-scale problem or with size-limited solvers
    solve(model)
    # solveDEA(model)

    return getvalue(Theta)
end

# return the distance to the inner boundary wrt to DMU t in an output-oriented analysis
# eq. (2) in Chen and Johnson (2010)
# inputs: input data; outputs: output data; t: DMU t under evaluation
function innerDistOnY(inputs, outputs, inputRHS, outputRHS)
    scale, inputNo = size(inputs)
    outputNo = size(outputs)[2]
    # modeling begins
    model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
    @variable(model, Lambda[1:scale] >= 0)
    @variable(model, Theta)
    @objective(model, Min, Theta)
    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) == inputRHS[i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) == Theta*outputRHS[j])
    @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)

    # use solveDEA(model) instead of solve(model) in the following for a large-scale problem or with size-limited solvers
    solve(model)
    # solveDEA(model)

  return getvalue(Theta)
end

# return the distance to the output boundary wrt to DMU t in an output-oriented analysis
# eqs. (11) and (12) in Chen and Johnson (2010)
# inputs: input data; outputs: output data; t: DMU t under evaluation
function outerDistOnX(inputs, outputs, inputRHS, outputRHS)
    scale, inputNo = size(inputs)
    outputNo = size(outputs)[2]
    # modeling begins
    model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
    @variable(model, Lambda[1:scale] >= 0)
    @variable(model, Theta)
    @objective(model, Max, Theta)
    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) == Theta*inputRHS[i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) == outputRHS[j])
    @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)

    # use solveDEA(model) instead of solve(model) in the following for a large-scale problem or with size-limited solvers
    solve(model)
    # solveDEA(model)

    return getvalue(Theta)
end

# return the distance to the output boundary wrt to DMU t in an output-oriented analysis
# eqs. (9) and (10) in Chen and Johnson (2010)
# inputs: input data; outputs: output data; t: DMU t under evaluation
function innerDistOnX(inputs, outputs, inputRHS, outputRHS)
    scale, inputNo = size(inputs)
    outputNo = size(outputs)[2]
    # modeling begins
    model = Model(solver = GurobiSolver(OutputFlag=0)) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
    @variable(model, Lambda[1:scale] >= 0)
    @variable(model, Theta)
    @objective(model, Min, Theta)
    @constraint(model, inputCon[i=1:inputNo], sum(Lambda[r]*inputs[r,i] for r = 1:scale) == Theta*inputRHS[i])
    @constraint(model, outputCon[j=1:outputNo], sum(Lambda[r]*outputs[r,j] for r = 1:scale) == outputRHS[j])
    @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)

    # use solveDEA(model) instead of solve(model) in the following for a large-scale problem or with size-limited solvers
    solve(model)
    # solveDEA(model)

    return getvalue(Theta)
end

# write result to file
function writeResult(file, result)
    # write
    f = open(file, "w")
      write(f, "DMU, tol delta_O, no of DMUs impacted (O), tol delta_I, no of DMUs impacted (I)\r\n")
    close(f)

    for r = 1 : size(result)[1]
        f = open(file, "a")

        for c = 1 : size(result)[2]
            write(f, "$(result[r,c]),")
        end
        write(f, "\r\n")
        close(f)
    end
end

outlierDetection()


text:
using CSV
using DataFrames
using DelimitedFiles

include("constants.jl")

#path_dat = "C://Users//rapiduser//Documents//GitHub//grace-arpae//UnitCommitment_BAU"
#cd(path_dat)

##
# Importing input data from the input spreadsheets
# Generators' specifications
DF_Generators = CSV.read(".//inputs//data_generators.csv", DataFrame);
# Generators location map: if a generator g is located in zone z Map_Gens[g,z]=1; and 0 otherwise
Map_Gens = readdlm(".//inputs//location_generators.csv", ','; header = true);

#Peaker Units' specifications and location
DF_Peakers = CSV.read(".//inputs//data_peakers.csv", DataFrame);
# Peakers location map: if a peaker p is located in zone z Map_Gens[p,z]=1; # and 0 otherwise
Map_Peakers = readdlm(".//inputs//location_peakers.csv", ','; header = true);

# Storage Units' specification and location
DF_Storage = CSV.read(".//inputs//data_storage.csv", DataFrame);
Map_Storage = readdlm(".//inputs//location_storage.csv", ','; header = true);

# Energy demand at each location
FUCR_Demands = readdlm(".//inputs//data_demand.csv", ','; header = true);
SUCR_Demands = readdlm(".//inputs//data_demand_updated.csv", ','; header = true);
BUCR_Demands = readdlm(".//inputs//data_demand_actual.csv", ','; header = true);

# solar generation data at each location
FUCR_SolarGs = readdlm(".//inputs//data_solar.csv", ','; header = true);
SUCR_SolarGs = readdlm(".//inputs//data_solar_updated.csv", ','; header = true);
BUCR_SolarGs = readdlm(".//inputs//data_solar_actual.csv", ','; header = true);

# wind energy data for each location
FUCR_WindGs = readdlm(".//inputs//data_wind.csv", ','; header = true);
SUCR_WindGs = readdlm(".//inputs//data_wind_updated.csv", ','; header = true);
BUCR_WindGs = readdlm(".//inputs//data_wind_actual.csv", ','; header = true);

#hydro generation data for each location
FUCR_HydroGs = readdlm(".//inputs//data_hydro.csv", ','; header = true);
SUCR_HydroGs = readdlm(".//inputs//data_hydro_updated.csv", ','; header = true);
BUCR_HydroGs = readdlm(".//inputs//data_hydro_actual.csv", ','; header = true);

#nuclear generation timeseries for each location
FUCR_NuclearGs = readdlm(".//inputs//data_nuclear.csv", ','; header = true);
SUCR_NuclearGs = readdlm(".//inputs//data_nuclear_updated.csv", ','; header = true);
BUCR_NuclearGs = readdlm(".//inputs//data_nuclear_actual.csv", ','; header = true);

#Cogenerators' generation timeseries for each location
FUCR_CogenGs = readdlm(".//inputs//data_cogen.csv", ','; header = true);
SUCR_CogenGs = readdlm(".//inputs//data_cogen_updated.csv", ','; header = true);
BUCR_CogenGs = readdlm(".//inputs//data_cogen_actual.csv", ','; header = true);

TranC = readdlm(".//inputs//LineCapacity.csv", ','; header = true);
TranS = readdlm(".//inputs//LineSusceptance.csv", ','; header = true);
Reserve_Reqs = readdlm(".//inputs//data_reserve_reqs.csv", ','; header = true);

FuelPrice = readdlm(".//inputs//data_fuel_price.csv", ','; header = true);
FuelPricePeakers = readdlm(".//inputs//data_fuel_price_peakers.csv", ','; header = true);

FuelPrice_head = FuelPrice[2];
FuelPrice = FuelPrice[1];
FuelPrice = Array{Float64}(FuelPrice[2:GENS+1, 4:368]);

FuelPricePeakers_head = FuelPricePeakers[2];
FuelPricePeakers = FuelPricePeakers[1];
FuelPricePeakers = Array{Float64}(FuelPricePeakers[2:PEAKERS+1, 4:368]);

# Reorganize data
FUCR_Demands_head = FUCR_Demands[2];
FUCR_Demands = FUCR_Demands[1];
FUCR_Demands = Array{Int64}(FUCR_Demands[1:60313, 4:(4+N_ZONES-1)]);

SUCR_Demands_head = SUCR_Demands[2];
SUCR_Demands = SUCR_Demands[1];
SUCR_Demands = Array{Int64}(SUCR_Demands[1:60313, 4:(4+N_ZONES-1)]);

BUCR_Demands_head = BUCR_Demands[2];
BUCR_Demands = BUCR_Demands[1];
BUCR_Demands = BUCR_Demands[1:8760, 3:(3+N_ZONES-1)];
BUCR_Demands = Array{Int64}(BUCR_Demands);

# Solar
FUCR_SolarGs_head = FUCR_SolarGs[2];
FUCR_SolarGs = FUCR_SolarGs[1];
FUCR_SolarGs = Array{Int64}(FUCR_SolarGs[1:8760, 3:(3+N_ZONES-1)]);

SUCR_SolarGs_head = SUCR_SolarGs[2];
SUCR_SolarGs = SUCR_SolarGs[1];
SUCR_SolarGs = SUCR_SolarGs[1:8760, 3:(3+N_ZONES-1)];
SUCR_SolarGs = Array{Int64}(SUCR_SolarGs);

BUCR_SolarGs_head = BUCR_SolarGs[2];
BUCR_SolarGs = BUCR_SolarGs[1];
BUCR_SolarGs = Array{Int64}(BUCR_SolarGs[1:8760, 3:(3+N_ZONES-1)]);

#Wind
FUCR_WindGs_head = FUCR_WindGs[2];
FUCR_WindGs = FUCR_WindGs[1];
FUCR_WindGs = Array{Int64}(FUCR_WindGs[1:8760, 3:(3+N_ZONES-1)]);

SUCR_WindGs_head = SUCR_WindGs[2];
SUCR_WindGs = SUCR_WindGs[1];
SUCR_WindGs = Array{Int64}(SUCR_WindGs[1:8760, 3:(3+N_ZONES-1)]);

BUCR_WindGs_head = BUCR_WindGs[2];
BUCR_WindGs = BUCR_WindGs[1];
BUCR_WindGs = Array{Int64}(BUCR_WindGs[1:8760, 3:(3+N_ZONES-1)]);

#Hydro
FUCR_HydroGs_head = FUCR_HydroGs[2];
FUCR_HydroGs = FUCR_HydroGs[1];
FUCR_HydroGs = Array{Int64}(FUCR_HydroGs[1:8760, 2:(2+N_ZONES-1)]);

SUCR_HydroGs_head = SUCR_HydroGs[2];
SUCR_HydroGs = SUCR_HydroGs[1];
SUCR_HydroGs = Array{Int64}(SUCR_HydroGs[1:8760, 2:(2+N_ZONES-1)]);

BUCR_HydroGs_head = BUCR_HydroGs[2];
BUCR_HydroGs = BUCR_HydroGs[1];
BUCR_HydroGs = Array{Int64}(BUCR_HydroGs[1:8760, 2:(2+N_ZONES-1)]);

#Nuclear
FUCR_NuclearGs_head = FUCR_NuclearGs[2];
FUCR_NuclearGs = FUCR_NuclearGs[1];
FUCR_NuclearGs = Array{Int64}(FUCR_NuclearGs[1:8760, 3:(3+N_ZONES-1)]);

SUCR_NuclearGs_head = SUCR_NuclearGs[2];
SUCR_NuclearGs = SUCR_NuclearGs[1];
SUCR_NuclearGs = Array{Int64}(SUCR_NuclearGs[1:8760, 3:(3+N_ZONES-1)]);

BUCR_NuclearGs_head = BUCR_NuclearGs[2];
BUCR_NuclearGs = BUCR_NuclearGs[1];
BUCR_NuclearGs = Array{Int64}(BUCR_NuclearGs[1:8760, 3:(3+N_ZONES-1)]);

#Cogen
FUCR_CogenGs_head = FUCR_CogenGs[2];
FUCR_CogenGs = FUCR_CogenGs[1];
FUCR_CogenGs = Array{Int64}(FUCR_CogenGs[1:8760, 3:(3+N_ZONES-1)]);

SUCR_CogenGs_head = SUCR_CogenGs[2];
SUCR_CogenGs = SUCR_CogenGs[1];
SUCR_CogenGs = Array{Int64}(SUCR_CogenGs[1:8760, 3:(3+N_ZONES-1)]);

BUCR_CogenGs_head = BUCR_CogenGs[2];
BUCR_CogenGs = BUCR_CogenGs[1];
BUCR_CogenGs = Array{Int64}(BUCR_CogenGs[1:8760, 3:(3+N_ZONES-1)]);

Map_Gens_head = Map_Gens[2];
Map_Gens = Map_Gens[1];
Map_Gens = Array{Int64}(Map_Gens[:,3:N_ZONES+2]);

Map_Peakers_head = Map_Peakers[2];
Map_Peakers = Map_Peakers[1];
Map_Peakers = Array{Int64}(Map_Peakers[:,3:N_ZONES+2]);

Map_Storage_head = Map_Storage[2];
Map_Storage = Map_Storage[1];
Map_Storage = Array{Int64}(Map_Storage[:,2:N_ZONES+1])

TranC_head = TranC[2];
TranC = TranC[1];
TranC = Array{Float64}(TranC[1:N_ZONES,2:(2+N_ZONES-1)]);

TranS_head = TranS[2];
TranS = TranS[1];
TranS = Array{Float64}(TranS[1:N_ZONES,2:(2+N_ZONES-1)]);

Reserve_Reqs_head = Reserve_Reqs[2];
Reserve_Reqs = Reserve_Reqs[1];
Reserve_Req_Up = Array{Float64}(Reserve_Reqs)

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
# test for mat_center function
@testset "Test mat_center function $mat" for mat in ["her", "rec"]
    # test for Hermitian matrix
    if mat == "her"
        for idx = 1: 10
            @test abs(tr(mat_center(randn(3,3); mat = mat))) <= 10*eps(Float64);
        end
    end

    # test for rectangular case
    if mat == "rec"
        for idx = 1: 10
            @test sum(vec(mat_center(randn(3,3); mat ="rec"))) <= 10*eps(Float64)
        end
    end
end

# test for free_whiten function
@testset "Test free_whiten function $mat" for mat in ["her", "rec"]
    # test for Hermitian matrix
    if mat == "her"
        for idx = 1: 10
            G = randn(10, 10)
            Z = [G + G', G*G']
            Y = free_whiten(Z; mat = "her")[1]
            @test maximum(abs.( tr.(Y*Y')./10 - I)) < 10*eps(Float64)
        end
    end

    # test for rectangular matrix
    if mat == "rec"
        for idx = 1: 10
            Z = [randn(10, 20) for i = 1: 3]
            Y = free_whiten(Z; mat = "rec")[1]
            @test maximum(abs.( tr.(Y*Y')./10 - I)) < 10*eps(Float64)
        end
    end
end


text:
export randomw
"""
function randomw(n,nsamp=1)

Samples from a uniform distribution in the set

{w in R^n |  sum(w)=1, diff(w).<=0, w.>=0}

Input:

    1. n        - number of strains

Optional argument:

    2. nsamp    - number of samples (default=1)

Output:

    1. wsamples - a matrix of size (n,nsamp) or a vector of length n if nsamp==1
"""
function randomw(n::Int64,nsamp::Int64=1)

    if n == 1
        return nsamp == 1 ? ones(1) : ones(1,nsamp)
    end

    wsamples = zeros(n, nsamp)

    for k = 1:nsamp
        w = [0; sort!(rand(n-1)); 1]
        w = sort(diff(w), rev=true)
        wsamples[:, k] = w;
    end

    if nsamp==1
        return vec(wsamples)
    else
        return wsamples
    end

end


text:
module Objective
export cost_minimization

using JuMP

function cost_minimization(model_dict)
    """
    Minimizes total system monetary cost.
    """
    backend_model= model_dict["backend_model"]
    cost = model_dict["variables"]["cost"]
    loc_techs_cost = model_dict["sets"]["loc_techs_cost"]
    @objective(backend_model, :Min,
        sum(cost["monetary", loc_tech] for loc_tech in loc_techs_cost)
    )
end
end

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


text:
import Base: append!, eye, size

export ConvergenceHistory, KrylovSubspace

*(f::Function, v::AbstractVector)=f(v) #Syntax to mimic matvec product

type Eigenpair{S,T}
    val::S
    vec::Vector{T}
end

type KrylovSubspace{T, OpT}
    A::OpT     #The linear operator that generates the subspace
    n::Int     #Dimension of problem
    order::Int #Order (maximum size) of subspace
    v::Vector{Vector{T}} #The Krylov vectors
    mvps::Int  #Count of matrix-vector products
end

KrylovSubspace(A, n::Int, order::Int, T::Type)=KrylovSubspace{T,typeof(A)}(A, n, order, Vector{T}[], 0)

KrylovSubspace{T}(A, n::Int, order::Int, v::Vector{Vector{T}}=Vector{T}[])=KrylovSubspace{T,typeof(A)}(A, n, order, v, 0)

KrylovSubspace{T}(A::AbstractMatrix{T}, order::Int, v::Vector{Vector{T}}=Vector{T}[])=
    KrylovSubspace{T,typeof(A)}(A, size(A,2), order, v, 0)

#Reset an existing KrylovSubspace
function KrylovSubspace{T}(A::KrylovSubspace{T}, order::Int=size(A,2), v::Vector{Vector{T}}=Vector{T}[])
    A.order = order
    A.v = v
    A.mvps = 0
end

lastvec(K::KrylovSubspace) = K.v[end]
function nextvec(K::KrylovSubspace)
    K.mvps += 1
    K.A*lastvec(K)
end

size(K::KrylovSubspace) = length(K.v)
function size(K::KrylovSubspace, n::Int)
    if isa(K.A, AbstractMatrix)
        return size(K.A, n)
    else
        return n==2 ? size(K) : error("Illegal dimension $n of KrylovSubspace")
    end
end
eye{T}(K::KrylovSubspace{T}) = isa(K.A, Function) ? x->x : eye(T, size(K.A)...)

function append!{T}(K::KrylovSubspace{T}, w::Vector{T})
    if size(K) == K.order; shift!(K.v) end
    push!(K.v, w)
end
appendunit!{T}(K::KrylovSubspace{T}, w::Vector{T}) = append!(K, w/norm(w))

#Initialize the KrylovSubspace K with a random unit vector
function initrand!{T}(K::KrylovSubspace{T})
    K.v = Vector{T}[initrand!(Array(T, K.n))]
end

#Initialize the KrylovSubspace K with a specified nonunit vector
#If nonzero, try to normalize it
function init!{T}(K::KrylovSubspace{T}, v::Vector{T})
#    K.v = Vector{T}[all(v.==zero(T)) ? v : v/norm(v)]
    K.v = Vector{T}[copy(v)]
end

#Orthogonalize a vector v against the last p basis vectors defined by the
#KrylovSubspace K
function orthogonalize{T}(v::Vector{T}, K::KrylovSubspace{T}, p::Int=K.order;
    method::Symbol=:ModifiedGramSchmidt, normalize::Bool=false)
    Kk = K.v[max(1,end-p+1):end]
    if method == :GramSchmidt
        cs = T[dot(e, v) for e in Kk]
        for (i, e) in enumerate(Kk)
            v -= cs[i] * e
        end
    elseif method == :ModifiedGramSchmidt || method== :Householder
        #The numerical equivalence of ModifiedGramSchmidt and Householder was
        #established in doi:10.1137/0613015
        cs = zeros(T, p)
        for (i, e) in enumerate(Kk)
            cs[i] = dot(e, v)
            v-= cs[i] * e
        end
    else
        error("Unsupported orthogonalization method: $(method)")
    end
    normalize && (v /= convert(T, norm(v)))
    v, cs #Return orthogonalized vector and its coefficients
end

#Return a vector in the Krylov basis, i.e. compute a linear combination of Krylov vectors
#Index 1 = oldest remembered vector
function *(a::AbstractVector, K::KrylovSubspace)
    N = length(a)
    N==length(K.v) && throw(DimensionMismatch("Tried to compute linear combination of $N Krylov vectors but only $(length(K.v)) are present."))
    w = a[N] * K.v[N]
    for j = (N-1):-1:1
        w += a[j] * K.v[j]
    end
    w
end


text:
using SchattenNorms
using Base.Test

@test isapprox(snorm([1 0; 0 -1]),sqrt(2.0))
@test isapprox(snorm([1 0; 0 -1],1.0),2.0)
@test isapprox(snorm([1 0; 0 -1],Inf),1.0)

for i=1:100
    for d=2:4
        R = randn(d,d)
        U,_,_ = svd(randn(d,d))
        V,_,_ = svd(randn(d,d))
        p = sort(abs(randn(3))+1.0)
        @test isapprox(snorm(R),snorm(R,2.0))
        @test isapprox(snorm(U*R*V,1),snorm(R,1))
        @test isapprox(snorm(U*R*V),snorm(R,2.0))
        @test isapprox(snorm(U*R*V,Inf),snorm(R,Inf))
        @test isapprox(snorm(U*R*V,p[1]),snorm(R,p[1]))
        @test snorm(R,1) >= snorm(R,2) >= snorm(R,Inf)
        @test snorm(R,p[1]) >= snorm(R,p[2]) >= snorm(R,p[3])
        @test snorm(R,1.0) == snorm(R,1) == trnorm(R) == nucnorm(R)
    end
end

# E
# U,V
# @test isapprox(cbnorm(E),dnorm(E'))
# @test isapprox(cbnorm(E),cbnorm(U*E*V))
# @test isapprox(dnorm(E),dnorm(U*E*V))

include("test-dnorm.jl")


text:
using Test, BenchmarkTools
using LinearAlgebra
using SparseArrays
using WriteVTK
using WaveProp
using WaveProp.Mesh
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Interpolation
using WaveProp.FEM
using WaveProp.IO

vtx = ((-1,0),(1,0),(0,1))
el  = LagrangeTriangle(vtx)
u   = LagrangeBasis{ReferenceTriangle,1}()
v   = LagrangeBasis{ReferenceTriangle,1}()
q   = Gauss{ReferenceTriangle,3}()

M   = zeros(3,3)
FEM.elementary_matrix!(M,el,u,v,q)
M2  = FEM.mass_matrix(el,u,v,q)
@test M ≈ M2
M3 = zeros(3,3)
FEM.mass_matrix_unrolled!(M3,el,u,v,q)
@test M3 ≈ M

@btime FEM.elementary_matrix!($M,$el,$u,$v,$q)
@btime FEM.mass_matrix($el,$u,$v,$q)
@btime FEM.mass_matrix_unrolled!($M3,$el,$u,$v,$q)

M   = zeros(3,3)
ak(u, v) = (i,j,x̂,el,x) -> (inv(transpose(jacobian(el,x̂))) * grad(u)(x̂)[j]) ⋅ (inv(transpose(jacobian(el,x̂))) *  grad(v)(x̂)[i])
FEM.elementary_matrix!(M,el,u,v,q;f=ak)
M2 = zeros(3,3)
FEM.stiffness_matrix_unrolled!(M2,el,u,v,q)
@test M2 ≈ M

@btime FEM.elementary_matrix!($M,$el,$u,$v,$q;f=$ak)
@btime FEM.stiffness_matrix_unrolled!($M2,$el,$u,$v,$q)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# Starting sorting example for Learning Julia

arr1 = [2, 8, -3, -15, 5, -7, 9, 0, 11]

# sort the data by a simple transformation
result1 = sort(arr1, by = abs) # sorted by the absolute values (i.e. ignore the sign)
println(result1, "\n")

# check to see if an array is already sorted
arrtup = [(1, 2), (2, 1), (3, 5), (4, 0)]
println("Is arrtup sorted by 1st element: ", issorted(arrtup, by = x -> x[1]))
# the by argument specifies to look for only the first item in a tuple for sorted order
println("Is arrtup sorted by 2nd element: ", issorted(arrtup, by = x -> x[2]), "\n")

# Sorting custom types
struct MyRect
    length::Int
    width::Int
end

arr2 = [MyRect(12, 22), MyRect(10, 18), MyRect(7, 12), MyRect(9, 30)]

# declare a custom sorting function
function compareArea(rect1::MyRect, rect2::MyRect)
    area1 = rect1.length * rect1.width
    area2 = rect2.length * rect2.width
    return area1 < area2
end

sort!(arr2, lt = compareArea) # lt takes a custom less than function
println("Sorted arr2 by area: ", arr2)

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
using ThreeBodyTB
using JLD2

fillist = readdir()

for f in fillist
    if occursin("jld2",f)
        if occursin("2body" ,f )
            dim = 2
        elseif occursin("3body", f)
            dim = 3
        else
            println("err $f")
        end

        t = @load f
        savename=f[1:end-5]
        c = ThreeBodyTB.CalcTB.make_coefs(dat.names, dim, datH=dat.datH, datS=dat.datS, min_dist = dat.min_dist, maxmin_val_train=dat.maxmin_val_train, dist_frontier=dat.dist_frontier)
        ThreeBodyTB.CalcTB.write_coefs("$savename.xml", c)

    end
end


text:
using Reexport

@reexport module UncertainDatasets
    using ..UncertainValues
    using Distributions
    using IntervalArithmetic
    using RecipesBase
    using StatsBase
    using StaticArrays
    using Statistics

    # The abstract type for all types of dataset holding uncertain values
    include("AbstractUncertainDataset.jl")

    # An abstract type for uncertain datasets containing uncertain values yielding scalar
    # values when resampled.
    include("AbstractUncertainValueDataset.jl")
    include("AbstractUncertainIndexDataset.jl")

    # One composite type for indices, another one for values. This distinction allows more
    # flexibility when applying sampling constraints (some constraints may be meaningful
    # only for indices, for example).
    include("UncertainValueDataset.jl")
    include("UncertainIndexDataset.jl")

    # A generic type with all the functionality of `AbstractUncertainValueDataset`, if you
    # can't be bothered with specifying
    include("UncertainDataset.jl")

    # An abstract type for datasets containing both indices and data values.
    include("AbstractUncertainIndexValueDataset.jl")

    # A composite type with two fields: `indices` and `values`. Both fields may be
    # any subtype of AbstractUncertainValueDataset.
    include("UncertainIndexValueDataset.jl")

    # Conversion and promotion
    include("conversions.jl")

    """
        UVAL_COLLECTION_TYPES = Union{UD, UV} where {
            UD <: AbstractUncertainValueDataset,
            UV <: AbstractVector{T} where {
                T <: AbstractUncertainValue}}

    A type union used to represent types of uncertain values.
    """
    const UVAL_COLLECTION_TYPES = Union{UD, UV} where {
        UD <: AbstractUncertainValueDataset,
        UV <: AbstractVector{T} where {
            T <: AbstractUncertainValue}}

    export UVAL_COLLECTION_TYPES

end # module

"""
	UncertainDatasets

A module defining uncertain datasets, which are collections of uncertain values
defined in the `UncertainValues` module.
"""
UncertainDatasets


text:
module RegionTrees

using StaticArrays
import StaticArrays: similar_type, Size
using IterTools: chain
import Base: show, size, getindex, convert, parent
export Cell,
       split!,
       isleaf,
       findleaf,
       children,
       parent,
       child_boundary,
       HyperRectangle,
       center,
       vertices,
       faces,
       body_and_face_centers,
       AbstractRefinery,
       needs_refinement,
       refine_data,
       initial_data,
       allcells,
       allleaves,
       allparents,
       adaptivesampling!

include("twosarray.jl")
include("hyperrectangle.jl")
include("cell.jl")
include("adaptivesampling.jl")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# macro wrappers for various reflection functions

import Base.typesof

separate_kwargs(args...; kwargs...) = (args, kwargs.data)

function gen_call_with_extracted_types(__module__, fcn, ex0)
    if isa(ex0, Expr)
        if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex0.args)
            return quote
                local arg1 = $(esc(ex0.args[1]))
                local args, kwargs = $separate_kwargs($(map(esc, ex0.args[2:end])...))
                $(fcn)(Core.kwfunc(arg1),
                       Tuple{typeof(kwargs), Core.Typeof(arg1), map(Core.Typeof, args)...})
            end
        elseif ex0.head == :call
            return Expr(:call, fcn, esc(ex0.args[1]),
                        Expr(:call, typesof, map(esc, ex0.args[2:end])...))
        elseif ex0.head == :(.)
            return Expr(:call, fcn, Base.getproperty,
                        Expr(:call, typesof, map(esc, ex0.args)...))
        elseif ex0.head == :(=) && length(ex0.args) == 2 && ex0.args[1].head == :(.)
            return Expr(:call, fcn, Base.setproperty!,
                        Expr(:call, typesof, map(esc, [ex0.args[1].args..., ex0.args[2]])...))
        end
    end
    if isa(ex0, Expr) && ex0.head == :macrocall # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return Expr(:call, fcn, esc(ex0.args[1]), Tuple{#=__source__=#LineNumberNode, #=__module__=#Module, Any[ Core.Typeof(a) for a in ex0.args[3:end] ]...})
    end
    ex = Meta.lower(__module__, ex0)
    exret = Expr(:none)
    if !isa(ex, Expr)
        exret = Expr(:call, :error, "expression is not a function call or symbol")
    elseif ex.head == :call
        if any(e->(isa(e, Expr) && e.head==:(...)), ex0.args) &&
            (ex.args[1] === GlobalRef(Core,:_apply) ||
             ex.args[1] === GlobalRef(Base,:_apply))
            # check for splatting
            exret = Expr(:call, ex.args[1], fcn,
                        Expr(:tuple, esc(ex.args[2]),
                            Expr(:call, typesof, map(esc, ex.args[3:end])...)))
        else
            exret = Expr(:call, fcn, esc(ex.args[1]),
                         Expr(:call, typesof, map(esc, ex.args[2:end])...))
        end
    elseif ex.head == :body
        a1 = ex.args[1]
        if isa(a1, Expr) && a1.head == :call
            a11 = a1.args[1]
            if a11 == :setindex!
                exret = Expr(:call, fcn, a11,
                             Expr(:call, typesof, map(esc, a1.args[2:end])...))
            end
        end
    end
    if ex.head == :thunk || exret.head == :none
        exret = Expr(:call, :error, "expression is not a function call, "
                                  * "or is too complex for @$fcn to analyze; "
                                  * "break it down to simpler parts if possible")
    end
    return exret
end

for fname in [:which, :less, :edit, :functionloc, :code_warntype,
              :code_llvm, :code_llvm_raw, :code_native]
    @eval begin
        macro ($fname)(ex0)
            gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
        end
    end
end
macro which(ex0::Symbol)
    ex0 = QuoteNode(ex0)
    return :(which($__module__, $ex0))
end

for fname in [:code_typed, :code_lowered]
    @eval begin
        macro ($fname)(ex0)
            thecall = gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
            quote
                results = $thecall
                length(results) == 1 ? results[1] : results
            end
        end
    end
end

"""
    @functionloc

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns a tuple `(filename,line)` giving the location for the method that would be called for those arguments.
It calls out to the `functionloc` function.
"""
:@functionloc

"""
    @which

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns the `Method` object for the method that would be called for those arguments. Applied
to a variable, it returns the module in which the variable was bound. It calls out to the
`which` function.
"""
:@which

"""
    @less

Evaluates the arguments to the function or macro call, determines their types, and calls the `less`
function on the resulting expression.
"""
:@less

"""
    @edit

Evaluates the arguments to the function or macro call, determines their types, and calls the `edit`
function on the resulting expression.
"""
:@edit

"""
    @code_typed

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_typed`](@ref) on the resulting expression.
"""
:@code_typed

"""
    @code_lowered

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_lowered`](@ref) on the resulting expression.
"""
:@code_lowered

"""
    @code_warntype

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_warntype`](@ref) on the resulting expression.
"""
:@code_warntype

"""
    @code_llvm

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_llvm`](@ref) on the resulting expression.
"""
:@code_llvm

"""
    @code_native

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_native`](@ref) on the resulting expression.
"""
:@code_native

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
module CSPSolver

include("solve.jl")
include("kenken.jl")
include("sudoku.jl")

export Constraint, CellConstraint, RowConstraint, ColumnConstraint
export GridCSP, KenKen, Sudoku
export CageConstraint  # KenKen
export BoxConstraint   # Sudoku
export consistent, complete, eliminate!, fill_in!, backtrack
export generate_random_sudoku, generate_random_kenken
export print_kenken

end

max_stars_repo_name: gindij/CSPSolver

max_issues_repo_name: gindij/CSPSolver

max_forks_repo_name: gindij/CSPSolver


text:
mutable struct SymbolDFG{T <: AbstractParams, V <: DFGNode, F <:DFGNode} <: AbstractDFG
    g::SymbolFactorGraph{V,F}
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    #NOTE does not exist
    # nodeCounter::Int64
    #NOTE does not exist
    # labelDict::Dict{Symbol, Int64}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
end

#TODO? do we not want props such as userId, robotId, sessionId, etc...
function SymbolDFG{T,V,F}(g::SymbolFactorGraph{V,F}=SymbolFactorGraph{V,F}();
                           description::String="LightGraphs.jl implementation",
                           userId::String="User ID",
                           robotId::String="Robot ID",
                           sessionId::String="Session ID",
                           userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           params::T=NoSolverParams()) where {T <: AbstractParams, V <:DFGNode, F<:DFGNode}

    SymbolDFG{T,V,F}(g, description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], params)
end

SymbolDFG{T}(g::SymbolFactorGraph{DFGVariable,DFGFactor}=SymbolFactorGraph{DFGVariable,DFGFactor}(); kwargs...) where T <: AbstractParams = SymbolDFG{T,DFGVariable,DFGFactor}(g; kwargs...)

Base.propertynames(x::SymbolDFG, private::Bool=false) =
    (:g, :description, :userId, :robotId, :sessionId, :nodeCounter, :labelDict, :addHistory, :solverParams)
        # (private ? fieldnames(typeof(x)) : ())...)

Base.getproperty(x::SymbolDFG,f::Symbol) = begin
    if f == :nodeCounter
        @error "Field nodeCounter deprecated. returning number of nodes"
        nv(x.g)
    elseif f == :labelDict
        @error "Field labelDict deprecated. Consider using exists(dfg,label) or getLabelDict(dfg) instead. Returning internals copy"
        copy(x.g.fadjdict)
    else
        getfield(x,f)
    end
end


text:
@info "Test state pattern..."

g = Gate()
println("The gate is open? ", isopen(g))
unlock!(g)
println("The gate is open? ", isopen(g))
open(g)


text:
function findpeaks(activation; filterlength=10, minh=0, threshold=-Inf, minw=0, maxw=Inf, doplot=false, kwargs...)
    activationf = filtfilt(ones(filterlength),[filterlength],activation)
    doplot && plot!(activationf, lab="Filtered Activation", kwargs...)

    peaks = allpeaks(activationf)
    doplot &&  scatter!(peaks, activationf[peaks], lab="All peaks", kwargs...)

    peaks = peaks[activationf[peaks] .>= minh] # Remove too small peaks
    doplot &&  scatter!(peaks, activationf[peaks]*1.1, lab="Below minh removed", kwargs...)

    peaks = remove_threshold(peaks,activationf,threshold)
    peaks = remove_width(peaks,activationf,minw,maxw)
    doplot &&  scatter!(peaks, activationf[peaks]*1.2, lab="Too wide or narrow removed" , m=(10,:xcross), kwargs...)
    return peaks
end

function remove_threshold(peaks,activation,threshold)
    base = max.(activation[peaks .- 1],activation[peaks .+ 1])
    peaks = peaks[activation[peaks].-base .>= threshold]
end

function remove_width(peaks,activation,minw,maxw)
    if minw <= 0 && maxw == Inf
        return peaks # Nothing needs to be done
    end
    i = 1
    while i < length(peaks)
        speaks = sortperm(activation[peaks], rev=true) # sorted peak locations (in peaks vector)

        if speaks[i] > 1
            wl = abs(peaks[speaks[i]] - peaks[speaks[i]-1]) # highest peak location - location to the left
            if wl < minw || wl > maxw
                deleteat!(peaks,speaks[i]-1)
                continue
            end
        end
        if speaks[i] < length(peaks)
            wr = abs(peaks[speaks[i]] - peaks[speaks[i]+1])
            if wr < minw || wr > maxw
                deleteat!(peaks,speaks[i]+1)
                continue
            end
        end
        i += 1
    end
    if !isempty(peaks)
        @assert minimum(diff(peaks)) >= minw
    end
    peaks
end

function allpeaks(activation)
    da = diff(activation)
    peaks = findall((diff(da) .< 0) .& (da[1:end-1] .> 0) .& (da[2:end] .< 0)) .+ 1
end


text:
module JSObjectLiteral

export @js

"""
@js expression

Tries to parse `expression` as a javascript object literal.

If `expression` is a valid Julia expression, and (almost) looks like a javascript
object literal, it is parsed and converted into a Julia json-like structure.

# Examples
```
e = 5.0
g = "gee!"
j = Dict("k": 1_000)
@js {
  a: 1,
  b: [2, 3],
  c : {
    d: "doubly-quoted string",
    e
  },
  f: g,
  h.i: j.k
}
```
"""
macro js(expr)
    return esc(js(expr))
end

## for expressions, dispatch according to the head
js(expr::Expr) = js(Val{expr.head}, expr)
js(x) = x

## assignment
function js(::Type{Val{:(=)}}, expr)
    if isa(expr.args[1], Expr) && expr.args[1].head == :braces ## { a, b } = dict
        ##  error("Assignment to braces expression not supported yet")
        all(isa(arg, Symbol) for arg in expr.args[1].args) || error("Braced expression as LHS must only contain keys")
        lhs = Expr(:tuple, expr.args[1].args...)
        keys = map(string, expr.args[1].args)
        dict = js(expr.args[2])
        rhs = :([$dict[key] for key in $keys])
        return Expr(:(=), lhs, :($rhs))
    else
        expr.args = map(js, expr.args)
        return expr
    end
end

## key index
function js(::Type{Val{:.}}, expr)
    return Expr(:ref, js(expr.args[1]), string(expr.args[2].value))
end

## array index
function js(::Type{Val{:ref}}, expr)
    expr.args = map(js, expr.args)
    return expr
end

## dictionary creation
function js(::Type{Val{:braces}}, expr)
    dict = :(Dict{String,Any}())
    for arg in expr.args
        push!(dict.args, pair(arg))
    end
    return dict
end

function pair(e::Expr)
    if e.head != :call || length(e.args) !=3 || e.args[1] != :(:)
        error("Expected : operator with 2 arguments")
    end
    return keyvalue(e.args[2], e.args[3])
end
pair(s::Symbol) = :($(string(s)) => $(js(s)))

## keyvalue deals with keys of the type a.b
function keyvalue(key::Expr, value)
    key.head == :. || error("Expected . operator")
    expr = Expr(:braces, Expr(:call, :(:), key.args[2].value, value))
    return keyvalue(key.args[1], expr)
end
keyvalue(key::Symbol, value) = :($(string(key)) => $(js(value)))
keyvalue(key::AbstractString, value) = :($key => $(js(value)))

## array creation
function js(::Type{Val{:vect}}, expr)
    return Expr(:ref, :Any, map(js, expr.args)...) ## make sure arrays are of type Any
end

## cover == amongst others
function js(::Type{Val{:call}}, expr)
    expr.args[2:end] = map(js, expr.args[2:end])
    return expr
end

js(::Type{Val{s}}, expr) where s = expr

function js(x, expr)
    @warn "catchall js(x, expr) " * string(typeof(x)) * " " * string(expr)
    dump(x)
    dump(expr)
    return expr
end

include("JSObject.jl")

end


text:
"""
    applyrule(rules::Chain, data, state, (i, j))

Chained rules. If a `Chain` of rules is passed to applyrule, run them sequentially for each
cell.  This can have much beter performance as no writes occur between rules, and they are
essentially compiled together into compound rules. This gives correct results only for
CellRule, or for a single NeighborhoodRule followed by CellRule.
"""
@inline applyrule(rules::Chain{<:Tuple{<:NeighborhoodRule,Vararg}}, data, state, index, buf) = begin
    state = applyrule(rules[1], data, state, index, buf)
    applyrule(tail(rules), data, state, index)
end
@inline applyrule(rules::Chain, data, state, index) = begin
    state == zero(state) && return state
    newstate = applyrule(rules[1], data, state, index)
    applyrule(tail(rules), data, newstate, index)
end
@inline applyrule(rules::Chain{Tuple{}}, data, state, index) = state


text:
# test various model functions
bigdata = joinpath(dirname(@__FILE__), "..", "data", "big.txt")

TMPDIR = "./tmp"
mkpath(TMPDIR)

CORPUS = bigdata
VOCAB_FILE = joinpath(TMPDIR, "vocab.txt")
COOCCURRENCE_FILE = joinpath(TMPDIR, "cooccurrence.bin")
COOCCURRENCE_SHUF_FILE = joinpath(TMPDIR, "cooccurrence.shuf.bin")
SAVE_FILE = joinpath(TMPDIR, "vectors")
VERBOSE = 0
MEMORY = 4.0
VOCAB_MIN_COUNT = 5
VECTOR_SIZE = 50
MAX_ITER = 2
WINDOW_SIZE = 5
BINARY= 0
NUM_THREADS = 1
X_MAX = 10.0
HEADER = 1

# Create and load model
vocab_count(CORPUS, VOCAB_FILE, min_count=VOCAB_MIN_COUNT, verbose=VERBOSE)
cooccur(CORPUS, VOCAB_FILE, COOCCURRENCE_FILE,
        memory=MEMORY, verbose=VERBOSE, window_size=WINDOW_SIZE)
shuffle(COOCCURRENCE_FILE, COOCCURRENCE_SHUF_FILE,
memory=MEMORY, verbose=VERBOSE)
glove(COOCCURRENCE_SHUF_FILE, VOCAB_FILE, SAVE_FILE,
      threads=NUM_THREADS, x_max=X_MAX, iter=MAX_ITER,
      vector_size=VECTOR_SIZE, binary=BINARY,
      write_header=HEADER, verbose=VERBOSE)
model_file = joinpath(SAVE_FILE) * ifelse(BINARY==1, ".bin", ".txt")
model = wordvectors(model_file, Float32, header=true, kind= :text)

# word vectors
println(model)

len_vecs, num_words = size(model)
wordvecs = model.vectors
@test size(wordvecs) == (len_vecs, num_words)

words = vocabulary(model)
word1 = words[rand(1:end)]
word2 = words[rand(1:end)]
word3 = words[rand(1:end)]

@test in_vocabulary(model, word1)

n = rand(1:100)

indxs, mes = cosine(model, word1, n)
@test words[indxs] == cosine_similar_words(model, word1, n)
w4_indx = indxs[rand(1:end)]
loc = findall(in(w4_indx), indxs)
word4 = words[w4_indx]
@test index(model, word4) == w4_indx

s = similarity(model, word1, word4)
@test mes[loc[1]] ≈ s

inx, mes = analogy(model, [word1, word2], [word3], n)
@test words[inx] == analogy_words(model, [word1, word2], [word3], n)

rm(TMPDIR, recursive=true, force=true)

println("model passed test...")


text:
"""
    radiative_transfer(atm::ModelAtmosphere, α, S, α_ref, mu_grid)

Returns the astrophysical flux at each wavelength.

inputs:
- `atm`: the model atmosphere.
- `α`: a matrix (atmospheric layers × wavelengths) containing the absoprtion coefficient
- `S`: the source fuction as a matrix of the same shape.
- `α_ref`: the continuum absoprtion coefficient at the reference wavelength (5000 Å). Used to
   rescale the total absorption to match the model atmosphere. This value should be calculated by
   Korg.
- `mu_grid`: (required if atm is a [`ShellAtmosphere`](@ref)) the values of μ at which to calculate
   the surface intensity, which is integrated to obtain the astrophysical flux.
"""
function radiative_transfer(atm::PlanarAtmosphere, α, S, α_ref, mu_grid=nothing)
    τ5 = [l.tau_5000 for l in atm.layers] #τ at 5000 Å according to model atmosphere
    planar_transfer(α, S, τ5, α_ref)
end
function radiative_transfer(atm::ShellAtmosphere, α, S, α_ref, mu_grid)
    τ5 = [l.tau_5000 for l in atm.layers] #τ at 5000 Å according to model atmosphere
    radii = [atm.R + l.z for l in atm.layers]
    photosphere_correction = radii[1]^2 / atm.R^2
    #discard I, take F only
    photosphere_correction * spherical_transfer(α, S, τ5, α_ref, radii, mu_grid)[2]
end

"""
    planar_transfer(α, S, τ_ref, α_ref)

Returns the astrophysical flux. See [`radiative_transfer`](@ref) for an explantion of the arguments.
"""
function planar_transfer(α, S, τ_ref, α_ref)
    log_τ_ref = log.(τ_ref)
    τ_ref_α_ref = τ_ref ./ α_ref
    I = Vector{eltype(α)}(undef, size(α, 2))
    τ_λ = Vector{eltype(α)}(undef, size(α, 1))
    for λ_ind in 1:size(α, 2)
        @inbounds cumulative_trapezoid_rule!(τ_λ, log_τ_ref, view(α, :, λ_ind) .* τ_ref_α_ref)
        I[λ_ind] = 2π * all_mu_transfer_integral(τ_λ, view(S, :, λ_ind))
    end
    I
end

"""
    spherical_transfer(α, S, τ_ref, α_ref, radii, μ_surface_grid)

Perform radiative transfer along rays emerging at the μ values in `μ_surface_grid` in a spherically
symmetric atmosphere resolved at radii `radii` [cm]. See [`radiative_transfer`](@ref) for an
explantion of the arguments. Note that `radii` should be in decreasing order.

Returns `(flux, intensity)`, where `flux` is the astrophysical flux, and `intensity`, a matrix of
shape (mu values × wavelengths), is the surface intensity as a function of μ.
"""
function spherical_transfer(α, S, τ_ref, α_ref, radii, μ_surface_grid)
    R, r0 = radii[1], radii[end] #lower bound of atmosphere

    #precompute for use in transfer integral
    log_τ_ref = log.(τ_ref)
    τ_ref_α_ref = τ_ref ./ α_ref

    #general element type with which to preallocate arrays
    el_type = typeof(promote(radii[1], α[1], S[1], μ_surface_grid[1])[1])

    #preallocations
    #geometric path-length correction, and other wavelength-indenpendent integrand factors
    integrand_factor = Matrix{el_type}(undef, size(α, 1), length(μ_surface_grid))
    #the index of the lowest atmospheric layer pierced by each ray
    lowest_layer_indices = Vector{Int}(undef, length(μ_surface_grid))
    for (μ_ind, μ_surface) in enumerate(μ_surface_grid)
        # impact parameter of ray
        b = R * sqrt(1 - μ_surface^2)

        #doing this with `findfirst` is messier at first and last index
        i = argmin(abs.(radii .- b))
        if radii[i] < b
            i -= 1
        end
        lowest_layer_indices[μ_ind] = i

        integrand_factor[1:i, μ_ind] = @. radii[1:i] ./ sqrt(radii[1:i]^2 - b^2) * τ_ref_α_ref[1:i]
    end

    #preallocations
    I = Matrix{el_type}(undef, length(μ_surface_grid), size(α, 2)) #the surface intensity
    integrand = Vector{el_type}(undef, size(α, 1))                 #integrand of τ integral
    τ_λ = Vector{el_type}(undef, size(α, 1))                       #optical depth at a particular λ
    #iterate over λ in the outer loop, μ in the inner loop
    for λ_ind in 1:size(α, 2), μ_ind in 1:length(μ_surface_grid)
        i = lowest_layer_indices[μ_ind]
        for k in 1:i #I can't figure out how to write this as a fast one-liner
            integrand[k] = α[k, λ_ind] * integrand_factor[k, μ_ind]
        end
        cumulative_trapezoid_rule!(τ_λ, log_τ_ref, integrand, i) #compute τ_λ
        I[μ_ind, λ_ind] = ray_transfer_integral(view(τ_λ,1:i), view(S,1:i,λ_ind))

        #this branch could be factored out of this loop, which might speed things up.
        if i < length(radii)
            #if the ray never leaves the model atmosphere, include the contribution from the
            #other side of the star.  Calculate the optical depths you get by reversing the τ_λ.
            #This is less accurate than actually integrating to find τ, but the effect is small.
            #This should probably by audited for off-by-one errors.  It's also inneficient.
            τ_prime = τ_λ[i] .+ [cumsum(reverse(diff(view(τ_λ,1:i)))) ; 0]
            I[μ_ind, λ_ind] += ray_transfer_integral(view(τ_prime, 1:i), view(S,1:i,λ_ind))
        else #otherwise assume I=S at atmosphere lower boundary.  This is a _tiny_ effect.
            I[μ_ind, λ_ind] += exp(-τ_λ[end]) * S[end, λ_ind]
        end
    end
    #calculate 2π∫μIdμ to get astrophysical flux
    F = 2π * [Korg.trapezoid_rule(μ_surface_grid, μ_surface_grid .* I) for I in eachcol(I)]
    I, F
end

"""
    all_mu_transfer_integral(τ, S)

Compute exactly the solution to the transfer integral obtained be linearly interpolating the source
function, `S` across optical depths `τ`, without approximating the factor of E₂(τ).
"""
function all_mu_transfer_integral(τ, S)
    I = 0
    for i in 1:length(τ)-1
        @inbounds m = (S[i+1] - S[i])/(τ[i+1] - τ[i])
        @inbounds b = S[i] - m*τ[i]
        @inbounds I += (_plane_parallel_approximate_transfer_integral(τ[i+1], m, b) -
                        _plane_parallel_approximate_transfer_integral(τ[i], m, b))
    end
    I
end

"""
    ray_transfer_integral(τ, S)

Compute exactly the solution to the transfer integral obtained be linearly interpolating the source
function, `S` across optical depths `τ`, without approximating the factor of exp(-τ).

This breaks the integral into the sum of integrals of the form
\$\\int (m\\tau + b) \\exp(-\\tau)\$ d\\tau\$ ,
which is equal to
\$ -\\exp(-\\tau) (m*\\tau + b + m)\$.
"""
function ray_transfer_integral(τ, S)
    if length(τ) == 1
        return 0.0
    end
    I = 0.0
    next_exp_negτ = exp(-τ[1])
    for i in 1:length(τ)-1
        @inbounds m = (S[i+1] - S[i])/(τ[i+1] - τ[i])
        cur_exp_negτ = next_exp_negτ
        @inbounds next_exp_negτ = exp(-τ[i+1])
        @inbounds I += (-next_exp_negτ * (S[i+1] + m) + cur_exp_negτ * (S[i] + m))
    end
    I
end

"""
    _plane_parallel_approximate_transfer_integral(τ, m, b)

The exact solution to \$\\int (m\\tau + b) E_2(\\tau)\$ d\\tau\$.

The exponential integral function, expint, captures the integral over the disk of the star to
get the emergent astrophysical flux. You can verify it by substituting the variable of integration
in the exponential integal, t, with mu=1/t.
"""
function _plane_parallel_approximate_transfer_integral(τ, m, b)
    1/6 * (τ*exponential_integral_2(τ)*(3b+2m*τ) - exp(-τ)*(3b + 2m*(τ+1)))
end

"""
    trapezoid_rule(xs, fs)

Approximate the integral f(x) with the trapezoid rule over x-values `xs` given f(x) values `fs`.
"""
function trapezoid_rule(xs, fs)
    Δs = diff(xs)
    weights = [0 ; Δs] + [Δs ; 0]
    sum(0.5 * weights .* fs)
end

"""
    cumulative_trapezoid_rule!(out, xs, fs, [len=length(xs)])

Approximate the closed integral of f(x) from the first element of `xs` to each element of `xs`
(up to `len`) with the trapezoid rule, given f values `fs`.

Assigns to the preallocated vector `out`.
"""
function cumulative_trapezoid_rule!(out, xs, fs, len=length(xs))
    out[1] = 0.0
    for i in 2:len
        out[i] = out[i-1] + 0.5*(fs[i]+fs[i-1])*(xs[i]-xs[i-1])
    end
end

"""
    exponential_integral_2(x)

Approximate second order exponential integral, E_2(x).  This stiches together several series
expansions to get an approximation which is accurate within 1% for all `x`.
"""
function exponential_integral_2(x)
    if x == 0
        0.0
    elseif x < 1.1
        _expint_small(x)
    elseif x < 2.5
        _expint_2(x)
    elseif x < 3.5
        _expint_3(x)
    elseif x < 4.5
        _expint_4(x)
    elseif x < 5.5
        _expint_5(x)
    elseif x < 6.5
        _expint_6(x)
    elseif x < 7.5
        _expint_7(x)
    elseif x < 9
        _expint_8(x)
    else
        _expint_large(x)
    end
end

function _expint_small(x)
    #euler mascheroni constant
    ℇ = 0.57721566490153286060651209008240243104215933593992
    1 + ((log(x) + ℇ - 1) + (-0.5 + (0.08333333333333333 + (-0.013888888888888888 +
                                                            0.0020833333333333333*x)*x)*x)*x)*x
end
function _expint_large(x)
    invx = 1/x
    exp(-x) * (1 + (-2 + (6 + (-24 + 120*invx)*invx)*invx)*invx)*invx
end
function _expint_2(x)
    x -= 2
    0.037534261820486914 + (-0.04890051070806112 + (0.033833820809153176 + (-0.016916910404576574 +
                                          (0.007048712668573576 -0.0026785108140579598*x)*x)*x)*x)*x
end
function _expint_3(x)
    x -= 3
    0.010641925085272673   + (-0.013048381094197039   + (0.008297844727977323   +
            (-0.003687930990212144   + (0.0013061422257001345  - 0.0003995258572729822*x)*x)*x)*x)*x
end
function _expint_4(x)
    x -= 4
    0.0031982292493385146  + (-0.0037793524098489054  + (0.0022894548610917728  +
            (-0.0009539395254549051  + (0.00031003034577284415 - 8.466213288412284e-5*x )*x)*x)*x)*x
end
function _expint_5(x)
    x -= 5
    0.000996469042708825   + (-0.0011482955912753257  + (0.0006737946999085467  +
            (-0.00026951787996341863 + (8.310134632205409e-5   - 2.1202073223788938e-5*x)*x)*x)*x)*x
end
function _expint_6(x)
    x -= 6
    0.0003182574636904001  + (-0.0003600824521626587  + (0.00020656268138886323 +
            (-8.032993165122457e-5   + (2.390771775334065e-5   - 5.8334831318151185e-6*x)*x)*x)*x)*x
end
function _expint_7(x)
    x -= 7
    0.00010350984428214624 + (-0.00011548173161033826 + (6.513442611103688e-5   +
            (-2.4813114708966427e-5  + (7.200234178941151e-6   - 1.7027366981408086e-6*x)*x)*x)*x)*x
end
function _expint_8(x)
    x -= 8
    3.413764515111217e-5   + (-3.76656228439249e-5    + (2.096641424390699e-5   +
            (-7.862405341465122e-6   + (2.2386015208338193e-6  - 5.173353514609864e-7*x )*x)*x)*x)*x
end


text:
import Pkg
Pkg.activate("..")

using Distributions: Exponential, DiscreteUniform, Poisson, Normal, DiscreteNonParametric, truncated
using CrimsonSkyline
using Random: seed!
using Logging
using StatsBase: mean

seed!(2021)

# source: http://people.reed.edu/~jones/141/Coal.html
const coal_mining_data = [
    4, 5, 4, 1, 0, 4, 3, 4, 0, 6, 3, 3,
    4, 0, 2, 6, 3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5, 2,
    2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0, 1, 0, 1, 1, 0, 0,
    3, 1, 0, 3, 2, 2, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0,
    0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 0, 1,
    4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1
]

function changepoint_model(t, data)
    rate_1 = sample(t, "rate_1", Exponential(5.0))
    rate_2 = sample(t, "rate_2", Exponential(5.0))
    changepoint = sample(t, "changepoint", DiscreteUniform(2, length(data) - 1))
    for (i, d) in enumerate(data)
        i < changepoint ? observe(t, "data $i", Poisson(rate_1), d) : observe(t, "data $i", Poisson(rate_2), d)
    end
end

rate_1_proposal(t0, t1, params...) = propose(t1, "rate_1", truncated(Normal(t0["rate_1"].value, 0.5), 0.01, Inf))
rate_2_proposal(t0, t1, params...) = propose(t1, "rate_2", truncated(Normal(t0["rate_2"].value, 0.5), 0.01, Inf))
function cp_proposal(t0, t1, data)
    ld = length(data) - 1
    pt = t0["changepoint"].value
    if pt == 2
        propose(t1, "changepoint", DiscreteNonParametric([2, 3], [0.5, 0.5]))
    elseif pt == ld
        propose(t1, "changepoint", DiscreteNonParametric([ld - 1, ld], [0.5, 0.5]))
    else
        propose(t1, "changepoint", DiscreteNonParametric([pt - 1, pt, pt + 1], [1/3.0, 1/3.0, 1/3.0]))
    end
end

function main()
    @time results = mh(
        changepoint_model, [rate_1_proposal, rate_2_proposal, cp_proposal];
        params = (coal_mining_data,), burn = 1000, thin = 100, num_iterations = 11000, inverse_verbosity = 1000
    )
    ests = hpdi(results, 0.9, ["rate_1", "rate_2", "changepoint"])
    @info "90% posterior HPDIs: $ests"
    plot_marginal(results, "changepoint", "plots/changepoint", "changepoint.png")
end

main()


text:
using Plots
using DelaunayMaps

m  = Map(20, 20, 3);
plot_map(m, window=:fit)
savefig("../images/example_3_points.png")
m  = Map(20, 20, 10);
plot_map(m, show_nums=false, window=:fit)
savefig("../images/example_10_points.png")
m  = Map(20, 20, 100);
plot_map(m, show_nums=false, window=:fit)
savefig("../images/example.png")


text:
module GMeansController
using DataFrames, Query, Combinatorics
using ..GMeans
using JSON

export call

"""
call(message)

Computes g-means from records.

# Arguments
- `message::Dict`: the arguments dict
  - `"records"`: the array of array of rows

# Return
the json string of result object

# Examples
```jldoctest
julia> GMeansController.call(Dict("records" => [[1.0, 2.0, 0.0], [2.0, 1.0, 2.0], [3.0, 0.0, 9.0]]))
"{\"centers\":[[3.0,0.0,9.0],[2.0,1.0,2.0],[1.0,2.0,0.0]],\"clusters\":[[[3.0,0.0,9.0]],[[2.0,1.0,2.0]],[[1.0,2.0,0.0]]]}"
```
"""
function call(message)
  gmeans_result = message["records"] |> collect |> GMeans.call

  Dict(
      :clusters => gmeans_result[:clusters],
      :centers => gmeans_result[:centers]
  ) |> JSON.json
end
end