text:
using Printf
using BenchmarkTools
function heat()
    N = 1001
    T0 = Matrix{Float64}(undef,N,N)
    T1 = Matrix{Float64}(undef,N,N)
    x  = Matrix{Float64}(undef,N,N)
    y  = Matrix{Float64}(undef,N,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    for j=1:N
        for i=1:N
            x[i,j] = (i-1)*dx
        end
    end
    for j=1:N
        for i=1:N
            y[i,j] = (j-1)*dx
        end
    end

    dt = dx
    tmax = (10-mod(10,dt))/dt
    k = 0.25*dx
    alpha = k*dt/dx^2
    for i=1:N
        T1[1,i] = 1.0
        T1[N,i] = 1.0
        T1[i,1] = 1.0
        T1[i,N] = 1.0
    end
    for i=1:N
        T0[1,i] = 1.0
        T0[N,i] = 1.0
        T0[i,1] = 1.0
        T0[i,N] = 1.0
    end
    t = 0*dt
    vel = open("Julia.dat", "w")
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"y\", \"T\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
    for j=1:N
        for i=1:N
            write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
        end
    end

    for k=1:tmax
        for j=2:N-1
            for i=2:N-1
                T1[i,j] = T0[i,j]+alpha*((T0[i+1,j]-2*T0[i,j]+T0[i-1,j])+(T0[i,j+1]-2*T0[i,j]+T0[i,j-1]))
            end
        end
        T0 .= T1
        if mod(k,500)==0
            write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
            write(vel, "variables = \"x\", \"y\", \"T\""," \n")
            write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
            for j=1:N
                for i=1:N
                    write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
                end
            end
        end
    end
    t = tmax*dt
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"y\", \"T\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
    for j=1:N
        for i=1:N
            write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
        end
    end
    close(vel)
end

function heat_Vec()
    N = 1001
    T0 = Matrix{Float64}(undef,N,N)
    T1 = Matrix{Float64}(undef,N,N)
    x  = Matrix{Float64}(undef,N,N)
    y  = Matrix{Float64}(undef,N,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x[:,1] = (0:N-1)*dx
    x[:,:] .= x[:,1]
    y[1,:] = (0:N-1)*dx
    y[:,:] .= y[1,:]

    dt = dx
    tmax = (10-mod(10,dt))/dt
    k = 0.25*dx
    alpha = k*dt/dx^2
    T1[1,:] .= 1.0
    T1[N,:] .= 1.0
    T1[:,1] .= 1.0
    T1[:,N] .= 1.0

    T0[1,:] .= 1.0
    T0[N,:] .= 1.0
    T0[:,1] .= 1.0
    T0[:,N] .= 1.0
    t = 0*dt
    vel = open("Julia.dat", "w")
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"y\", \"T\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
    for j=1:N
        for i=1:N
            write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
        end
    end

    for k=1:tmax
        T1[2:N-1,2:N-1] .= T0[2:N-1,2:N-1].+alpha.*((T0[3:N,2:N-1].-2 .*T0[2:N-1,2:N-1].+T0[1:N-2,2:N-1]).+
                                                    (T0[2:N-1,3:N].-2 .*T0[2:N-1,2:N-1].+T0[2:N-1,1:N-2]))
        T0 .= T1
        if mod(k,500)==0
            write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
            write(vel, "variables = \"x\", \"y\", \"T\""," \n")
            write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
            for j=1:N
                for i=1:N
                    write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
                end
            end
        end
    end
    t = tmax*dt
    write(vel, "title =\"ZoneTime_",string(t),"\""," \n")
    write(vel, "variables = \"x\", \"y\", \"T\""," \n")
    write(vel, "zone T=\"Zone_"*string(t)*"\" i=",@sprintf("%d",N)," j=",@sprintf("%d",N)," \n")
    for j=1:N
        for i=1:N
            write(vel, @sprintf("%1.9e",x[i,j])," ",@sprintf("%1.9e",y[i,j])," ",@sprintf("%1.9e",T0[i,j])," \n")
        end
    end
    close(vel)
end

function heat_wof()
    N = 1001
    T0 = Matrix{Float64}(undef,N,N)
    T1 = Matrix{Float64}(undef,N,N)
    x  = Matrix{Float64}(undef,N,N)
    y  = Matrix{Float64}(undef,N,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    for j=1:N
        for i=1:N
            x[i,j] = (i-1)*dx
        end
    end
    for j=1:N
        for i=1:N
            y[i,j] = (j-1)*dx
        end
    end

    dt = dx
    tmax = (10-mod(10,dt))/dt
    k = 0.25*dx
    alpha = k*dt/dx^2
    for i=1:N
        T1[1,i] = 1.0
        T1[N,i] = 1.0
        T1[i,1] = 1.0
        T1[i,N] = 1.0
    end
    for i=1:N
        T0[1,i] = 1.0
        T0[N,i] = 1.0
        T0[i,1] = 1.0
        T0[i,N] = 1.0
    end

    for k=1:tmax
        for j=2:N-1
            for i=2:N-1
                T1[i,j] = T0[i,j]+alpha*((T0[i+1,j]-2*T0[i,j]+T0[i-1,j])+(T0[i,j+1]-2*T0[i,j]+T0[i,j-1]))
            end
        end
        T0 .= T1
    end
end

function heat_Vec_wof()
    N = 1001
    T0 = Matrix{Float64}(undef,N,N)
    T1 = Matrix{Float64}(undef,N,N)
    x  = Matrix{Float64}(undef,N,N)
    y  = Matrix{Float64}(undef,N,N)
    a = 0
    b = π;
    dx = (b-a)/(N-1)
    x[:,1] = (0:N-1)*dx
    x[:,:] .= x[:,1]
    y[1,:] = (0:N-1)*dx
    y[:,:] .= y[1,:]

    dt = dx
    tmax = (10-mod(10,dt))/dt
    k = 0.25*dx
    alpha = k*dt/dx^2
    T1[1,:] .= 1.0
    T1[N,:] .= 1.0
    T1[:,1] .= 1.0
    T1[:,N] .= 1.0

    T0[1,:] .= 1.0
    T0[N,:] .= 1.0
    T0[:,1] .= 1.0
    T0[:,N] .= 1.0

    for k=1:tmax
        T1[2:N-1,2:N-1] .= T0[2:N-1,2:N-1].+alpha.*((T0[3:N,2:N-1].-2 .*T0[2:N-1,2:N-1].+T0[1:N-2,2:N-1]).+
                                                    (T0[2:N-1,3:N].-2 .*T0[2:N-1,2:N-1].+T0[2:N-1,1:N-2]))
        T0 .= T1
    end

end

results1 = @benchmark heat()
results2 = @benchmark heat_Vec()
results3 = @benchmark heat_wof()
results4 = @benchmark heat_Vec_wof()


text:
# Load the packages we need
using JLD2, Plots, DataFrames

# Load the parameters data frame
@load "params.jld2" params

# Create empty arrays that will hold the R2 values
R2_train = []
R2_val = []

# Let's load all the results we want, and put them into the arrays
for i = 1:nrow(params)
	# load the results of experiments i

end

# Add the R2 values to the dataframe

# Subset the dataframe to only some of the results

# Make the plot

Plots.savefig("plot.pdf")

ext: jl


text:
using Test
using LDPCStorage

@testset "LDPCStorage" begin

    all_tests = Dict{String, Array{String, 1}}(
        "File Formats" => ["test_alist.jl",
                           "test_cscmat.jl",
                           ],
        )

    for (testsetname, test_files) in all_tests
        @testset "$testsetname"  begin
            @info "Running tests for `$testsetname`\n$(repeat("=", 60))\n"
            for source in test_files
                @testset "$source" begin
                    @info "Running tests in `$source`..."
                    include(source)
                end
            end
        end
    end

end  # @testset "LDPCStorage"


text:
function funSimSimpleLDPP(eigenVectL,eigenValL)

eigenVectK = (eigenVectL);
eigenValK = eigenValL./(1 .+eigenValL); #eigenvalues of K

#Bernoulli trials (ie coin flips) to determine number of points
booleEigen = (rand(length(eigenValK)) .<= eigenValK);

#number of points in the DPP realization
numbPointsDPP=sum(booleEigen);
#retrieve eigenvectors corresponding to successful Bernoulli trials
global spaceV = eigenVectK[:,booleEigen]; #subspace V
indexDPP = zeros(Int8,numbPointsDPP); #index for final DPP configuration

global numbPointsRemain=numbPointsDPP;
#Loop through for all points
for ii=1:numbPointsDPP
    #Compute probabilities for each point i
    global Prob_i=vec(sum(spaceV.^2,dims=2)); #sum across rows
    Prob_i=Prob_i./sum(Prob_i); #normalize

    #Choose a point (from 1 to numbPoints) using (prob mass function) Prob_i
    uRand=rand(1);
    indexCurrent= findall(x -> x> uRand[1], cumsum(Prob_i))[1];
    indexDPP[ii]=indexCurrent; #update index

    if ii<numbPointsDPP
        #Choose a vector to remove
        jj=findall(x-> x>0, abs.(spaceV[indexCurrent,:]))[1];
        columnVj=spaceV[:,jj]; #j-th column of V
        booleKeep_j=(jj).!=(1:numbPointsRemain);
        global spaceV=spaceV[:,booleKeep_j]; #remove column
        #probably a better way for removing columns

        #Update matrix V by removing Vj component from the space
        spaceV=spaceV-kron(columnVj, transpose(spaceV[indexCurrent,:]./columnVj[indexCurrent]));

        #Orthonormalize
        decompQR=qr(spaceV);
        spaceV=Matrix(decompQR.Q);
        #need Matrix function for properly sized matrix
    end

    global numbPointsRemain= numbPointsRemain-1; #update
end
return indexDPP=sort(indexDPP);
end


text:
# https://github.com/nolta/Winston.jl

# credit goes to https://github.com/jverzani for contributing to the first draft of this backend implementation

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

## dictionaries for conversion of Plots.jl names to Winston ones.
@compat const winston_linestyle =  Dict(:solid=>"solid",
                                :dash=>"dash",
                                :dot=>"dotted",
                                :dashdot=>"dotdashed"
                               )

@compat const winston_marker = Dict(:none=>".",
                            :rect => "square",
                            :ellipse=>"circle",
                            :diamond=>"diamond",
                            :utriangle=>"triangle",
                            :dtriangle=>"down-triangle",
                            :cross => "plus",
                            :xcross => "cross",
                            :star5 => "asterisk"
                           )

function _before_add_series(plt::Plot{WinstonPackage})
  Winston.ghf(plt.o)
end

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

function _create_plot(pkg::WinstonPackage; kw...)
  d = Dict(kw)
  wplt = Winston.FramedPlot(title = d[:title], xlabel = d[:xlabel], ylabel = d[:ylabel])

  Plot(wplt, pkg, 0, d, Dict[])
end

copy_remove(d::Dict, s::Symbol) = delete!(copy(d), s)

function addRegressionLineWinston(d::Dict, wplt)
  xs, ys = regressionXY(d[:x], d[:y])
  Winston.add(wplt, Winston.Curve(xs, ys, kind="dotted"))
end

function getWinstonItems(plt::Plot)
  if isa(plt.o, Winston.FramedPlot)
    wplt = plt.o
    window, canvas = nothing, nothing
  else
    window, canvas, wplt = plt.o
  end
  window, canvas, wplt
end

function _add_series(::WinstonPackage, plt::Plot; kw...)
  d = Dict(kw)

  window, canvas, wplt = getWinstonItems(plt)

  # until we call it normally, do the hack
  if d[:linetype] == :bar
    d = barHack(;d...)
  end

  e = Dict()
  e[:color] = getColor(d[:linecolor])
  e[:linewidth] = d[:linewidth]
  e[:kind] = winston_linestyle[d[:linestyle]]
  e[:symbolkind] = winston_marker[d[:markershape]]
  # markercolor     # same choices as `color`, or :match will set the color to be the same as `color`
  e[:symbolsize] = d[:markersize] / 5

  # pos             # (Int,Int), move the enclosing window to this position
  # screen          # Integer, move enclosing window to this screen number (for multiscreen desktops)

  ## lintype :path, :step, :stepinverted, :sticks, :dots, :none, :heatmap, :hexbin, :hist, :bar
  if d[:linetype] == :none
    Winston.add(wplt, Winston.Points(d[:x], d[:y]; copy_remove(e, :kind)..., color=getColor(d[:markercolor])))

  elseif d[:linetype] == :path
    x, y = d[:x], d[:y]
    Winston.add(wplt, Winston.Curve(x, y; e...))

    fillrange = d[:fillrange]
    if fillrange != nothing
      if isa(fillrange, AbstractVector)
        y2 = fillrange
      else
        y2 = Float64[fillrange for yi in y]
      end
      Winston.add(wplt, Winston.FillBetween(x, y, x, y2, fillcolor=getColor(d[:fillcolor])))
    end

  elseif d[:linetype] == :scatter
    if d[:markershape] == :none
      d[:markershape] = :ellipse
    end

  # elseif d[:linetype] == :step
  #     fn = Winston.XXX

  # elseif d[:linetype] == :stepinverted
  #     fn = Winston.XXX

  elseif d[:linetype] == :sticks
      Winston.add(wplt, Winston.Stems(d[:x], d[:y]; e...))

  # elseif d[:linetype] == :dots
  #     fn = Winston.XXX

  # elseif d[:linetype] == :heatmap
  #     fn = Winston.XXX

  # elseif d[:linetype] == :hexbin
  #     fn = Winston.XXX

  elseif d[:linetype] == :hist
      hst = hist(d[:y], d[:nbins])
      Winston.add(wplt, Winston.Histogram(hst...; copy_remove(e, :nbins)...))

  # elseif d[:linetype] == :bar
  #     # fn = Winston.XXX

  else
    error("linetype $(d[:linetype]) not supported by Winston.")

  end

  # markershape
  if d[:markershape] != :none
    Winston.add(wplt, Winston.Points(d[:x], d[:y]; copy_remove(e, :kind)..., color=getColor(d[:markercolor])))
  end

  # optionally add a regression line
  d[:smooth] && d[:linetype] != :hist && addRegressionLineWinston(d, wplt)

  push!(plt.seriesargs, d)
  plt
end

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

@compat const _winstonNames = Dict(
    :xlims => :xrange,
    :ylims => :yrange,
    :xscale => :xlog,
    :yscale => :ylog,
  )

function _update_plot(plt::Plot{WinstonPackage}, d::Dict)
  window, canvas, wplt = getWinstonItems(plt)
  for k in (:xlabel, :ylabel, :title, :xlims, :ylims)
    if haskey(d, k)
      Winston.setattr(wplt, string(get(_winstonNames, k, k)), d[k])
    end
  end

  for k in (:xscale, :yscale)
    if haskey(d, k)
      islogscale = d[k] == :log10
      Winston.setattr(wplt, (k == :xscale ? :xlog : :ylog), islogscale)
    end
  end

end

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

function createWinstonAnnotationObject(plt::Plot{WinstonPackage}, x, y, val::@compat(AbstractString))
  Winston.text(x, y, val)
end

function _add_annotations{X,Y,V}(plt::Plot{WinstonPackage}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    createWinstonAnnotationObject(plt, ann...)
  end
end

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

function _create_subplot(subplt::Subplot{WinstonPackage}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
end

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

function addWinstonLegend(plt::Plot, wplt)
  if plt.plotargs[:legend]
    Winston.legend(wplt, [sd[:label] for sd in plt.seriesargs])
  end
end

function Base.writemime(io::IO, ::MIME"image/png", plt::PlottingObject{WinstonPackage})
  window, canvas, wplt = getWinstonItems(plt)
  addWinstonLegend(plt, wplt)
  writemime(io, "image/png", wplt)
end

function Base.display(::PlotsDisplay, plt::Plot{WinstonPackage})

  window, canvas, wplt = getWinstonItems(plt)

  if window == nothing
    if Winston.output_surface != :gtk
      error("Gtk is the only supported display for Winston in Plots.  Set `output_surface = gtk` in src/Winston.ini")
    end
    # initialize window
    w,h = plt.plotargs[:size]
    canvas = Gtk.GtkCanvasLeaf()
    window = Gtk.GtkWindowLeaf(canvas, plt.plotargs[:windowtitle], w, h)
    plt.o = (window, canvas, wplt)
  end

  addWinstonLegend(plt, wplt)

  Winston.display(canvas, wplt)
  Gtk.showall(window)
end

function Base.display(::PlotsDisplay, subplt::Subplot{WinstonPackage})
  # TODO: display/show the Subplot object
end


text:
# IO
# --

function readbytes!(source::IO, buffer::AbstractArray{UInt8}, from::Int, to::Int)
    return Base.readbytes!(source, view(buffer, from:to), to - from + 1)
end

function writebytes(sink::IO, buffer::AbstractArray{UInt8}, n::Int, eof::Bool)
    return write(sink, view(buffer, 1:n))
end

# IOStream
# --------

function readbytes!(source::IOStream, buffer::AbstractArray{UInt8}, from::Int, to::Int)
    return Int(ccall(:ios_readall, UInt, (Ptr{Void}, Ptr{Void}, UInt), source.ios,
                     pointer(buffer, from), to - from + 1))
end

# TODO: using ios_write, but look into _os_write_all
#function Base.write(source::IOStream, buffer::Vector{UInt8}, n::Int)
#end


text:
# Code for video 3_3:

Pkg.add("RDatasets")
using RDatasets
data = dataset("datasets", "iris")
head(data)
# 6x5 DataFrames.DataFrame
# | Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species  |
# |-----|-------------|------------|-------------|------------|----------|
# | 1   | 5.1         | 3.5        | 1.4         | 0.2        | "setosa" |
# | 2   | 4.9         | 3.0        | 1.4         | 0.2        | "setosa" |
# | 3   | 4.7         | 3.2        | 1.3         | 0.2        | "setosa" |
# | 4   | 4.6         | 3.1        | 1.5         | 0.2        | "setosa" |
# | 5   | 5.0         | 3.6        | 1.4         | 0.2        | "setosa" |
# | 6   | 5.4         | 3.9        | 1.7         | 0.4        | "setosa" |

summary(data)
"150x5 DataFrames.DataFrame"

features = convert(Array, data[:, 1:4]);
show(features)
150x4 Array{Float64,2}:
 5.1  3.5  1.4  0.2
 4.9  3.0  1.4  0.2
 4.7  3.2  1.3  0.2

labels = convert(Array, data[:, 5]);
labels
150-element Array{ASCIIString,1}:
 "setosa"
 "setosa"
 ...

data[:SepalWidth]
# 150-element DataArrays.DataArray{Float64,1}:
#  3.5
#  3.0
#  3.2
#  ...

julia> eltype(data[:SepalWidth])
Float64

data[2]
150-element DataArrays.DataArray{Float64,1}:
 3.5
 3.0
 3.2
...

data_sub = data[end-9:end,2:3]
10x2 DataFrames.DataFrame
| Row | SepalWidth | PetalLength |
|-----|------------|-------------|
| 1   | 3.1        | 5.6         |
| 2   | 3.1        | 5.1         |
| 3   | 2.7        | 5.1         |
| 4   | 3.2        | 5.9         |
| 5   | 3.3        | 5.7         |
| 6   | 3.0        | 5.2         |
| 7   | 2.5        | 5.0         |
| 8   | 3.0        | 5.2         |
| 9   | 3.4        | 5.4         |
| 10  | 3.0        | 5.1         |

data_matrix = convert(Matrix, data)
150x5 Array{Any,2}:
 5.1  3.5  1.4  0.2  "setosa"
 4.9  3.0  1.4  0.2  "setosa"
 4.7  3.2  1.3  0.2  "setosa"
 4.6  3.1  1.5  0.2  "setosa"
 5.0  3.6  1.4  0.2  "setosa"
 5.4  3.9  1.7  0.4  "setosa"
 4.6  3.4  1.4  0.3  "setosa"
 5.0  3.4  1.5  0.2  "setosa"
 4.4  2.9  1.4  0.2  "setosa"

data_matrix = convert(Matrix, data[1:2])
150x2 Array{Float64,2}:
 5.1  3.5
 4.9  3.0
 4.7  3.2
 4.6  3.1
 5.0  3.6
 5.4  3.9
 4.6  3.4

sort(unique(data[:SepalWidth]))
23-element DataArrays.DataArray{Float64,1}:
 2.0
 2.2
 2.3
 2.4
 2.5
 ...
 2.6
 2.7
 2.8
 2.9
 3.0
 3.1
 3.2
 3.3
 3.4
 3.5
 3.6
 3.7
 3.8
 3.9
 4.0
 4.1
 4.2
 4.4

data[ data[:SepalWidth].==3.2, : ]
# 13x5 DataFrames.DataFrame
# | Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species      |
# |-----|-------------|------------|-------------|------------|--------------|
# | 1   | 4.7         | 3.2        | 1.3         | 0.2        | "setosa"     |
# | 2   | 4.7         | 3.2        | 1.6         | 0.2        | "setosa"     |
# | 3   | 5.0         | 3.2        | 1.2         | 0.2        | "setosa"     |
# | 4   | 4.4         | 3.2        | 1.3         | 0.2        | "setosa"     |
# | 5   | 4.6         | 3.2        | 1.4         | 0.2        | "setosa"     |
# | 6   | 7.0         | 3.2        | 4.7         | 1.4        | "versicolor" |
# | 7   | 6.4         | 3.2        | 4.5         | 1.5        | "versicolor" |
# | 8   | 5.9         | 3.2        | 4.8         | 1.8        | "versicolor" |
# | 9   | 6.5         | 3.2        | 5.1         | 2.0        | "virginica"  |
# | 10  | 6.4         | 3.2        | 5.3         | 2.3        | "virginica"  |
# | 11  | 6.9         | 3.2        | 5.7         | 2.3        | "virginica"  |
# | 12  | 7.2         | 3.2        | 6.0         | 1.8        | "virginica"  |
# | 13  | 6.8         | 3.2        | 5.9         | 2.3        | "virginica"  |

colwise(mean,data[1:4])
4-element Array{Any,1}:
 [5.843333333333334]
 [3.0573333333333332]
 [3.7579999999999996]
 [1.1993333333333331]

for row in eachrow(data[1:4])
    println(mean(convert(Array,row)))
end
2.55
2.375
2.35
2.3499999999999996
...

by(data, :Species, size)
3x2 DataFrames.DataFrame
| Row | Species      | x1     |
|-----|--------------|--------|
| 1   | "setosa"     | (50,5) |
| 2   | "versicolor" | (50,5) |
| 3   | "virginica"  | (50,5) |

by(data, :Species, df -> mean(df[:SepalLength]))
3x2 DataFrames.DataFrame
| Row | Species      | x1    |
|-----|--------------|-------|
| 1   | "setosa"     | 5.006 |
| 2   | "versicolor" | 5.936 |
| 3   | "virginica"  | 6.588 |

julia> by(data, :Species, df -> DataFrame(μ = mean(df[:SepalLength])))
3x2 DataFrames.DataFrame
| Row | Species      | μ     |
|-----|--------------|-------|
| 1   | "setosa"     | 5.006 |
| 2   | "versicolor" | 5.936 |
| 3   | "virginica"  | 6.588 |

julia> by(data, :Species) do df
           DataFrame(μ = mean(df[:SepalLength]), σ² = var(df[:SepalLength]))
       end
3x3 DataFrames.DataFrame
| Row | Species      | μ     | σ²       |
|-----|--------------|-------|----------|
| 1   | "setosa"     | 5.006 | 0.124249 |
| 2   | "versicolor" | 5.936 | 0.266433 |
| 3   | "virginica"  | 6.588 | 0.404343 |

aggregate(data, :Species, mean)
3x5 DataFrames.DataFrame
| Row | Species      | SepalLength_mean | SepalWidth_mean | PetalLength_mean |
|-----|--------------|------------------|-----------------|------------------|
| 1   | "setosa"     | 5.006            | 3.428           | 1.462            |
| 2   | "versicolor" | 5.936            | 2.77            | 4.26             |
| 3   | "virginica"  | 6.588            | 2.974           | 5.552            |

| Row | PetalWidth_mean |
|-----|-----------------|
| 1   | 0.246           |
| 2   | 1.326           |
| 3   | 2.026           |

for subdf in groupby(data, :Species)
    println(size(subdf, 1))
    @printf "%6.3f\t" mean(subdf[:SepalLength])
    @printf "%6.3f\t" mean(subdf[:SepalWidth])
    @printf "%6.3f\t" mean(subdf[:PetalLength])
    @printf "%6.3f\n" mean(subdf[:PetalWidth])
end
50
 5.006   3.428   1.462   0.246
50
 5.936   2.770   4.260   1.326
50
 6.588   2.974   5.552   2.026

d = stack(data, :SepalLength, :Species)
150x3 DataFrames.DataFrame
| Row | variable    | value | Species     |
|-----|-------------|-------|-------------|
| 1   | SepalLength | 5.1   | "setosa"    |
| 2   | SepalLength | 4.9   | "setosa"    |
| 3   | SepalLength | 4.7   | "setosa"    |
| 4   | SepalLength | 4.6   | "setosa"    |
| 5   | SepalLength | 5.0   | "setosa"    |
| 6   | SepalLength | 5.4   | "setosa"    |
| 7   | SepalLength | 4.6   | "setosa"    |
| 8   | SepalLength | 5.0   | "setosa"    |
| 9   | SepalLength | 4.4   | "setosa"    |
| 10  | SepalLength | 4.9   | "setosa"    |
| 11  | SepalLength | 5.4   | "setosa"    |
⋮
| 139 | SepalLength | 6.0   | "virginica" |

sort!(data, cols = [:SepalWidth])
150x5 DataFrames.DataFrame
| Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species      |
|-----|-------------|------------|-------------|------------|--------------|
| 1   | 5.0         | 2.0        | 3.5         | 1.0        | "versicolor" |
| 2   | 6.0         | 2.2        | 5.0         | 1.5        | "virginica"  |
| 3   | 6.0         | 2.2        | 4.0         | 1.0        | "versicolor" |
| 4   | 6.2         | 2.2        | 4.5         | 1.5        | "versicolor" |
| 5   | 4.5         | 2.3        | 1.3         | 0.3        | "setosa"     |
| 6   | 5.0         | 2.3        | 3.3         | 1.0        | "versicolor" |
| 7   | 5.5         | 2.3        | 4.0         | 1.3        | "versicolor" |
| 8   | 6.3         | 2.3        | 4.4         | 1.3        | "versicolor" |

sort!(data)
150x5 DataFrames.DataFrame
| Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species     |
|-----|-------------|------------|-------------|------------|-------------|
| 1   | 4.3         | 3.0        | 1.1         | 0.1        | "setosa"    |
| 2   | 4.4         | 2.9        | 1.4         | 0.2        | "setosa"    |
| 3   | 4.4         | 3.0        | 1.3         | 0.2        | "setosa"    |
| 4   | 4.4         | 3.2        | 1.3         | 0.2        | "setosa"    |
| 5   | 4.5         | 2.3        | 1.3         | 0.3        | "setosa"    |
| 6   | 4.6         | 3.1        | 1.5         | 0.2        | "setosa"    |
| 7   | 4.6         | 3.2        | 1.4         | 0.2        | "setosa"    |
| 8   | 4.6         | 3.4        | 1.4         | 0.3        | "setosa"    |
| 9   | 4.6         | 3.6        | 1.0         | 0.2        | "setosa"    |
| 10  | 4.7         | 3.2        | 1.3         | 0.2        | "setosa"    |
| 11  | 4.7         | 3.2        | 1.6         | 0.2        | "setosa"    |
⋮
| 139 | 7.1         | 3.0        | 5.9         | 2.1        | "virginica" |
| 140 | 7.2         | 3.0        | 5.8         | 1.6        | "virginica" |
| 141 | 7.2         | 3.2        | 6.0         | 1.8        | "virginica" |
| 142 | 7.2         | 3.6        | 6.1         | 2.5        | "virginica" |

sort!(data, rev=true)
150x5 DataFrames.DataFrame
| Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species     |
|-----|-------------|------------|-------------|------------|-------------|
| 1   | 7.9         | 3.8        | 6.4         | 2.0        | "virginica" |
| 2   | 7.7         | 3.8        | 6.7         | 2.2        | "virginica" |
| 3   | 7.7         | 3.0        | 6.1         | 2.3        | "virginica" |
| 4   | 7.7         | 2.8        | 6.7         | 2.0        | "virginica" |
| 5   | 7.7         | 2.6        | 6.9         | 2.3        | "virginica" |
| 6   | 7.6         | 3.0        | 6.6         | 2.1        | "virginica" |
| 7   | 7.4         | 2.8        | 6.1         | 1.9        | "virginica" |
| 8   | 7.3         | 2.9        | 6.3         | 1.8        | "virginica" |
| 9   | 7.2         | 3.6        | 6.1         | 2.5        | "virginica" |
| 10  | 7.2         | 3.2        | 6.0         | 1.8        | "virginica" |
| 11  | 7.2         | 3.0        | 5.8         | 1.6        | "virginica" |
⋮
| 139 | 4.8         | 3.0        | 1.4         | 0.1        | "setosa"    |
| 140 | 4.7         | 3.2        | 1.6         | 0.2        | "setosa"    |
| 141 | 4.7         | 3.2        | 1.3         | 0.2        | "setosa"    |

sort!(data, cols = [:Species, :SepalLength, :SepalWidth],
            rev = [true, false, false])
150x5 DataFrames.DataFrame
| Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species     |
|-----|-------------|------------|-------------|------------|-------------|
| 1   | 4.9         | 2.5        | 4.5         | 1.7        | "virginica" |
| 2   | 5.6         | 2.8        | 4.9         | 2.0        | "virginica" |
| 3   | 5.7         | 2.5        | 5.0         | 2.0        | "virginica" |
| 4   | 5.8         | 2.7        | 5.1         | 1.9        | "virginica" |
| 5   | 5.8         | 2.7        | 5.1         | 1.9        | "virginica" |
| 6   | 5.8         | 2.8        | 5.1         | 2.4        | "virginica" |
| 7   | 5.9         | 3.0        | 5.1         | 1.8        | "virginica" |

max_stars_repo_path: Section 3/3_3.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 3/3_3.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 3/3_3.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
####
#### Numerical methods for simulating Hamiltonian trajectory.
####

# TODO: The type `<:Tuple{Integer,Bool}` is introduced to address
# https://github.com/TuringLang/Turing.jl/pull/941#issuecomment-549191813
# We might want to simplify it to `Tuple{Int,Bool}` when we figured out
# why the it behaves unexpected on Windos 32.

abstract type AbstractIntegrator end

NesterovDualAveraging(δ::AbstractFloat, i::AbstractIntegrator) = NesterovDualAveraging(δ, nom_step_size(i))

stat(::AbstractIntegrator) = NamedTuple()

"""
    nom_step_size(::AbstractIntegrator)

Get the nominal integration step size. The current integration step size may
differ from this, for example if the step size is jittered. Nominal step size is
usually used in adaptation.
"""
nom_step_size(i::AbstractIntegrator) = step_size(i)

"""
    step_size(::AbstractIntegrator)

Get the current integration step size.
"""
function step_size end

abstract type AbstractLeapfrog{T} <: AbstractIntegrator end

step_size(lf::AbstractLeapfrog) = lf.ϵ
jitter(::Union{AbstractRNG,AbstractVector{<:AbstractRNG}}, lf::AbstractLeapfrog) = lf
temper(lf::AbstractLeapfrog, r, ::NamedTuple{(:i, :is_half),<:Tuple{Integer,Bool}}, ::Int) = r
stat(lf::AbstractLeapfrog) = (step_size=step_size(lf), nom_step_size=nom_step_size(lf))

function step(
    lf::AbstractLeapfrog{T},
    h::Hamiltonian,
    z::P,
    n_steps::Int=1;
    fwd::Bool=n_steps > 0,  # simulate hamiltonian backward when n_steps < 0
    res::Union{Vector{P}, P}=z
) where {T<:AbstractScalarOrVec{<:AbstractFloat}, P<:PhasePoint}
    n_steps = abs(n_steps)  # to support `n_steps < 0` cases

    ϵ = fwd ? step_size(lf) : -step_size(lf)
    ϵ = ϵ'

    @unpack θ, r = z
    @unpack value, gradient = ∂H∂θ(h, θ)
    for i = 1:n_steps
        # Tempering
        r = temper(lf, r, (i=i, is_half=true), n_steps)
        # Take a half leapfrog step for momentum variable
        r = r - ϵ / 2 .* gradient
        # Take a full leapfrog step for position variable
        ∇r = ∂H∂r(h, r)
        θ = θ + ϵ .* ∇r
        # Take a half leapfrog step for momentum variable
        @unpack value, gradient = ∂H∂θ(h, θ)
        r = r - ϵ / 2 .* gradient
        # Tempering
        r = temper(lf, r, (i=i, is_half=false), n_steps)
        # Create a new phase point by caching the logdensity and gradient
        z = phasepoint(h, θ, r; ℓπ=DualValue(value, gradient))
        !isfinite(z) && break
        if res isa Vector
            res[i] = z
        else
            res = z
        end
    end
    res
end

struct Leapfrog{T<:AbstractScalarOrVec{<:AbstractFloat}} <: AbstractLeapfrog{T}
    ϵ       ::  T
end
Base.show(io::IO, l::Leapfrog) = print(io, "Leapfrog(ϵ=$(round(l.ϵ; sigdigits=3)))")

### Jittering

struct JitteredLeapfrog{FT<:AbstractFloat,T<:AbstractScalarOrVec{FT}} <: AbstractLeapfrog{T}
    ϵ0      ::  T  # nominal (unjittered) step size
    jitter  ::  FT
    ϵ       ::  T  # current (jittered) step size
end

JitteredLeapfrog(ϵ0, jitter) = JitteredLeapfrog(ϵ0, jitter, ϵ0)

function Base.show(io::IO, l::JitteredLeapfrog)
    print(io, "JitteredLeapfrog(ϵ0=$(round(l.ϵ0; sigdigits=3)), jitter=$(round(l.jitter; sigdigits=3)), ϵ=$(round(l.ϵ; sigdigits=3)))")
end

nom_step_size(lf::JitteredLeapfrog) = lf.ϵ0

# Jitter step size; ref: https://github.com/stan-dev/stan/blob/1bb054027b01326e66ec610e95ef9b2a60aa6bec/src/stan/mcmc/hmc/base_hmc.hpp#L177-L178
function _jitter(
    rng::Union{AbstractRNG,AbstractVector{<:AbstractRNG}},
    lf::JitteredLeapfrog{FT,T}
) where {FT<:AbstractFloat,T<:AbstractScalarOrVec{FT}}
    ϵ = lf.ϵ0 .* (1 .+ lf.jitter .* (2 .* rand(rng) .- 1))
    return reconstruct(lf, ϵ=ϵ)
end

jitter(rng::AbstractRNG, lf::JitteredLeapfrog) = _jitter(rng, lf)

jitter(
    rng::AbstractVector{<:AbstractRNG},
    lf::JitteredLeapfrog{FT,T}
) where {FT<:AbstractFloat,T<:AbstractScalarOrVec{FT}} = _jitter(rng, lf)

### Tempering

struct TemperedLeapfrog{FT<:AbstractFloat,T<:AbstractScalarOrVec{FT}} <: AbstractLeapfrog{T}
    ϵ       ::  T
    α       ::  FT
end

function Base.show(io::IO, l::TemperedLeapfrog)
    print(io, "TemperedLeapfrog(ϵ=$(round(l.ϵ; sigdigits=3)), α=$(round(l.α; sigdigits=3)))")
end

"""
    temper(lf::TemperedLeapfrog, r, step::NamedTuple{(:i, :is_half),<:Tuple{Integer,Bool}}, n_steps::Int)

Tempering step. `step` is a named tuple with
- `i` being the current leapfrog iteration and
- `is_half` indicating whether or not it's (the first) half momentum/tempering step
"""
function temper(lf::TemperedLeapfrog, r, step::NamedTuple{(:i, :is_half),<:Tuple{Integer,Bool}}, n_steps::Int)
    if step.i > n_steps
        throw(BoundsError("Current leapfrog iteration exceeds the total number of steps."))
    end
    i_temper = 2(step.i - 1) + 1 + !step.is_half    # counter for half temper steps
    return i_temper <= n_steps ? r * sqrt(lf.α) : r / sqrt(lf.α)
end


text:
module HydroPowerModels

using JuMP, PowerModels, SDDP
import Reexport

mutable struct HydroPowerModel
    policygraph::SDDP.PolicyGraph
    alldata::Array{Dict{Any,Any}}
    params::Dict
end

include("variable.jl")
include("constraint.jl")
include("utilities.jl")
include("IO.jl")
include("simulate.jl")
include("train.jl")
include("visualization/visualize_data.jl")
include("objective.jl")
include("build_model.jl")

export  hydrothermaloperation, create_param,
        set_active_demand!, flat_dict, signif_dict

Reexport.@reexport using PowerModels, SDDP

"""
    hydrothermaloperation(alldata::Array{Dict{Any,Any}}, params::Dict)

Create a hydrothermal power operation model containing the policygraph the system data and the planning parameters.

Required parameters are:
-   alldata is a vector of dicts with information of the problem's stages.
-   param is a dict containing solution parameters.
"""
function hydrothermaloperation( alldata::Array{Dict{Any,Any}}, params::Dict;
                                build_model::Function=HydroPowerModels.build_opf_powermodels,
                                build_graph::Function=HydroPowerModels.build_graph)
    # verbose
    if !params["verbose"]
        PowerModels.silence()
    end

    # if set silence the solver
    # related to https://github.com/JuliaOpt/JuMP.jl/pull/1921
    if !params["verbose"]
        try
            MOI.set(JuMP.backend(Model(params["optimizer"])), MOI.Silent(), true)
        catch
            @info "Silent() attribute not implemented by the optimizer."
        end
    end

    # Model Definition
    graph = build_graph(params)
    policygraph = SDDP.PolicyGraph(graph,
                    sense       = :Min,
                    optimizer   = params["optimizer"],
                    #optimizer_forward = params["optimizer_forward"],
                    #optimizer_backward = params["optimizer_backward"],
                    lower_bound = 0.0,
                    direct_mode = false
                                            ) do sp, t #, isforward

        # if set silence the solver
        # related to https://github.com/JuliaOpt/JuMP.jl/pull/1921
        if !params["verbose"]
            try
                MOI.set(JuMP.backend(sp), MOI.Silent(), true)
            catch
                #@info "Silent() attribute not implemented by the optimizer."
            end
        end
        # Extract current data
        data = alldata[min(t,size(alldata,1))]

        # gather useful information from data
        gatherusefulinfo!(data)

        # build eletric grid model using PowerModels
        pm = build_opf_powermodels(sp,data,params)

        #if isforward # NOT YET IMPLEMENTED
            #pm = PowerModels.build_model(data["powersystem"], params["model_constructor_grid_forward"],
                #params["post_method"], jump_model=sp, setting = params["setting"])
        #else
            #pm = PowerModels.build_model(data["powersystem"], params["model_constructor_grid_backward"],
            #    params["post_method"], jump_model=sp, setting = params["setting"])
        #end

        # create reference to variables
        createvarrefs!(sp,pm)

        # save AbstractPowerModel and Data
        sp.ext[:pm] = pm
        sp.ext[:data] = data

        # save lower_bound
        sp.ext[:lower_bound] = 0.0

        # resevoir variables
        variable_volume(sp, data)

        # outflow and spillage variables
        variable_outflow(sp, data)
        variable_spillage(sp, data)

        # hydro balance
        variable_inflow(sp, data)
        rainfall_noises(sp, data, params, cidx(t,data["hydro"]["size_inflow"][1]))
        constraint_hydro_balance(sp, data, params)

        # hydro_generation
        constraint_hydro_generation(sp, data, pm)

        # deficit
        variable_deficit(sp, data, pm)
        constraint_mod_deficit(sp, data, pm)

        # costs stage
        variable_cost(sp, data)
        add_gen_cost(sp, data)
        add_spill_cost(sp, data)
        add_deficit_cost(sp, data)

        # Stage objective
        set_objective(sp, data)

        # # variable primal start
        # JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), JuMP.all_variables(sp), NaN)
        JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), sp[:deficit], 0.0)
        JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), sp[:inflow], 0.0)
        JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), sp[:outflow], 0.0)
        JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), sp[:spill], 0.0)
        for r in sp[:reservoir]
            JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), r.in, 0.0)
        end
        for r in sp[:reservoir]
            JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), r.out, 0.0)
        end

    end

    # save data
    m = HydroPowerModel(policygraph,alldata,params)

    return m
end

end


text:
# TODO Cleanup
struct Variable
	label::String
	parents::Vector{String}
	cardinality::Integer
	sample_size::Integer
	states_possible::Integer
	states_missing::Integer
	states_parents::Integer
	cpt::DataFrame

	function Variable(label::String, dag::DAG, data::DataFrame)
		#: sample size
		sample_size = size(data)[1]
		parents = inneighbors(dag, label)
		#: conditional state counting
		cpt = combine(groupby(data, [parents; label]), nrow => :count)
		transform!(groupby(cpt, parents), :count => sum => :total)
		#: extracting marginal states
		states = [unique(cpt[!, name]) for name in [parents; label]]
		#: cardinality of the variable
		cardinality = size(last(states))[1]
		#: number of possible conditional states
		states_possible = prod(size(state)[1] for state in states)
		#: number missing states
		states_missing = states_possible - size(cpt)[1]
		#: number of parents states
		states_parents = states_possible/cardinality
		return new(label, parents, cardinality, sample_size, states_possible, states_missing, states_parents, cpt)
	end
end

gammaln(x) = logabsgamma(x)[1]

function score_bdeu(variables::Vector{Variable}, equivalent_sample_size=10)
	score::Float64 = 0.0
	for variable in variables
		alpha = equivalent_sample_size / variable.states_parents
		beta = equivalent_sample_size / variable.states_possible
		log_gamma_counts_beta = sum(gammaln.(variable.cpt.count .+ beta)) + variable.states_missing * gammaln(beta)
		state_counts_agg = groupby(variable.cpt, variable.parents)
		log_gamma_cond = sum(gammaln.([sum(x.count) + alpha for x in state_counts_agg])) + (variable.states_parents - length(state_counts_agg)) * gammaln(alpha)
		local_score::Float64 = log_gamma_counts_beta - log_gamma_cond + variable.states_parents * gammaln(alpha) - variable.states_possible * gammaln(beta)
		score += local_score
	end
	return score
end

function score_bic(variables::Vector{Variable})
	score::Float64 = 0.0
	for variable in variables
		log_likelihoods = log.(variable.cpt.count)
		log_conditionals = log.(variable.cpt.total)
		log_likelihoods -= log_conditionals
		log_likelihoods .*= variable.cpt.count
		local_score::Float64 = sum(log_likelihoods)
		local_score -= 0.5 * log(variable.sample_size) * variable.states_parents * (variable.cardinality - 1)
		score += local_score
	end
	return score
end


text:
using Plots
using Random
using VegaLite
using DataFrames
using LinearAlgebra
using TripletEmbeddings

Random.seed!(4)

# We define a random Embedding in ℜ²
d = 2
n = 100
Y = rand(d, n)

triplets = Triplets(Y)
loss = STE(σ = 1/sqrt(2))
X = Embedding(size(Y))

@time violations = fit!(loss, triplets, X; verbose=true, max_iterations=200)
procrustes!(X, Y)

scatter(Y[1,:], Y[2,:]); scatter!(X[1,:], X[2,:])

dfX = DataFrame(embedding = "X", time = 1:n, x = X[1,:], y = X[2,:])
dfY = DataFrame(embedding = "Y", time = 1:n, x = Y[1,:], y = Y[2,:])

vcat(dfX, dfY) |> @vlplot(:point, x = :x, y = :y, color=:embedding, width=600, height=400)


text:
mutable struct LinearRegression{T} <: AbstractMLModel
    β0::T
    β::Vector{T}
    A::Matrix{T}
end
LinearRegression(T = Float64) = LinearRegression(zero(T), zeros(T, 0), zeros(T, 0, 0))

algorithm(::LinearRegression{T}) where {T} = Sweep(T)

mutable struct Sweep{T} <: ClosedFormSolver
    A::Matrix{T}  # [x 1 y]' * [x 1 y] =
end
Sweep(T=Float64) = Sweep(zeros(T, 0, 0))

# function update!(model::LinearRegression{T}, strat::Sweep{T}, data::MLData{X,Y,Nothing,:rows}) where {T,X,Y}
#     x, y = xy(data)
#     n = nobs(data)
#     p = npredictors(data)
#     strat.A = zeros(T, p + 2, p + 2)
#     model.β = zeros(T, p)
#     strat.A[1:end-2, 1:end-2] = x'x / n
#     strat.A[1:end-2, end-1] = vec(mean(x, 1))
#     strat.A[end-1, ]
#     strat.A[end] = one(T)
# end

# function learn!(model, strat::LearningStrategy, data)
#     setup!(strat, model[, data])
#     for (i, item) in enumerate(data)
#         update!(model, strat[, i], item)
#         hook(strat, model[, data], i)
#         finished(strat, model[, data], i) && break
#     end
#     cleanup!(strat, model)
#     model
# end


text:
using LULESH
using MPI
using Enzyme
using LLVM
using Printf

LLVM.clopts("-memdep-block-scan-limit=70000")
LLVM.clopts("-dse-memoryssa-walklimit=10000")
LLVM.clopts("-attributor-max-iterations=128")
LLVM.clopts("-capture-tracking-max-uses-to-explore=256")

# Enzyme.API.printperf!(true)
# Enzyme.API.printall!(true)
# Enzyme.API.instname!(true)

Enzyme.API.inlineall!(true)
# Size of Domain is around 1024
Enzyme.API.maxtypeoffset!(32)
isdefined(Enzyme.API, :strictAliasing!) && Enzyme.API.strictAliasing!(true)
isdefined(Enzyme.API, :typeWarning!) &&  Enzyme.API.typeWarning!(false)
Enzyme.API.looseTypeAnalysis!(true)

function main(nx, structured, num_iters, mpi, enzyme, verification)
    # TODO: change default nr to 11
    nr = 1
    balance = 1
    cost = 1
    floattype = Float64
    devicetype = Vector

    if mpi
        !MPI.Initialized() && MPI.Init()
        comm = MPI.COMM_WORLD
    else
        comm = nothing
    end

    if getMyRank(comm) == 0
        @info "Constructing LuleshProblem" num_iters structured nx nr balance cost ranks=getNumRanks(comm)
    end

    prob = LuleshProblem(num_iters, structured, nx, nr, balance, cost, devicetype, floattype, comm)

    # Set up the mesh and decompose. Assumes regular cubes for now
    # TODO: modify this constructor to account for new fields
    # TODO: setup communication buffers

    domain = Domain(prob)
    if enzyme
        shadowDomain = Domain(prob)
        if verification
            @time Enzyme.autodiff(doubleFrog, Duplicated(domain, shadowDomain))
        else
            @time Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
        end
        domain = Domain(prob)
    elseif verification
        # not enzyme and verification means we need perturbed domain for finite differences
        perturbedDomain = Domain(prob)
    end

    # getnodalMass = nodalMass(domain)

    # Initial domain boundary communication
    commRecv(domain, MSG_COMM_SBN, 1,
                domain.sizeX + 1, domain.sizeY + 1, domain.sizeZ + 1,
                true, false)
    fields = (domain.nodalMass,)
    commSend(domain, MSG_COMM_SBN, fields,
             domain.sizeX + 1, domain.sizeY + 1, domain.sizeZ + 1,
             true, false)
    commSBN(domain, fields)

    # End initialization
    if mpi
        MPI.Barrier(prob.comm::MPI.Comm)
    end

    if getMyRank(prob.comm) == 0
        if (structured)
            @info "Running" until=domain.stoptime domain=(nx,nx,nx)
        else
            @info "Running" until=domain.stoptime domain=domain.numElem
            @warn "Unstructured setup not supported"
        end
    end

    eps = 1e-6
    if verification
        num_iters = 1
        if enzyme
            println("seeding gradients")
            for f in fieldnames(Domain)
                prop = getproperty(shadowDomain, f)
                if(typeof(prop) == Float64)
                    setproperty!(shadowDomain,f,0.0)
                    @printf("set %s to %f\n",f,sum(getproperty(shadowDomain,f)))
                elseif(typeof(prop) == Vector{Float64})
                    fill!(prop, 0.0)
                    @printf("filled %s with %f\n",f,sum(getproperty(shadowDomain,f)))
                end
            end
            fill!(shadowDomain.e, 1.0)
        else
            # finite difference initialization
            println("seeding FD")
            domain.e[1] += eps
            perturbedDomain.e[1] -= eps
        end
    end

    # timestep to solution
    start = getWtime(prob.comm)
    while domain.time < domain.stoptime
        # this has been moved after computation of volume forces to hide launch latencies
        timeIncrement!(domain)
        if enzyme
            if verification
                Enzyme.autodiff(doubleFrog, Duplicated(domain, shadowDomain))
            else
                Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
            end
        else
            if verification
                doubleFrog(domain)
                timeIncrement!(perturbedDomain)
                doubleFrog(perturbedDomain)
            else
                lagrangeLeapFrog(domain)
            end
        end

        # checkErrors(domain, its, myRank)
        if getMyRank(prob.comm) == 0
            @info "Completed" cycle=domain.cycle time=domain.time dt=domain.deltatime
        end
        if domain.cycle >= num_iters
            break
        end
    end

    # Use reduced max elapsed time
    elapsed_time = getWtime(prob.comm) - start
    elapsed_timeG = comm_max(elapsed_time, prob.comm)

    if verification
        if enzyme
            println("harvesting gradients")
            @printf("gradient checksum: %f\n", shadowDomain.e[1])
        else
            @printf("FD checksum: %f\n",sum(domain.e - perturbedDomain.e) / (2*eps))
        end
    end

    if getMyRank(prob.comm) == 0
        verifyAndWriteFinalOutput(elapsed_timeG, domain, nx, getNumRanks(prob.comm))
    end

    if mpi
        MPI.Finalize()
    end
end

if !isinteractive()
    args = LULESH.parse_cmd()

    # assume cube subdomain geometry for now (nx)
    nx = args["N"]
    structured = args["s"]
    num_iters = args["num_iters"]
    mpi = args["mpi"]
    enzyme = args["enzyme"]
    verification = args["verification"]
    main(nx, structured, num_iters, mpi, enzyme, verification)
end


text:
# Test the show functions on different kinds of graphs

using Graphs

adjlst = simple_adjlist(3)
add_edge!(adjlst, 1, 2)
add_edge!(adjlst, 1, 3)
add_edge!(adjlst, 2, 3)

println("simple_adjlist")
println("====================")
println(adjlst)
println()

inclst = simple_inclist(4, is_directed=false)
add_edge!(inclst, 1, 2)
add_edge!(inclst, 1, 3)
add_edge!(inclst, 2, 4)
add_edge!(inclst, 3, 4)

println("simple_inclist")
println("====================")
println(inclst)
println()

sg = simple_graph(4)
add_edge!(sg, 1, 2)
add_edge!(sg, 1, 3)
add_edge!(sg, 2, 4)
add_edge!(sg, 3, 4)

println("simple_graph")
println("====================")
println(sg)
println()

eg = graph(ExVertex, ExEdge{ExVertex})
v1 = add_vertex!(eg, "a")
v2 = add_vertex!(eg, "b")
v3 = add_vertex!(eg, "c")
add_edge!(eg, v1, v2)
add_edge!(eg, v1, v3)
add_edge!(eg, v2, v3)

println("extended_graph")
println("====================")
println(eg)
println()


text:
using LinearAlgebra

export blocknewton

default_block_errmeasure(nep::NEP) = (S,X) -> opnorm(compute_MM(nep,S,X))

"""
    (S,X)=blocknewton(nep [S,] [X,] [errmeasure,] [tol,] [maxit,] [armijo_factor,] [armijo_max,] [logger])

Applies the block Newton method to `nep::AbstractSPMF`. The method computes
an invariant pair `(S,X)` using the block Newton approach of Kressner.
The variables `S`,`X`
correspond to starting approximations.
The function `errmeasure` shoule be defined for errmeasure(S,X)
and meausures the error in the pair `(S,X)`.
See [`augnewton`](@ref) for other parameters.

# Example
The example shows that `compute_MM()` becomes zero when a solution
has been computed.
```julia-repl
julia> nep=nep_gallery("dep0",3);
julia> (S,X)= blocknewton(nep)
julia> compute_MM(nep,S,X)
3×2 Array{Complex{Float64},2}:
 -1.11022e-16+0.0im  1.11022e-16+0.0im
          0.0+0.0im          0.0+0.0im
  1.38778e-17+0.0im  2.77556e-17+0.0im
```
This example solves the `gun` problem from the Berlin-Manchester collection
```julia-repl
julia> using NonlinearEigenproblems.Gallery
julia> nep=nep_gallery("nlevp_native_gun");
julia> II=[1.0 0; 0 1]; S=150^2*II; V=[II;zeros(size(nep,1)-2,2)];
julia> (Z,X)=blocknewton(nep,S=S,X=V,logger=1,armijo_factor=0.5,maxit=20)
Iteration 1: Error: 6.081316e+03
Iteration 2: Error: 1.701970e-02 Armijo scaling=0.031250
Iteration 3: Error: 1.814887e-02 Armijo scaling=0.250000
...
Iteration 13: Error: 6.257442e-09
Iteration 14: Error: 2.525942e-15
```

# References
* D. Kressner A block Newton method for nonlinear eigenvalue problems, Numer. Math., 114 (2) (2009), pp. 355-372
"""
function blocknewton(nep::AbstractSPMF;
                     S::Matrix=zeros(2,2),
                     X::Matrix=Matrix(1.0I, size(nep,1), 2),
                     errmeasure::Function =  default_block_errmeasure(nep::NEP),
                     tol::Real=eps(real(eltype(S)))*100,
                     maxit::Integer=10,
                     logger=0,
                     armijo_factor::Real=1,
                     armijo_max::Integer=5)

    @parse_logger_param!(logger)

    T=complex(eltype(S))
    # This implementation is for complex arithmetic only
    # since the original paper is based on the Schur form (not real Schur form)
    # This could potentially be relaxed to real problems if the
    # eigenvalue approximations remain real.
    S=complex(S);
    X=complex(X);

    n=size(nep,1);
    p=size(S,1);

    # Construct the orthogonalization matrix
    W=Vl(X,S);
    # reshape W to WW (3D matrix)
    WW=zeros(T,n,p,p)
    for j=1:size(W,2)
        WW[:,:,j] = W[(j-1)*n .+ (1:n), :]
    end

    local err0=0
    # Main loop
    for k=1:maxit
        err0=errmeasure(S,X)
        push_iteration_info!(logger,k,err=err0,continues=true);
        if (err0<tol)
            push_iteration_info!(logger,"");
            return S,X
        end

        Res= compute_MM(nep,S,X)
        # Solve the linear system by first transforming S to
        # upper triangular form and then doing a backward substitution
        # and finally reversing the backward substitution
        RR,QQ=schur(complex(S))
        dSt,dXt = newtonstep_linsys(T,nep,
                                    RR, X*QQ,
                                    WW,  # Orthogonalization matrix
                                    Res*QQ, # RT
                                    zeros(T,p,p),  #RV=0
                                    logger);
        dX=dXt*QQ';
        dS=QQ*dSt*QQ';

        local j=1
        # Update the approximations
        if (armijo_factor<1)
            (ΔS,ΔV,j,scaling)=  armijo_rule_block(nep,errmeasure,
                                                  err0,
                                                  S,X,-dS,-dX,
                                                  armijo_factor,armijo_max)

            St=S+ΔS
            Xt=X+ΔV
        else
            j = 0
            St=S-dS
            Xt=X-dX

        end
        if (j>0)
            push_info!(logger," Armijo scaling=$scaling")
        else
            push_info!(logger,"");
        end

        # Carry out the orthogonalization
        W,R = qr(Vl(Xt,St))
        W = Matrix(W)
        # reshape W to WW (3D matrix)
        for j=1:size(W,2)
            WW[:,:,j] = W[(j-1)*n .+ (1:n), :]
        end
        X[:]=Xt/R; S[:]=(R*St)/R;
    end
    msg="Number of iterations exceeded. maxit=$(maxit)."
    throw(NoConvergenceException(S,X,err0,msg))

end

# Solve one linear system for a triangular S
# (S,X) Schur pair approx
# W tensor with orthogonalization coeffs
# RT,RV the right-hand side of linear system
function newtonstep_linsys(::Type{T},nep::AbstractSPMF,S, X, W, RT, RV, logger) where {T}

    n = size(nep,1);
    p = size(X,2); l = size(W,3);
    dX = zeros(T, n,p); dS = zeros(T,p,p);

    # Fetch the SPMF
    fv=get_fv(nep);
    Av=get_Av(nep)

    m = size(fv,1);

    # Initialize a tensor with a f_i(S)
    fS = zeros(T,p,p,m);
    for j = 1:m
        fS[:,:,j] = fv[j](S);
    end

    # Work column by-column
    for i = 1:p,
        s = S[i,i];

        # Construct the matrix in equation (20) in Kressner Num. Math.
        T11=compute_Mder(nep,s)
        S_expanded = [S I; zero(S) s*I]
        # T12 = compute_MM(nep,[0*X X],S_expanded) # Can maybe be computed like this? Would avoid the explicit use of Av and and fv
        T12 = zeros(T,n,p);
        for j = 1:m
            DF = fv[j](S_expanded)
            DF1=DF[1:p,p+1:2*p];
            T12 = T12 + Av[j]*X*DF1;
        end
        T21 = W[:,:,1]';
        for j = 2:l
            T21 = T21 + s^(j-1) * W[:,:,j]';
        end
        DS = Matrix(1.0I, p, p)
        T22 = zeros(T,p,p);
        for j = 2:l
            T22 = T22 + W[:,:,j]'*X*DS;
            DS = s*DS + S^(j-2)
        end
        TT=[T11 T12; T21 T22]; # This can maybe be optimized with Schur complement
        sol =  TT \ [RT[:,i];RV[:,i]];
        # compute the dS and dX
        dX[:,i] = sol[1:n];
        dS[:,i] = sol[n+1:end];

        # Update RHS (equation (21)-(22) in Kressner Nummath)
        if (i<p) # Updated RHS needed in the next for loop
            Z = zeros(T,p,p);
            Z[:,i] = dS[:,i]; DS = Z;
            S2_expanded = [S Z; zero(S) S]
            for j = 1:m  # Update\tilde{R}_{T2} according to equation (21)
                Za=dX[:,i]*transpose(fS[i,i+1:p,j]) # First term
                DF = fv[j](S2_expanded)
                Zb=X*DF[1:p,p+i+1:2*p]; # Second term
                RT[:,i+1:p] +=  - Av[j] * (Za+Zb);
            end
            for j = 2:l # Update \tilde{R}_{V2} according to equation (22)
                Za=dX[:,i] * transpose((S^(j-2))[i,i+1:p])
                Zb=X*DS[:,i+1:p] ;
                RV[:,i+1:p] +=  -W[:,:,j]' * ( Za +Zb );
                DS = DS*S + S^(j-2)*DS;
            end
        end
    end
    return dS,dX
end

# Construct the orthogonalization matrix
# V=[X;X*S;XS^2;...]
function Vl(X,S)
    p=size(S,1); n=size(X,1);
    V=zeros(eltype(S),n*p,p);
    for j=1:p
        V[(j-1)*n .+ (1:n), :] = X * (S^(j-1))
    end
    return V
end

function armijo_rule_block(nep,errmeasure,err0,S,V,ΔS,ΔV,armijo_factor,armijo_max)
    j=0
    if (armijo_factor<1)
        # take smaller and smaller steps until errmeasure is decreasing
        while (errmeasure(S+ΔS,V+ΔV)>err0 && j<armijo_max)
            j=j+1;
            ΔS=ΔS*armijo_factor;
            ΔV=ΔV*armijo_factor;
        end
    end
    return  (ΔS,ΔV,j,armijo_factor^j)
end


text:
const PETLION_VERSION = (Meta.parse.(split(String(Symbol(@PkgVersion.Version)),"."))...,)
const options = Dict{Symbol,Any}(
    :SAVE_SYMBOLIC_FUNCTIONS => true,
    :FILE_DIRECTORY => nothing,
    :FACTORIZATION_METHOD => :KLU, # :KLU or :LU
)

# https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units
const const_Faradays  = 96485.3321233
const const_Ideal_Gas = 8.31446261815324

abstract type AbstractJacobian <: Function end
abstract type AbstractModel{T<:AbstractJacobian,temp,solid_diff,Fickian,age} end
abstract type AbstractMethod end
abstract type AbstractRun{T<:AbstractMethod,input<:Any} end

struct method_I   <: AbstractMethod end
struct method_V   <: AbstractMethod end
struct method_P   <: AbstractMethod end

struct discretizations_per_section
    p::Int64
    s::Int64
    n::Int64
    a::Int64
    z::Int64
    r_p::Int64
    r_n::Int64
    diff::Int64
    alg::Int64
    tot::Int64
end
discretizations_per_section(p,s,n,a,z,r_p,r_n) = discretizations_per_section(p,s,n,a,z,r_p,r_n,-1,-1,-1)

Base.@kwdef mutable struct run_info
    exit_reason::String = ""
    flag::Int64 = -1
    iterations::Int64 = 1
end

struct run_constant{T<:AbstractMethod,in<:Union{Number,Symbol,Function}} <: AbstractRun{T,in}
    input::in
    value::Vector{Float64}
    method::T
    t0::Float64
    tf::Float64
    name::String
    info::run_info
end
struct run_function{T<:AbstractMethod,func<:Function} <: AbstractRun{T,func}
    func::func
    value::Vector{Float64}
    method::T
    t0::Float64
    tf::Float64
    name::String
    info::run_info
end
@inline value(run::AbstractRun) = @inbounds run.value[1]

Base.@kwdef struct index_state <: AbstractUnitRange{Int64}
    start::Int64 = 0
    stop::Int64 = 0
    a::UnitRange{Int64} = 0:0
    p::UnitRange{Int64} = 0:0
    s::UnitRange{Int64} = 0:0
    n::UnitRange{Int64} = 0:0
    z::UnitRange{Int64} = 0:0
    sections::Tuple = ()
    var_type::Symbol = :NA
end

struct jacobian_symbolic{T<:Function} <: AbstractJacobian
    func::T
    sp::SparseMatrixCSC{Float64,Int64}
end
(jac::jacobian_symbolic{<:Function})(x...) = jac.func(x...)

struct res_FD{T<:Function} <: Function
    f!::T
    Y_cache::Vector{Float64}
    YP_cache::Vector{Float64}
    θ_tot::Vector{Float64}
    N::discretizations_per_section
end
function (res_FD::res_FD{T})(res::AbstractVector{<:Number}, Y::AbstractVector{<:Number}) where T<:Function
    if length(Y) == res_FD.N.tot
        Y_new = Y
    else
        Y_new = zeros(eltype(Y), length(res_FD.Y_cache))
        @inbounds @views Y_new[1:res_FD.N.diff] .= res_FD.Y_cache[1:res_FD.N.diff]
        @inbounds @views Y_new[res_FD.N.diff+1:end] .= Y
    end
    res_FD.f!(res, 0.0, Y_new, res_FD.YP_cache, res_FD.θ_tot)
end

struct jacobian_AD{T<:Function} <: AbstractJacobian
    f!::res_FD{T}
    sp::SparseMatrixCSC{Float64,Int64}
    jac_cache::ForwardColorJacCache
end
@inline function (jac::jacobian_AD{T})(t,Y,YP,γ::Float64,p::P,run) where {T<:Function,P<:AbstractModel}
    J = jac.sp
    forwarddiff_color_jacobian!(J, jac.f!, Y, jac.jac_cache)
    if size(J) == (p.N.tot-1,p.N.tot)
        @inbounds for i in 1:p.N.diff
            J[i,i] += -γ
        end
    end
    return nothing
end

struct jacobian_combined{
    T1<:Function,
    T2<:Union{SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false},SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}},
    T3<:Function,
    T4<:LinearAlgebra.Factorization{Float64},
    }
    sp::SparseMatrixCSC{Float64,Int64}
    base_func::T1
    J_base::T2
    scalar_func::T3
    J_scalar::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}
    θ_tot::Vector{Float64}
    θ_keys::Vector{Symbol}
    factor::T4
end
Base.getindex(J::jacobian_combined,i...) = getindex(J.sp,i...)
Base.axes(J::jacobian_combined,i...) = axes(J.sp,i...)

struct residual_combined{
    T1<:Function,
    T2<:Function,
    T3<:Function,
    }
    f_diff!::T1
    f_alg!::T2
    f_scalar!::T3
    ind_diff::UnitRange{Int64}
    ind_alg::UnitRange{Int64}
    θ_tot::Vector{Float64}
    θ_keys::Vector{Symbol}
end
struct Jac_and_res{T<:Sundials.IDAIntegrator}
    J_full::jacobian_combined
    R_full::residual_combined
    J_alg::jacobian_combined
    R_diff::residual_combined
    R_alg::residual_combined
    int::Vector{T}
end

Base.@kwdef mutable struct boundary_stop_prev_values
    t_final_interp_frac::Float64 = +1.0
    V::Float64 = -1.0
    SOC::Float64 = -1.0
    T::Float64 = -1.0
    c_s_n::Float64 = -1.0
    I::Float64 = -1.0
    η_plating::Float64 = -1.0
    c_e_min::Float64 = -1.0
    dfilm::Float64 = -1.0
end

remove_module_name(x::String) = replace(x, "$(@__MODULE__)."=>"")
function create_immutable_version(structure::DataType; str_replacements=(""=>"",), conv_replacements=(""=>"",))
    """
    Immutable versions of struct can be beneficial for performance. This function
    creates an immutable version of mutable structures with "_immutable" appended to the end.
    `str_replacements` uses the `replace` function to substitute string values in the evaluated
    struct.
    """

    name = String(Symbol(structure))
    fields = String.(fieldnames(structure))
    types = String.(Symbol.(fieldtypes(structure)))

    args = join([field * "::" * type for (field,type) in zip(fields,types)], "\n")

    super = supertype(structure)

    # string to create new struct
    str_immutable = "struct $(name)_immutable <: $super
    $(args)
    end"

    if str_replacements isa Pair{String, String}
        str_replacements = (str_replacements,)
    end
    if conv_replacements isa Pair{String, String}
        conv_replacements = (conv_replacements,)
    end

    # string to convert the mutable struct to the immutable struct
    conversion = "$(name)_immutable(x::$(name)) = $(name)_immutable($(join(map(field->"x.$field", fields), ",")))"

    for replacement in str_replacements
        str_immutable = replace(str_immutable, replacement)
    end
    for replacement in conv_replacements
        conversion = replace(conversion, replacement)
    end

    # The names cannot be prepended by "PETLION.", otherwise it creates errors
    str_immutable = remove_module_name(str_immutable)
    conversion = remove_module_name(conversion)

    eval.(Meta.parse.((str_immutable,conversion)))

    return nothing
end

abstract type AbstractStopConditions end
Base.@kwdef mutable struct boundary_stop_conditions <: AbstractStopConditions
    V_max::Float64 = -1.0
    V_min::Float64 = -1.0
    SOC_max::Float64 = -1.0
    SOC_min::Float64 = -1.0
    T_max::Float64 = -1.0
    c_s_n_max::Float64 = -1.0
    I_max::Float64 = NaN
    I_min::Float64 = NaN
    η_plating_min::Float64 = NaN
    c_e_min::Float64 = NaN
    dfilm_max::Float64 = NaN
    prev::boundary_stop_prev_values = boundary_stop_prev_values()
end
create_immutable_version(boundary_stop_conditions; conv_replacements=", x.prev" => ", x.$(fieldtypes(boundary_stop_conditions)[end])()")

Base.@kwdef mutable struct _funcs_numerical
    rxn_p::Function = emptyfunc
    rxn_n::Function = emptyfunc
    OCV_p::Function = emptyfunc
    OCV_n::Function = emptyfunc
    D_s_eff::Function = emptyfunc
    D_eff::Function = emptyfunc
    K_eff::Function = emptyfunc
    thermodynamic_factor::Function = emptyfunc
end

struct options_numerical{temp,solid_diff,Fickian,age}
    temperature::Bool
    solid_diffusion::Symbol
    Fickian_method::Symbol
    aging::Union{Bool,Symbol}
    cathode::Function
    anode::Function
    rxn_p::Function
    rxn_n::Function
    OCV_p::Function
    OCV_n::Function
    D_s_eff::Function
    rxn_rate::Function
    D_eff::Function
    K_eff::Function
    thermodynamic_factor::Function
    jacobian::Symbol
end
options_numerical(temp,solid_diff,Fickian,age,x...) =
    options_numerical{temp,solid_diff,Fickian,age}(temp,solid_diff,Fickian,age,x...)

const states_logic = model_states{
    Bool,
    Bool,
    <:Tuple
}

const indices_states = model_states{
    index_state,
    index_state,
    Tuple,
}

abstract type AbstractOptionsModel end
Base.@kwdef mutable struct options_simulation <: AbstractOptionsModel
    outputs = (:t, :V)
    SOC::Number = 1.0
    abstol::Float64 = 1e-6
    reltol::Float64 = 1e-6
    abstol_init::Float64 = abstol
    reltol_init::Float64 = reltol
    maxiters::Int64 = 10_000
    check_bounds::Bool = true
    reinit::Bool = true
    verbose::Bool = false
    interp_final::Bool = true
    tstops::Vector{<:Number} = Float64[]
    tdiscon::Vector{<:Number} = Float64[]
    interp_bc::Symbol = :interpolate
    save_start::Bool = false
    var_keep::states_logic = model_states_logic(outputs)[1]
    stop_function::Function = (x...) -> nothing
    calc_integrator::Bool = false
end

create_immutable_version(options_simulation; str_replacements=(
    "_immutable" =>"_immutable{T<:Function}",
    "stop_function::Function" => "stop_function::T",
    "outputs::Any" => "outputs::Tuple",
    ))

struct save_start_info{T<:AbstractMethod}
    method::T
    SOC::Float64 # rounded to 3rd decimal place
    I::Float64 # rounded to 3rd decimal place
end

struct cache_run
    θ_tot::Vector{Float64}
    θ_keys::Vector{Symbol}
    state_labels::Vector{Symbol}
    vars::Tuple
    save_start_dict::Dict{save_start_info,Vector{Float64}}
    Y0::Vector{Float64}
    YP0::Vector{Float64}
    res::Vector{Float64}
    Y_alg::Vector{Float64}
    id::Vector{Int64}
end

struct model_funcs{T1<:Function,T2<:Function,T3<:Function,T4<:AbstractJacobian,T5<:AbstractJacobian}
    initial_guess!::T1
    f_diff!::T2
    f_alg!::T3
    J_y!::T4
    J_y_alg!::T5
    Dict_constant::Dict{DataType,Jac_and_res}
    Dict_function::Dict{DataType,Dict{DataType,Jac_and_res}}
    Dict_residual::Dict{DataType,Jac_and_res}
end
model_funcs(x...) = model_funcs(x...,
    Dict{DataType,Jac_and_res}(),
    Dict{DataType,Dict{DataType,Jac_and_res}}(),
    Dict{DataType,Jac_and_res}()
    )

struct model{T<:AbstractJacobian,temp,solid_diff,Fickian,age} <: AbstractModel{T,temp,solid_diff,Fickian,age}
    θ::Dict{Symbol,Float64}
    numerics::options_numerical{temp,solid_diff,Fickian,age}
    N::discretizations_per_section
    ind::indices_states
    opts::options_simulation
    bounds::boundary_stop_conditions
    cache::cache_run
    funcs::model_funcs{<:Function,<:Function,<:Function,T,<:AbstractJacobian}
end

const model_jac{jac}               = model{jac,<:Any,<:Any,<:Any,<:Any}
const model_temp{temp}             = model{<:AbstractJacobian,temp,<:Any,<:Any,<:Any}
const model_solid_diff{solid_diff} = model{<:AbstractJacobian,<:Any,solid_diff,<:Any,<:Any}
const model_Fickian{Fickian}       = model{<:AbstractJacobian,<:Any,<:Any,Fickian,<:Any}
const model_age{age}               = model{<:AbstractJacobian,<:Any,<:Any,<:Any,age}

const AbstractModelJac{jac}              = AbstractModel{jac,<:Any,<:Any,<:Any,<:Any}
const AbstractModelTemp{temp}            = AbstractModel{<:AbstractJacobian,temp,<:Any,<:Any,<:Any}
const AbstractModelSolidDiff{solid_diff} = AbstractModel{<:AbstractJacobian,<:Any,solid_diff,<:Any,<:Any}
const AbstractModelFickian{Fickian}      = AbstractModel{<:AbstractJacobian,<:Any,<:Any,Fickian,<:Any}
const AbstractModelAge{age}              = AbstractModel{<:AbstractJacobian,<:Any,<:Any,<:Any,age}

struct model_skeleton{temp,solid_diff,Fickian,age} <: AbstractModel{AbstractJacobian,temp,solid_diff,Fickian,age}
    θ::Dict{Symbol,Any}
    numerics::options_numerical{temp,solid_diff,Fickian,age}
    N::discretizations_per_section
    ind::indices_states
    opts::options_simulation
    bounds::boundary_stop_conditions
    cache::cache_run
end

struct run_results{T<:AbstractRun}
    run::T
    tspan::Tuple{Float64,Float64}
    info::run_info
    run_index::UnitRange{Int64}
    int::Sundials.IDAIntegrator
    opts::options_simulation_immutable
    bounds::boundary_stop_conditions_immutable
    N::discretizations_per_section
    numerics::options_numerical
    p::model
end

const solution = model_states{
    Array{Float64,1},
    VectorOfArray{Float64,2,Array{Array{Float64,1},1}},
    Array{run_results,1},
}

Base.length(sol::solution) = length(sol.results)
Base.isempty(sol::solution) = isempty(sol.results)
function Base.getindex(sol::T, i1::Int) where T<:solution
    ind = (sol.results[i1].run_index) .+ (1-sol.results[1].run_index[1])
    T([fields == :results ? [sol.results[i1]] : (x = getproperty(sol, fields); length(x) > 1 ? x[ind] : x) for fields in fieldnames(T)]...)
end
function Base.getindex(sol::T, i::UnitRange{Int64}) where T<:solution
    ind = ((sol.results[i[1]].run_index[1]):(sol.results[i[end]].run_index[end])) .+ (1-sol.results[1].run_index[1])
    T([fields == :results ? sol.results[i] : (x = getproperty(sol, fields); length(x) > 1 ? x[ind] : x) for fields in fieldnames(T)]...)
end
Base.lastindex(sol::T) where T<:solution = length(sol)
Base.firstindex(::T) where T<:solution = 1

Base.empty!(f::model_funcs) = ([empty!(getproperty(f,field)) for (_type,field) in zip(fieldtypes(model_funcs),fieldnames(model_funcs)) if _type <: Dict];nothing)
Base.empty!(p::model) = empty!(p.funcs)

const STATE_NAMES = Dict{Symbol,String}(
    :c_e => "Electrolyte Conc. (mol/m³)",
    :c_s_avg => "Solid Conc. (mol/m³)",
    :T => "Temperature (K)",
    :film => "Li Plating Thickness (m)",
    :Q => "Q",
    :j => "Molar Ionic Flux (mol/m²⋅s)",
    :j_s => "Side Reaction Flux (mol/m²⋅s)",
    :Φ_e => "Electrolyte Potential (V)",
    :Φ_s => "Solid-phase Potential (V)",
    :I => "Current (C-rate)",
    :V => "Voltage (V)",
    :P => "Power (W)",
    :SOC => "State-of-Charge (-)",
    :SOH => "State-of-Health (-)",
)

## Modifying Base functions
@recipe function plot(sol::solution, x_name::Symbol=:V;linewidth=2,legend=false)
    x = getproperty(sol, x_name)
    if x isa AbstractMatrix
        x = x'
    end

    if haskey(STATE_NAMES, x_name)
        ylabel = STATE_NAMES[x_name]
    else
        ylabel = "$x_name"
    end

    if length(sol.t) ≠ length(x) error("$x_name is not in `outputs`") end

    time_unit, time_scale = time_units(sol.t[end])[2:3]

    legend --> legend
    yguide --> ylabel
    xguide --> "Time ($(time_unit))"
    linewidth --> linewidth
    sol.t./time_scale, x
end

function time_units(t::Number)
    if     t < 3600
        time_scale = 1.0
        time_unit = "s"
    elseif t < 3600*24
        time_scale = 3600.0
        time_unit = "hr"
    else
        time_scale = 3600.0*24
        time_unit = "days"
    end
    t /= time_scale
    return t, time_unit, time_scale
end

function C_rate_string(I::Number;digits::Int64=4)
    I_rat = rationalize(Float64(I))
    num = I_rat.num
    den = I_rat.den

    if den > 100 || (abs(num) > 10 && den > 10)
        return "$(round(I;digits=digits))C"
    end

    str = I < 0 ? "-" : ""

    if isone(abs(num)) && !isone(abs(den))
        str *= "C"
    else
        str *= "$(abs(num))C"
    end

    if !isone(den)
        str *= "/$den"
    end

    return str
end

function Base.show(io::IO, p::AbstractModel)
    function show_bounds(title, min, max, units="")
        if isnan(min) && isnan(max)
            return ""
        end

        pad = 19
        str = "  $title "
        if isnan(max)
            str = rpad(str * "min:",pad)*"$min$units\n"
        elseif isnan(min)
            str = rpad(str * "max:",pad)*"$max$units\n"
        else
            str = rpad(str * "bounds:",pad)*"[$min$units, $max$units]\n"
        end

        return str

    end

    # spacing
    sp = p.numerics.solid_diffusion == :Fickian ? "  " : ""

    # create the header for model
    if p isa model_skeleton
        header = "model_skeleton"
    else
        header = [x for x in remove_module_name(summary(p))]
        deleteat!(header, findall('{' .== header)[2]:length(header)-1)
        header = join(header)
    end
    header = "$(@__MODULE__) $header"

    str = string(
    "$header:\n",
    "  Cathode: $(p.numerics.cathode), $(p.numerics.rxn_p), & $(p.numerics.OCV_p)\n",
    "  Anode:   $(p.numerics.anode), $(p.numerics.rxn_n), & $(p.numerics.OCV_n)\n",
    "  System:  $(p.numerics.D_s_eff), $(p.numerics.rxn_rate), $(p.numerics.D_eff), $(p.numerics.K_eff), & $(p.numerics.thermodynamic_factor)\n",
    :model_methods ∈ fieldnames(typeof(p)) && !isempty(p.model_methods) ?
    "  Methods: $(join(p.model_methods, ", "))\n" : "",
    "  --------\n",
    "  Temperature:     $(p.numerics.temperature)\n",
    "  Solid diffusion: $(p.numerics.solid_diffusion)",
    p.numerics.solid_diffusion == :Fickian ?
    ", $(p.numerics.Fickian_method)\n" : "\n",
    "  Aging:           $(p.numerics.aging)\n",
    show_bounds("Voltage", p.bounds.V_min, p.bounds.V_max, " V"),
    show_bounds("SOC", p.bounds.SOC_min, p.bounds.SOC_max),
    show_bounds("Current", p.bounds.I_min, p.bounds.I_max, "C"),
    p.numerics.temperature ?
    show_bounds("Temperature", NaN, p.bounds.T_max-273.15, " °C") : "",
    show_bounds("Anode sat.", NaN, p.bounds.c_s_n_max),
    "  --------\n",
    p.numerics.temperature ?
    "  N.a: $sp$(p.N.a)\n" : "",
    "  N.p: $sp$(p.N.p)\n",
    p.numerics.solid_diffusion == :Fickian ?
    "  N.r_p: $(p.N.r_p)\n" : "",
    "  N.s: $sp$(p.N.s)\n",
    "  N.n: $sp$(p.N.n)\n",
    p.numerics.solid_diffusion == :Fickian ?
    "  N.r_n: $(p.N.r_n)\n" : "",
    p.numerics.temperature ?
    "  N.z: $sp$(p.N.z)\n" : "",
    )

    print(io, str[1:end-1])
end

function Base.show(io::IO, ind::indices_states)

    outputs_tot = Symbol[]
    @inbounds for (name,_type) in zip(fieldnames(typeof(ind)), fieldtypes(typeof(ind)))
        if !(_type <: Tuple)
            push!(outputs_tot, name)
        end
    end
    outputs_tot = (outputs_tot...,)

    vars = Symbol[]
    tot = Int64[]
    types = Symbol[]
    indices = UnitRange{Int64}[]

    for field in outputs_tot
        ind_var = getproperty(ind, field)
        if ind_var.var_type ∈ (:differential, :algebraic)
            push!(vars, field)
            push!(tot, ind_var.start)
            push!(types, ind_var.var_type)
            push!(indices, ind_var.start:ind_var.stop)
        end
    end
    vars    .= vars[sortperm(tot)]
    indices .= indices[sortperm(tot)]

    pad = maximum(length.(String.(vars)))+2

    str = [
        "indices_states:";
        ["  " * rpad("$(var): ", pad) * "$(length(index) > 1 ? index : index[1]), $(_type)" for (index,var,_type) in zip(indices,vars,types)]
    ]

    print(io, join(str, "\n"))
end
function Base.show(io::IO, bounds::T) where T<:AbstractStopConditions
    fields = fieldnames(T)[findall(fieldtypes(T) .<: Number)]
    vals = [getproperty(bounds,field) for field in fields]
    ind_remove = findall(.!isnan.(vals))
    fields = fields[ind_remove]
    vals = vals[ind_remove]

    pad = maximum(length.(String.(fields)))+2

    str = [
        "$T:";
        ["  " * rpad("$(field): ", pad) * "$(val)" for (field,val) in zip(fields,vals)]
    ]

    print(io, join(str, "\n"))
end
function Base.show(io::IO, opts::T) where T<:options_numerical
    fields = fieldnames(T)
    vals = [getproperty(opts,field) for field in fields]

    pad = maximum(length.(String.(fields)))+2

    str = [
        "$T:";
        ["  " * rpad("$(field): ", pad) * "$(val)" for (field,val) in zip(fields,vals)]
    ]

    print(io, join(str, "\n"))
end
function Base.show(io::IO, opts::T) where T<:AbstractOptionsModel
    fields = fieldnames(T)
    vals = [getproperty(opts,field) for field in fields]

    pad = maximum(length.(String.(fields)))+2

    str = [
        "$T:";
        ["  " * rpad("$(field): ", pad) * "$(val)" for (field,val) in zip(fields,vals)]
    ]

    print(io, join(str, "\n"))
end
method_symbol(::Type{method_I}) = :I
method_symbol(::Type{method_V}) = :V
method_symbol(::Type{method_P}) = :P

method_name(::run_constant{method_I,<:Any};        shorthand::Bool=false) = shorthand ? "I"      : "current"
method_name(::run_constant{method_V,<:Any};        shorthand::Bool=false) = shorthand ? "V"      : "voltage"
method_name(::run_constant{method_P,<:Any};        shorthand::Bool=false) = shorthand ? "P"      : "power"
method_name(::run_function{method_I,<:Function};   shorthand::Bool=false) = shorthand ? "I func"   : "current function"
method_name(::run_function{method_V,<:Function};   shorthand::Bool=false) = shorthand ? "V func"   : "voltage function"
method_name(::run_function{method_P,<:Function};   shorthand::Bool=false) = shorthand ? "P func"   : "power function"

method_string(run::run_constant{method_I,<:Any};     kw...) = method_name(run;kw...) * " = $(C_rate_string(value(run);digits=2))"
method_string(run::run_constant{method_V,<:Any};     kw...) = method_name(run;kw...) * " = $(round(value(run);digits=2)) V"
method_string(run::run_constant{method_P,<:Any};     kw...) = method_name(run;kw...) * " = $(round(value(run);digits=2)) W"

function Base.show(io::IO, result::run_results{T}) where {T<:AbstractRun}
    run = result.run
    tspan = result.tspan

    fix(x, digits=2) = round(x, digits=digits)
    str = method_string(run)
    str *= " from ($(fix(tspan[1])) s, $(fix(tspan[2])) s)"

    print(io, str)
end

function Base.show(io::IO, run::T) where {T<:AbstractRun}
    str = "Run for $(method_string(run;shorthand=true))"
    if !isempty(run.info.exit_reason)
        str *= " with exit: $(run.info.exit_reason)"
    end

    print(io, str)
end
(funcs::model_funcs)(sol::solution) = (@assert !isempty(sol); (@inbounds funcs(sol.results[end].run)))
Base.show(io::IO, ::model_funcs) = println(io,"$(@__MODULE__) model functions")
function Base.show(io::IO, sol::solution)
    results = sol.results
    if !isempty(sol)
        p = results[1].p
        Y = @views @inbounds sol.Y[end]
        t, time_unit = time_units(sol.t[end])
    end
    function str_runs()

        str = length(results) == 1 ? "  Run: " : "  Runs:"
        str *= " "^4

        methods = method_name.([result.run for result in results];shorthand=true)

        counts = ones(Int64, length(methods))

        @inbounds for i in length(counts)-1:-1:1
            if methods[i] == methods[i+1]
                counts[i] += counts[i+1]
                deleteat!(methods, i+1)
                deleteat!(counts, i+1)
            end
        end

        max_methods = 6
        show_final = 3
        methods .= @inbounds [(counts[i] == 1 ? methods[i] : "$(counts[i]) $(methods[i])") for i in 1:length(methods)]

        if length(methods) > (max_methods+1)
            first_range = 1:(max_methods-show_final)
            final_range = (length(counts) - (show_final-1)):length(counts)
            count_not_showed = sum(counts) - sum(counts[[first_range;final_range]])

            methods = [
                methods[first_range];
                "…$(count_not_showed)…";
                methods[final_range];
            ]
        end
        str *= join(methods, " → ") * "\n"

        return str
    end

    title = "$(@__MODULE__) simulation"
    if !isempty(sol)
        str = @views @inbounds string(
                "$title\n",
                "  --------\n",
                str_runs(),
                "  Time:    $(round(t;                  digits = 2)) $time_unit\n",
                "  Current: $(C_rate_string(calc_I(Y,p);digits = 4))\n",
                "  Voltage: $(round(calc_V(Y,p);        digits = 4)) V\n",
                "  Power:   $(round(calc_P(Y,p);        digits = 4)) W\n",
                "  SOC:     $(round(sol.SOC[end];     digits = 4))\n",
                !(p.numerics.aging == false) ?
                "  SOH:     $(round(sol.SOH[end];     digits = 4))\n"
                : "",
                !(p.numerics.temperature == false) ?
                "  Temp.:   $(round(temperature_weighting(calc_T(Y,p),p)-273.15; digits = 4)) °C\n"
                : "",
                "  Exit:    $(sol.results[end].info.exit_reason)",
            )
    else
        str = "$title: empty"
    end

    print(io, str)
end

Base.show(io::IO, ind::states_logic) = print(io, remove_module_name("states_logic using $(ind.results)"))

Base.deleteat!(a::VectorOfArray, i::Integer) = (Base._deleteat!(a.u, i, 1); a.u)

function emptyfunc end


text:
function recursive_replace(dir, (pattern, substitute))
    for (root, dirs, files) in walkdir(dir)
        for f in files
            if !isnothing(match(pattern, f))
                cp(
                    joinpath(root, f),
                    joinpath(root, replace(f, pattern => substitute))
                    ;force=true
                )
            end
        end
    end
end

#=
```@blog_meta
last_update="2021-08-13"
create="2021-08-12"
download_link="../x.txt"
tags=["Archive", "记录", "Test"]
```
=#

using Badges

function blog_meta(;
    create=nothing,
    last_update = nothing,
    download_link = nothing,
    tags = []
)
    html = ""

    if !isnothing(last_update)
        html *= Badge(
            label="Last Update",
            message=last_update,
            color="#97C40F"
        ) |> Badges.render
    end

    if !isnothing(create)
        html *= Badge(
            label="Create",
            message=create
        ) |> Badges.render
    end

    if !isnothing(download_link)
        html *= Badge(
            label="⇩",
            message=basename(download_link),
            rightLink=download_link
        ) |> Badges.render
    end

    for t in tags
        html *= Badge(
            label= "#",
            message=t,
            color="#0F80C1",
            rightLink="/search/?q=$t"
        ) |> Badges.render
    end

    "
$html
"
end

using Documenter: Utilities, Expanders, Documents
using Documenter.Utilities: Selectors
using Documenter.Expanders: ExpanderPipeline, iscode

abstract type BlogMetaBlocks <: ExpanderPipeline end
Selectors.order(::Type{BlogMetaBlocks}) = 20.0
Selectors.matcher(::Type{BlogMetaBlocks}, node, page, doc) = iscode(node, r"^@blog_meta")

function Selectors.runner(::Type{BlogMetaBlocks}, x, page, doc)
    curmod = get(page.globals.meta, :CurrentModule, Main)
    fields = Dict{Symbol, Any}()
    lines = Utilities.find_block_in_file(x.code, page.source)
    @debug "Evaluating @blog_meta block:\n$(x.code)"

    for (ex, str) in Utilities.parseblock(x.code, doc, page)
        if Utilities.isassign(ex)
            try
                fields[ex.args[1]] = Core.eval(curmod, ex.args[2])
            catch err
                push!(doc.internal.errors, :autodocs_block)
                @warn("""
                    failed to evaluate `$(strip(str))` in `@blog_meta` block in $(Utilities.locrepr(page.source, lines))
                    ```$(x.language)
                    $(x.code)
                    ```
                    """, exception = err)
            end
        end
    end

    if haskey(fields, :create)
        page.mapping[x] = blog_meta(
            ;create=get(fields, :create, nothing),
            last_update=get(fields, :last_update, nothing),
            download_link=get(fields, :download_link, nothing),
            tags=get(fields, :tags, [])
        ) |> Documents.RawHTML
    end
end

#=
```@inline ../path/to/your.html
```
=#

abstract type InlineBlocks <: ExpanderPipeline end
Selectors.order(::Type{InlineBlocks}) = 21.0
Selectors.matcher(::Type{InlineBlocks}, node, page, doc) = iscode(node, r"^@inline")

function Selectors.runner(::Type{InlineBlocks}, x, page, doc)
    m = split(x.language)
    if length(m) == 2
        f = m[2]
        page.mapping[x] = Documents.RawHTML("")
    end
end

#=
```@comment
```
=#

const COMMENT_TEMPLATE = """

    """

abstract type CommentBlocks <: ExpanderPipeline end
Selectors.order(::Type{CommentBlocks}) = 22.0
Selectors.matcher(::Type{CommentBlocks}, node, page, doc) = iscode(node, r"^@comment")

function Selectors.runner(::Type{CommentBlocks}, x, page, doc)
    page.mapping[x] = Documents.RawHTML(COMMENT_TEMPLATE)
end

#=
```@embed https://github.com/findmyway/TianJun.jl/blob/ca409281346785f779fc88b770ee8a10f5d07ea2/LICENSE#L1-L21
```
=#

using Markdown
using URIs
using Downloads

abstract type EmbedBlocks <: ExpanderPipeline end
Selectors.order(::Type{EmbedBlocks}) = 23.0
Selectors.matcher(::Type{EmbedBlocks}, node, page, doc) = iscode(node, r"^@embed")

function Selectors.runner(::Type{EmbedBlocks}, x, page, doc)
    m = split(x.language)
    if length(m) == 2
        L = URI(m[2])
        (user, repo, _, branch, filepath...) = URIs.splitpath(L)
        lines = Downloads.download(
            "$(L.scheme)://raw.githubusercontent.com/$user/$repo/$branch/$(join(filepath, '/'))";
        ) |> readlines

        m = match(r"L(?<start>\d+)(-L(?<end>\d+))?", L.fragment)
        if isnothing(m)
            s, e = 1, length(lines)
        else
            s, e = m[:start], m[:end]
            s = isnothing(s) ? 1 : parse(Int, s)
            e = isnothing(e) ? length(lines) : parse(Int, e)
        end

        lang = lstrip(splitext(filepath[end])[end], '.')
        if lang == "jl"
            lang = "julia"
        end

        page.mapping[x] = Documents.MultiOutput(
            [
                Markdown.Code(lang, join(lines[s:e], "\n")),
                Documents.RawHTML("""

                    ❤
                    Source: $(filepath[end])
                    ❤

                    """)
            ]
        )
    end
end

max_stars_repo_name: pilgrimygy/Pilgrimygy

max_issues_repo_name: pilgrimygy/Pilgrimygy

max_forks_repo_name: pilgrimygy/Pilgrimygy


text:
using Healpix
using PowerSpectra
using Plots

# retrieve the downgraded planck maps
nside = 256
map1 = 1e6 * PowerSpectra.planck256_map("100", "hm1", "I_STOKES")
map2 = 1e6 * PowerSpectra.planck256_map("100", "hm2", "I_STOKES")
plot(map1, clim=(-200,200))

mask1 = PowerSpectra.planck256_mask("100", "hm1", :T)
mask2 = PowerSpectra.planck256_mask("100", "hm2", :T)
plot(mask1)

# compute the cross-spectra of the masked maps
masked_alm_1 = map2alm(map1 * mask1)
masked_alm_2 = map2alm(map2 * mask2)
pCl = SpectralVector(alm2cl(masked_alm_1, masked_alm_2))

# compute the mode coupling matrix
M_TT = mcm(:TT, map2alm(mask1), map2alm(mask2))

# perform a linear solve to undo the effects of mode-coupling
Cl = M_TT \ pCl

# get lmax from nyquist frequency
lmax = nside2lmax(nside)

# get the planck instrumental beam and pixel window
Wl = PowerSpectra.planck_beam_Wl("100", "hm1", "100", "hm2", :TT, :TT; lmax=lmax)
pixwinT = SpectralVector(pixwin(nside)[1:(lmax+1)])

# plot our spectra
ell = eachindex(Cl)
prefactor = ell .* (ell .+ 1) ./ (2π)
plot(prefactor .*  Cl ./ (Wl .* pixwinT.^2), label="\$D_{\\ell}\$", xlim=(0,2nside) )

# compare it to theory
theory = PowerSpectra.planck_theory_Dl()  # returns a Dict of Dl indexed with :TT, :EE, ...
plot!(theory[:TT], label="theory TT")

##
using Healpix
using PowerSpectra
using Plots

nside = 256
m₁ = PowerSpectra.planck256_polmap("100", "hm1")
m₂ = PowerSpectra.planck256_polmap("100", "hm2")
maskT₁ = PowerSpectra.planck256_mask("100", "hm1", :T)
maskP₁ = PowerSpectra.planck256_mask("100", "hm1", :P)
maskT₂ = PowerSpectra.planck256_mask("100", "hm2", :T)
maskP₂ = PowerSpectra.planck256_mask("100", "hm2", :P)

# convert to μK
scale!(m₁, 1e6)
scale!(m₂, 1e6)

lmax = nside2lmax(256)
# utility function for doing TEB mode decoupling
Cl = master(m₁, maskT₁, maskP₁,
            m₂, maskT₂, maskP₂)
lmax = nside2lmax(256)
print(keys(Cl))
spec = :TE
Wl = PowerSpectra.planck_beam_Wl("100", "hm1", "100", "hm2", spec, spec; lmax=lmax)
pixwinT, pixwinP = pixwin(nside; pol=true)
pixwinT = SpectralVector(pixwinT[1:(lmax+1)])
pixwinP = SpectralVector(pixwinP[1:(lmax+1)])

ℓ = eachindex(Wl)
plot( (ℓ.^2 / (2π)) .*  Cl[spec] ./ (Wl .* pixwinT .* pixwinP),
    label="\$D_{\\ell}\$", xlim=(0,2nside) )
theory = PowerSpectra.planck_theory_Dl()
plot!(theory[spec], label="theory $(spec)")
##

spec = :EE
Wl = PowerSpectra.planck_beam_Wl("100", "hm1", "100", "hm2", spec, spec; lmax=lmax)
ℓ = eachindex(Wl)
plot( (ℓ.^2 / (2π)) .*  Cl[spec] ./ (Wl .* pixwinP.^2),
    label="\$D_{\\ell}\$", xlim=(0,2nside) )
theory = PowerSpectra.planck_theory_Dl()
plot!(theory[spec], label="theory $(spec)", ylim=(-10,50))
##

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
# VWN parameterization of the exchange correlation energy
function excVWN( n::Array{Float64,1} )
    # Constants
    X1 = 0.75*(3.0/(2.0*pi))^(2.0/3.0)
    A = 0.0310907
    x0 = -0.10498
    b = 3.72744
    c = 12.9352
    Q = sqrt(4*c-b*b)
    X0 = x0*x0 + b*x0 + c
    Npoints = size(n)[1]
    out = zeros(Float64,Npoints)
    for ip = 1:Npoints
        rs = (4.0*pi/3.0*n[ip])^(-1.0/3.0) # Added internal conversion to rs
        x = sqrt(rs)
        X = x*x + b*x + c
        out[ip] = -X1/rs + A*( log(x^2/X) + 2*b/Q*atan(Q/(2*x + b))
                  -(b*x0)/X0*( log((x-x0)*(x-x0)/X) + 2*(2*x0+b)/Q*atan(Q/(2*x + b)) ))
    end
    return out
end

# d/dn deriv of VWN parameterization of the exchange correlation energy
function excpVWN( n::Array{Float64,1} )
    # Constants
    X1 = 0.75*(3.0/(2.0*pi))^(2.0/3.0)
    A = 0.0310907
    x0 = -0.10498
    b = 3.72744
    c = 12.9352
    Q = sqrt(4.0*c - b*b)
    X0 = x0*x0 + b*x0 + c
    #
    Npoints = size(n)[1]
    out = zeros(Float64,Npoints)
    for ip = 1:Npoints
        rs = (4.0*pi/3.0*n[ip])^(-1.0/3.0) # Added internal conversion to rs
        x = sqrt(rs)
        X = x*x + b*x + c
        dx = 0.5/x # Chain rule needs dx/drho!
        out[ip] = dx*( 2.0*X1 / (rs*x)
                  + A*( 2.0/x - (2*x+b)/X - 4*b/(Q*Q + (2*x+b)*(2*x+b))
                  - (b*x0)/X0*( 2.0/(x-x0) - (2*x+b)/X - 4*(2*x0+b) /
                    (Q*Q + (2*x+b) * (2*x+b)) ) ))
        out[ip] = (-rs/(3*n[ip]))*out[ip] # Added d(rs)/dn from chain rule from rs to n conv
    end
    return out
end


text:
"""
    read_idxhead(io::IO)

Read and verify the header of a git packfile index. Only supports v2 indices.
"""
function read_idxhead(io::IO)
    magic = ntoh(read(io, UInt32))
    @assert magic == 0xff_74_4f_63
    version = ntoh(read(io, UInt32))
    @assert version == 2
end

"""
    lookup_idx(filename, hash::SHA1)

Lookup the id `hash` is in the packfile index file `filename`. If contained in the index,
return the absolute offset in the packfile, otherwise return -1.
"""
function lookup_idx(filename, hash::SHA1)
    byte1 = hash.bytes[1]
    open(filename) do io
        # 1. Header
        read_idxhead(io)

        # 2. Fan
        if byte1 == 0
            i = zero(UInt32)
        else
            skip(io, 4*(byte1-1))
            i = ntoh(read(io, UInt32))
        end
        skip(io, 4*(255-byte1))
        i_max = ntoh(read(io, UInt32)) # final element of fan

        # 3. SHA1
        skip(io, 20*i) # start point
        while true
            if i >= i_max
                # hash is not in index
                return Int64(-1)
            end
            i_hash = read(io, SHA1)
            if i_hash == hash
                break
            elseif i_hash > hash
                # hash is not in index
                return Int64(-1)
            end
            i += 1
        end
        skip(io, 20*(i_max - i - 1))

        # 4. CRC
        skip(io, 4*i_max)

        # 5. offset
        skip(io, 4*i)
        j = ntoh(read(io, UInt32))
        if j < 0x8000_0000
            # offsets < 2GB
            return Int64(j)
        end
        skip(io, 4*(i_max - i - 1))

        # 6. large files only
        j -= 0x8000_0000
        skip(io, 8*j)
        return Int64(ntoh(read(io, UInt64)))
    end
end

function read_idx(filename)
    byte1 = hash.bytes[1]
    open(filename) do f
        # 1. Header
        read_idxhead(io)

        # 2. Fan
        skip(io, 4*255)
        i_max = ntoh(read(io, UInt32))

        # 3. SHA1
        hashes = read(io, SHA1, i_max)

        # 4. offset
        offsets = Int64.(ntoh.(read(io, UInt32, i_max)))
        m = maximum(offsets)

        if m >= 0x8000_0000
            bigoffsets = ntoh.(read(io, Int64, m - 0x8000_0000 + 1))
            for (i,o) in enumerate(offsets)
                if o >= 0x8000_0000
                    offsets[i] = bigoffsets[o - 0x8000_0000 + 1]
                end
            end
        end

        return Dict(zip(hashes, offsets))
    end
end

"""
    read_packhead(io::IO)

Read and verify the header of a git packfile, returning the number of objects
contained. Only supports v2 format.
"""
function read_packhead(io::IO)
    magic = ntoh(read(io, UInt32))
    @assert magic == 0x50_41_43_4b # "PACK"
    version = ntoh(read(io, UInt32))
    @assert version == 2
    nobjs = ntoh(read(io, UInt32))
    return nobjs
end

"""
    t,n = read_packobjhead(io::IO)

Read and verify the individual object header from a git packfile, returning the object type code `t` and size `n`. `t` is one of
 - `OBJ_COMMIT`
 - `OBJ_TREE`
 - `OBJ_BLOB`
 - `OBJ_TAG`
 - `OBJ_OFS_DELTA`
 - `OBJ_REF_DELTA`
"""
function read_packobjhead(io::IO)
    b = read(io,UInt8)
    ismore = b >= 0b1000_0000       # more to follow
    t = ObjCode((b >> 4) & 0b0111)  # object type
    n = Int64(b & 0b1111)           # lowest order bits
    shft = 4
    while ismore
        b = read(io,UInt8)
        ismore = b >= 0b1000_0000
        n = (Int64(b & 0b0111_1111) << shft) | n
        shft += 7
    end
    return t,n
end

"""
    read_offset_varint(io::IO)

Read the variable-length integer format used to encode the offset of `OBJ_OFS_DELTA`
objects.
"""
function read_offset_varint(io::IO)
    b = read(io,UInt8)
    ismore = b >= 0b1000_0000
    n = Int64(b & 0b0111_1111)
    while ismore
        b = read(io,UInt8)
        ismore = b >= 0b1000_0000
        n = ((n+1) << 7) | (b & 0b0111_1111)
    end
    return n
end

"""
    read_delta_varint(io::IO)

Read the variable-length integer format used _inside_ the delta objects to store length of
source and target.
"""
function read_delta_varint(io::IO)
    b = read(io,UInt8)
    ismore = b >= 0b1000_0000
    shft = 7
    n = Int64(b & 0b0111_1111)
    while ismore
        b = read(io,UInt8)
        ismore = b >= 0b1000_0000
        n |= (Int64(b & 0b0111_1111) << shft)
        shft += 7
    end
    return n
end

"""
    getobjdata_pack(repo::GitRepo, packfile::String, offset::Integer)

Return the type code `t` and contents `data` of the object located at `offset` in
`packfile`. `repo` is required to lookup the source for objects stored as `OBJ_REF_DELTA`.
"""
function getobjdata_pack(repo::GitRepo, packfile::String, offset::Integer)
    open(packfile) do io
        read_packhead(io)
        seek(io, offset)
        getobjdata_pack(repo, packfile, offset, io)
    end
end

function getobjdata_pack(repo::GitRepo, packfile::String, offset::Integer, io::IO)
    t, len = read_packobjhead(io)
    if 0x01 <= UInt8(t) <= 0x04
        zio = ZlibDecompressorStream(io, stop_on_end=true)
        return GitRawObject(t, read(zio, len))
    else
        if t == Obj.offset_delta
            reloffset = read_offset_varint(io)
            rawsrc = getobjdata_pack(repo, packfile, offset - reloffset)
        elseif t == Obj.ref_delta
            srchash = read(io, SHA1)
            rawsrc = getobjdata(repo, refhash)
        else
            error("Invalid tag $t")
        end

        src = IOBuffer(rawsrc.data)
        delta = ZlibDecompressorStream(io, stop_on_end=true)
        sourcelen = read_delta_varint(delta)
        targetlen = read_delta_varint(delta)

        out = IOBuffer(targetlen)
        while position(out) < targetlen
            op = read(delta, UInt8)
            if op >= 0b1000_0000
                # copy operation
                u = UInt32(0) # position
                for i = 0:3
                    if op & (0b0000_0001 << i) != 0
                        u |= UInt32(read(delta, UInt8)) << (8*i)
                    end
                end

                k = UInt32(0) # number of bytes to copy
                for i = 0:2
                    if op & (0b0001_0000 << i) != 0
                        k |= UInt32(read(delta, UInt8)) << (8*i)
                    end
                end

                seek(src, Int(u))
                write(out, read(src, Int(k)))
            else
                # insert operation
                write(out, read(delta, Int(op)))
            end
        end
        return GitRawObject(tt, take!(out))
    end
end

function idxname(packfile::String)
    endswith(packfile, ".pack") || error("invalid packfile name $packfile")
    packfile[1:end-4]*"idx"
end

function buildindex(repo, packfile, idxfile=idxname(packfile))
    open(packfile) do pack
        open(idxfile,"w") do idx
            d = SortedDict{SHA1, Tuple{Int64, UInt32}}()

            nobjs = read_packhead(pack)
            for i = 1:nobjs
                raw = getobjdata_pack(repo, packfile, position(pack), pack)
                hash = oid(raw)
                crc  = crc32c(raw.data)
                d[hash] => (length(data), crc)
            end
        end
    end
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX