text:
# Figure 5.3a in "Infinite GMRES for parameterized linear systems"

using SparseArrays
using LinearAlgebra
using PyPlot
using MAT
pygui(true)
include("iar.jl");
include("infgmres.jl");

# Problem represented with explicitly available derivatives
struct GeneralProblem
    sigma
    Aders
    Afun
    b
end

function compute_Mlincomb(prob::GeneralProblem,sigma,Y)
    if (prob.sigma != sigma)
        @show sigma
        @show prob.sigma
        error("Wrong sigma");
    end
    z=zeros(eltype(Y),size(Y,1));
    for i=1:size(Y,2)
        z += prob.Aders[i]*Y[:,i];
    end
    return z;
end
function compute_Mder(prob::GeneralProblem,mu)
    return prob.Afun(mu);
end

## Define problem
data=matread("delay_problem_matrices.mat")
A0=data["A0"];
A1=data["A1"];
b=vec(data["b"]);
n=length(b);

# Identity matrix
II=Matrix{Float64}(I,n,n);

Afun= s-> -s*II+A0+A1*exp(-s);
sigma=0;
N=100;
ders=Vector{Matrix{Float64}}(undef,N);

for i=1:N
    ders[i]=((-1)^(i+1))*A1;
end
ders[1] += A0;
ders[2] += -II;

myproblem=GeneralProblem(sigma,ders,Afun,b);

## Run Infinite GMRES
(f,hist1)=infgmres(myproblem,30,method=:iar);
mu_list = range(1.e-5,stop=1.e-1,length=20);
mu_list = mu_list[2:end];
how_many = zeros(length(mu_list),1);
tol = 1.e-12;

for i=1:length(mu_list)
    global mu
    mu = mu_list[i];
    A_of_mu = myproblem.Afun(mu);
    res = zeros(100,1);
    for m=1:200
        global f
        (xx)=f(mu,m)
        res[m] = norm(A_of_mu*xx - myproblem.b)/norm(myproblem.b);
        if res[m] < tol
            how_many[i] = m;
            break;
        end
    end
end

figure(1)
ax = PyPlot.gca()
xlabel("mu")
ylabel("its")
plot(mu_list,how_many,"ro")

ext: jl

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
function sumlengths(points)
    T = eltype(points[1])
    result = zeros(T, length(points))
    for i=1:length(points)
        i0 = max(i-1,1)
        p1, p2 = points[i0], points[i]
        if !(any(map(isnan, p1)) || any(map(isnan, p2)))
            result[i] = result[i0] + norm(p1-p2)
        else
            result[i] = result[i0]
        end
    end
    result
end

intensity_convert(intensity, verts) = intensity
function intensity_convert(intensity::VecOrSignal{T}, verts) where T
    if length(to_value(intensity)) == length(to_value(verts))
        GLBuffer(intensity)
    else
        Texture(intensity)
    end
end
function intensity_convert_tex(intensity::VecOrSignal{T}, verts) where T
    if length(to_value(intensity)) == length(to_value(verts))
        TextureBuffer(intensity)
    else
        Texture(intensity)
    end
end
#TODO NaNMath.min/max?
dist(a, b) = abs(a-b)
mindist(x, a, b) = min(dist(a, x), dist(b, x))
function gappy(x, ps)
    n = length(ps)
    x <= first(ps) && return first(ps) - x
    for j=1:(n-1)
        p0 = ps[j]
        p1 = ps[min(j+1, n)]
        if p0 <= x && p1 >= x
            return mindist(x, p0, p1) * (isodd(j) ? 1 : -1)
        end
    end
    return last(ps) - x
end
function ticks(points, resolution)
    Float16[gappy(x, points) for x = range(first(points), stop=last(points), length=resolution)]
end

# ambigious signature
function _default(position::VectorTypes{<: Point}, s::style"lines", data::Dict)
    line_visualization(position, data)
end
function _default(position::MatTypes{<:Point}, s::style"lines", data::Dict)
    line_visualization(position, data)
end
function line_visualization(position::Union{VectorTypes{T}, MatTypes{T}}, data::Dict) where T<:Point
    p_vec = if isa(position, GPUArray)
        position
    else
        const_lift(vec, position)
    end

    @gen_defaults! data begin
        dims::Vec{2, Int32} = const_lift(position) do p
            sz = ndims(p) == 1 ? (length(p), 1) : size(p)
            Vec{2, Int32}(sz)
        end
        vertex              = p_vec => GLBuffer
        intensity           = nothing
        color_map           = nothing => Texture
        color_norm          = nothing
        color               = (color_map == nothing ? default(RGBA, s) : nothing) => GLBuffer
        thickness::Float32  = 2f0
        pattern             = nothing
        fxaa                = false
        preferred_camera    = :orthographic_pixel
        # Duplicate the vertex indices on the ends of the line, as our geometry
        # shader in `layout(lines_adjacency)` mode requires each rendered
        # segment to have neighbouring vertices.
        indices             = const_lift((p)->[1; 1:length(p); length(p)], p_vec) => to_index_buffer
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "lines.vert", "lines.geom", "lines.frag")
        gl_primitive        = GL_LINE_STRIP_ADJACENCY
        valid_vertex        = const_lift(p_vec) do pv
            map(p-> Float32(all(isfinite, p)), pv)
        end => GLBuffer
    end
    if pattern != nothing
        if !isa(pattern, Texture)
            if !isa(pattern, Vector)
                error("Pattern needs to be a Vector of floats. Found: $(typeof(pattern))")
            end
            tex = GLAbstraction.Texture(ticks(pattern, 100), x_repeat = :repeat)
            data[:pattern] = tex
        end
        @gen_defaults! data begin
            pattern_length = Float32(last(pattern))
            lastlen   = const_lift(sumlengths, p_vec) => GLBuffer
            maxlength = const_lift(last, lastlen)
        end
    end
    data[:intensity] = intensity_convert(intensity, vertex)
    data
end

to_points(x::Vector{LineSegment{T}}) where {T} = reinterpret(T, x, (length(x)*2,))

_default(positions::VectorTypes{LineSegment{T}}, s::Style, data::Dict) where {T <: Point} =
    _default(const_lift(to_points, positions), style"linesegment"(), data)

function _default(positions::VectorTypes{T}, s::style"linesegment", data::Dict) where T <: Point
    @gen_defaults! data begin
        vertex              = positions           => GLBuffer
        color               = default(RGBA, s, 1) => GLBuffer
        color_map           = nothing => Texture
        color_norm          = nothing
        thickness           = 2f0                 => GLBuffer
        shape               = RECTANGLE
        pattern             = nothing
        fxaa                = false
        indices             = const_lift(length, positions) => to_index_buffer
        # TODO update boundingbox
        shader              = GLVisualizeShader("fragment_output.frag", "util.vert", "line_segment.vert", "line_segment.geom", "lines.frag")
        gl_primitive        = GL_LINES
    end
    if !isa(pattern, Texture) && pattern != nothing
        if !isa(pattern, Vector)
            error("Pattern needs to be a Vector of floats")
        end
        tex = GLAbstraction.Texture(ticks(pattern, 100), x_repeat = :repeat)
        data[:pattern] = tex
        data[:pattern_length] = Float32(last(pattern))
    end
    data
end

function _default(positions::Vector{T}, range::AbstractRange, s::style"lines", data::Dict) where T <: AbstractFloat
    length(positions) != length(range) && throw(
        DimensionMismatsch("length of $(typeof(positions)) $(length(positions)) and $(typeof(range)) $(length(range)) must match")
    )
    _default(points2f0(positions, range), s, data)
end

function line_indices(array)
    len = length(array)
    result = Array(GLuint, len*2)
    idx = 1
    for i=0:(len-3), j=0:1
        result[idx] = i+j
        idx += 1
    end
    result
    #GLuint[i+j for i=0:(len-3) for j=0:1] # on 0.5
end
"""
Fast, non anti aliased lines
"""
function _default(position::VectorTypes{T}, s::style"speedlines", data::Dict) where T <: Point
    @gen_defaults! data begin
        vertex       = position => GLBuffer
        color_map    = nothing  => Vec2f0
        indices      = const_lift(line_indices, position) => to_index_buffer
        color        = (color_map == nothing ? default(RGBA{Float32}, s) : nothing) => GLBuffer
        color_norm   = nothing  => Vec2f0
        intensity    = nothing  => GLBuffer
        shader       = GLVisualizeShader("fragment_output.frag", "dots.vert", "dots.frag")
        gl_primitive = GL_LINES
    end
end


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2012
# Define operations for Array spaces
#--------------------------------------------------------------------

"""
    Extract a particular field from the array of Tuples  (AoT)
    Also handles exceptions for when the array index does not contain an array
    but a product space.
    Input: Array of Tuples (AoT)
    Output: Array of Field (AoF)
"""
function extract(AoT::Array{T, 2}, ID::Int)::Array{Field, 2} where T<:Union{ProductSpace, NTuple{N, Field}} where {N}
    AoF = Array{Field, 2}(undef, size(AoT))
    for index in CartesianIndices(AoT)
        try
            AoF[index] = AoT[index][ID]
        catch
            AoF[index] = Field(AoT[index], Array{Missing, 2}(missing, size(AoT[index])))
        end
    end
    return AoF
end

"""
 Find the minimum of a field over all the patches
 Input: Array of Field (AoF)
 Output: minimum (scalar)
"""
function Base. minimum(AoF::Array{Field, 2})::Number
    VoF = collect(Iterators.flatten(value.(AoF)))
    return minimum(skipmissing(VoF))
end

"""
 Find the maximum of a field over all the patches
 Input: Array of Field (AoF)
 Output: maximum (scalar)
"""
function Base. maximum(AoF::Array{Field, 2})::Number
    VoF = collect(Iterators.flatten(value.(AoF)))
    return maximum(skipmissing(VoF))
end

"""
    Takes in an array of spaces, and returns the value of the function on all of those spaces
    Input: Array of spaces (AoS)
    Output: Array of fields (AoF)
"""
function Field(AoS::Array{Union{ProductSpace, NTuple{N,Field}}, 2}, map::Function)  where {N}
    AoF = Array{Field, 2}(undef, size(AoS))
    for index in CartesianIndices(AoF)
        AoF[index] = Field(AoS[index], map)
    end
    return AoF
end

"""
Compute the root-mean-square error over the whole grid.
For each individual patch use spectral integration, and for the whole
grid use a standard L2 norm.
"""
function rmse(AoF::Matrix{Field})::Real
    return norm(norm.(AoF) / size(AoF)[1])
end

function Base. *(λ::Number, U::NTuple{N, Field})::NTuple{N, Field} where {N}
    return map(u->λ*u, U)
end

function Base. +(U::NTuple{N, Field{S}}, V::NTuple{N, Field{S}})::NTuple{N, Field{S}} where {N, S}
    return U .+ V
end

function Base. -(U::NTuple{N, Field{S}}, V::NTuple{N, Field{S}})::NTuple{N, Field{S}} where {N, S}
    return U .- V
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 14a4cc28-ab71-11eb-0021-1506fbc8ae2f
using Plots, LaTeXStrings

# ╔═╡ 639c60aa-be9e-4973-849d-345a61234663
md"# Chapter 1 Figures for rates"

# ╔═╡ d47fe2e9-07a1-492a-a6cd-e4a9f92cb5aa
begin
	mma(x) = 1/(1+x)
	mmr(x) = 1/(1+1/x)
	hfr(n) = x -> 1/(1 + x^n)
	hfa(n) = x -> 1/(1 + (1/x)^n)
	hs(λ, n) = x-> hfr(n)(x) + λ * hfa(n)(x)
end

# ╔═╡ 260fcd7d-ebf4-4859-8121-6455ff1e3c81
x = 0:0.01:2.0

# ╔═╡ 0dd7b136-8ee0-4d60-b8e5-c12094c5ac71
m2 = plot(mma, x, legend=:none, ylimits = (0. ,1.), xformatter=_->"", yformatter=_->"", title="Michaelis-Menten Repressor", yticks=[0, 0.5, 1.0])

# ╔═╡ c08d7b09-a6e4-4107-a9a8-54484f4f5829
m1 = plot(mmr, x, legend=:none, ylimits = (0. ,1.),
	xformatter=_->"",
	title="Michaelis-Menten Activator", ylabel="Normalized Promoter Activity", yticks=[0, 0.5, 1.0])

# ╔═╡ 47b79606-b710-4928-88e3-902dda002560
begin
	plot()
	for n = 2:15
		plot!(hfa(n), x, label=L"n = %$n")
	end
	m3 = plot!(legend=:none, xlabel=L"\frac{S_X}{K_X}", ylimits = (0. ,1.), title="Hill Activator", ylabel="Normalized Promoter Activity", yticks=[0, 0.5, 1.0])
end

# ╔═╡ 28d6f3cd-af64-43e0-b951-a104955e89b2
begin
	plot()
	for n = 2:15
		plot!(hfr(n), x)
	end
	m4 = plot!(legend=:none, xlabel=L"\frac{S_X}{K_X}", ylimits = (0. ,1.),
		yformatter=_->"",
		title="Hill Repressor", yticks=[0, 0.5, 1.0])
end

# ╔═╡ 3b79d3b7-7a58-4b88-82cf-507776972aca
begin
	plot()
	plot!(hs(.5, 3), x, label=L"\lambda = 0.5")
	plot!(hs(1, 3), x, label=L"\lambda = 1")
	plot!(hs(2, 3), x, label=L"\lambda = 2")
	shil = plot!(xlabel=L"\frac{S_X}{K_X}",ylabel="Normalized Promoter Activity",title="Shifted Hill", yticks=[0, 0.5, 1.0, 1.5, 2.0])
	savefig(shil, "ch1_shifted.pdf")
	shil
end

# ╔═╡ 573d16f2-b276-4a34-9ec1-c12e3b59af9f
begin
	figactivity = plot(m1,m2,m3,m4, layout=4, size=(750,750))
	figactivity
end

# ╔═╡ a7ef3696-ff92-4264-9dcc-4a38f21d9852
# savefig(figactivity, "ch1_activity.pdf")

# ╔═╡ Cell order:
# ╠═14a4cc28-ab71-11eb-0021-1506fbc8ae2f
# ╟─639c60aa-be9e-4973-849d-345a61234663
# ╠═d47fe2e9-07a1-492a-a6cd-e4a9f92cb5aa
# ╠═260fcd7d-ebf4-4859-8121-6455ff1e3c81
# ╠═0dd7b136-8ee0-4d60-b8e5-c12094c5ac71
# ╠═c08d7b09-a6e4-4107-a9a8-54484f4f5829
# ╠═47b79606-b710-4928-88e3-902dda002560
# ╠═28d6f3cd-af64-43e0-b951-a104955e89b2
# ╠═3b79d3b7-7a58-4b88-82cf-507776972aca
# ╠═573d16f2-b276-4a34-9ec1-c12e3b59af9f
# ╠═a7ef3696-ff92-4264-9dcc-4a38f21d9852


text:
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile

#########
# setup #
#########

# some objective function to work with
f(a, b) = sum(a' * b + a * b')

# pre-record a GradientTape for `f` using inputs of shape 100x100 with Float64 elements
const f_tape = GradientTape(f, (rand(100, 100), rand(100, 100)))

# compile `f_tape` into a more optimized representation
const compiled_f_tape = compile(f_tape)

# some inputs and work buffers to play around with
a, b = rand(100, 100), rand(100, 100)
inputs = (a, b)
results = (similar(a), similar(b))
all_results = map(DiffResults.GradientResult, results)
cfg = GradientConfig(inputs)

####################
# taking gradients #
####################

# with pre-recorded/compiled tapes (generated in the setup above) #
#-----------------------------------------------------------------#

# this should be the fastest method, and non-allocating
gradient!(results, compiled_f_tape, inputs)

# the same as the above, but in addition to calculating the gradients, the value `f(a, b)`
# is loaded into the the provided `DiffResult` instances (see DiffResults.jl documentation).
gradient!(all_results, compiled_f_tape, inputs)

# this should be the second fastest method, and also non-allocating
gradient!(results, f_tape, inputs)

# you can also make your own function if you want to abstract away the tape
∇f!(results, inputs) = gradient!(results, compiled_f_tape, inputs)

# with a pre-allocated GradientConfig #
#-------------------------------------#
# these methods are more flexible than a pre-recorded tape, but can be
# wasteful since the tape will be re-recorded for every call.

gradient!(results, f, inputs, cfg)

gradient(f, inputs, cfg)

# without a pre-allocated GradientConfig #
#----------------------------------------#
# convenient, but pretty wasteful since it has to allocate the GradientConfig itself

gradient!(results, f, inputs)

gradient(f, inputs)


text:
using Meshing
using GeometryTypes
using Juno
using Profile

# algo type to profile
algo = MarchingTetrahedra()

# set if using function or SDF variant
fn_mesh = false

## Constructions
torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)),0.05) do v
    (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus
end

torus_fn(v) = (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus

# clearout profile
Profile.clear()

# warmup and profile
if fn_mesh
    HomogenousMesh(torus_fn,HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), algo) # use default 128 size for warmup
    Juno.@profiler HomogenousMesh(torus_fn,HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (512,512,512), algo)
else
    HomogenousMesh(torus, algo)
    Juno.@profiler HomogenousMesh(torus, algo)
end

Juno.profiler()


text:
using KNITRO, JuMP, Compat.Test

m = Model(solver=KnitroSolver())
@variable(m, x)
setvalue(x, 5)
@NLobjective(m, Min, -log(x)+x^2)
solve(m)
@test isapprox(getobjectivevalue(m), 0.84657, atol=1.0e-5)
@test isapprox(getvalue(x), 0.707107, atol=1.0e-5)

ktrmod = internalmodel(m)
MathProgBase.freemodel!(ktrmod)


text:
using Opcode
using Base.Test

using PlyIO

function imwritepnm(x::Array{T,2},fname) where {T <: Real}
    width, height = size(x)
    xmax = maximum(x)
    open(fname,"w") do stream
        write(stream, "P6\n$width $height\n255\n")
        for j in 1:height, i in 1:width, k = 1:3
            write(stream,round(UInt8,255*x[i,j]/xmax))
        end
    end
end

function readply(filename)
    mesh = load_ply(filename)
    verts = vcat(Float64.(mesh["vertex"]["x"]).', Float64.(mesh["vertex"]["y"]).', Float64.(mesh["vertex"]["z"]).')
    faces = hcat(collect(Int.(y).+1 for y in mesh["face"]["vertex_indices"])...)
    (verts,faces)
end

verts,faces = readply("chopper.ply")
tree = AABBTree(verts,faces)

function imgit(tree, verts)
    imgsz = 512;
    bmin = minimum(verts[[2,3],:],2)
    bmax = maximum(verts[[2,3],:],2)
    mindist = 1e9;
    maxdist = -1e9;

    y = 1:imgsz;
    x = 1:imgsz;

    # view from the side
    fracy = (imgsz-y)./imgsz;
    cy = (bmax[1] - bmin[1]) .* fracy .+ bmin[1];
    fracz = x./imgsz;
    cz = (bmax[2] - bmin[2]) .* fracz .+ bmin[2];

    fy = repmat(cy.', length(cz), 1)
    fz = repmat(cz, 1, length(cy))
    #[fz,fy] = meshgrid(cy,cz); christ no meshgrid??

    from = vcat(-1000*ones(size(fy[:]))', fy[:]', fz[:]');
    to = vcat(1000*ones(size(fy[:]))', fy[:]', fz[:]');

    hits,dists,faceinds,barys = raycast(tree, from, to, mode=1)
    img = reshape(dists,size(fz))
end

img = imgit(tree,verts)
imwritepnm(img, "x1.ppm")

verts .= verts .+ rand(size(verts)).*20
deform!(tree, verts)

img = imgit(tree,verts)
imwritepnm(clamp.(img,0,1), "x2.ppm")

if false
verts = rand(3,20); faces = hcat(randperm(size(verts,2)), randperm(size(verts,2)),randperm(size(verts,2)))';
a = Opcode.aabb_tree_create(verts,size(verts,2),faces,size(faces,2));
origins = rand(3,10); dirs = rand(3,10); hits = Vector{Int}(); dists = Vector{Float64}(); faceinds = Vector{Int}(); barys = Vector{Float64}(0);
Opcode.aabb_tree_intersect(a, origins, dirs, hits, dists, faceinds, barys, size(origins,2))
end


text:
using DifferentialEquations
using ITensors
using ITensorTDVP

function ode_solver(
  H::TimeDependentSum,
  time_step,
  ψ₀;
  current_time=0.0,
  outputlevel=0,
  solver_alg=Tsit5(),
  kwargs...,
)
  if outputlevel ≥ 3
    println("    In ODE solver, current_time = $current_time, time_step = $time_step")
  end

  time_span = (current_time, current_time + time_step)
  u₀, ITensor_from_vec = to_vec(ψ₀)
  f(ψ::ITensor, p, t) = H(t)(ψ)
  f(u::Vector, p, t) = to_vec(f(ITensor_from_vec(u), p, t))[1]
  prob = ODEProblem(f, u₀, time_span)
  sol = solve(prob, solver_alg; kwargs...)
  uₜ = sol.u[end]
  return ITensor_from_vec(uₜ), nothing
end

function krylov_solver(
  H::TimeDependentSum, time_step, ψ₀; current_time=0.0, outputlevel=0, kwargs...
)
  if outputlevel ≥ 3
    println("    In Krylov solver, current_time = $current_time, time_step = $time_step")
  end
  ψₜ, info = exponentiate(H(current_time), time_step, ψ₀; kwargs...)
  return ψₜ, info
end


text:
@testset "Lagrangian lower bound tests" begin
    begin
        α = 0.5
        l = 2 * (1 - α)
        @test relaxedIncrementalProblem([2 3; 1 4], [2 2; 3 1], [5 4; 3 4], 6, getAssignmentConstraints(2), 0.5, l) == (11, :Optimal)
    end

    begin
        α = 0.5
        l = 2 * (1 - α)
        @test lagrangianLowerBound([2 3; 1 4], [2 2; 3 1], [5 4; 3 4], 6, getAssignmentConstraints(2), l, AssignmentProblemDescriptor(2)) == 11
    end
end


text:
include("../src/Graph.jl")
include("../src/Tools.jl")
include("../src/Algorithm.jl")
include("../src/Sampling.jl")
include("../src/Greedy.jl")

using StatsBase
using JSON
using JLD2

"""
Set Seed Nodes & change innate opinions
Set parameters.
k = 50
epsilon : s + epsilon
delta : edge probability parameter
beta: threshold
json struct js
"""

epsilon = 0.1
delta = 0.5
beta = 0.5
tags = []

jsdict = Dict{String, Any}()

# Read Graph
buf = split(ARGS[1], ',')
fileName = string("../data/all/", buf[1], ".jld2")

println(fileName)

(name, n, m, V, E) = load(fileName, "name", "n", "m", "V", "E")
G = Graph(n, m, V, E)

merge!(jsdict, Dict("name"=>name, "n"=>n, "m"=>m))

if ARGS[2] == "1"
    #Weighted Cascade
    push!(tags, "WC")
    E_p = load(fileName, "E_WC")
elseif ARGS[2] == "2"
    #Trivalency
    push!(tags, "TM") # Trivalency Model
    E_p = load(fileName, "E_TM")
end
G2 = Graph(n, m, V, E_p)

# Before Making changes
if ARGS[3] == "1"
    # Uniform distribution
    push!(tags, "Uniform")
    s = load(fileName, "s_Uni")
elseif ARGS[3] == "2"
    # Exponential distribution
    push!(tags, "Exponential")
    s = load(fileName, "s_Exp")
elseif ARGS[3] == "3"
    # Power-law distribution
    push!(tags, "Power-law")
    s = load(fileName, "s_Pow")
end

# Different parameters to select k
if ARGS[4] == "1"
    #push!(tags, "0.005")
    k = 1
elseif ARGS[4] == "3"
    #push!(tags, "0.015")
    k = 3
elseif ARGS[4] == "5"
    #push!(tags, "0.02")
    k = 5
end

println("n is ", n)
println("m is ", m)
println("k is ", k)

merge!(jsdict, Dict("k"=>k))
merge!(jsdict, Dict("model"=>tags[1], "s_dist"=>tags[2]))

name_greedy = ["aci", "ad", "ap", "aidc"]
name = ["aci", "ad", "ap", "aidc", "aidc"]

##################### Above are configs
#####################
# The only difference for marketing setting and backfire setting is that Δs is different

if ARGS[5] == "1"
    Δs = delta_s_marketing(s, epsilon)
elseif ARGS[5] == "2"
    Δs = delta_s_backfire(s, epsilon, beta)
end

##################### Below are related to algorithms
######################

#println("s is ", Δs)
Δaci, Δad, Δap, Δaidc, Δaidc_upp = assign_node_weights(G, s, Δs) # linear part weight
lin_qua = [Δaci, Δad, Δap, Δaidc]
lin = [Δaci, Δad, Δap, Δaidc, Δaidc_upp]

Δci, Δd, Δp, Δidc = assign_pair_weights(G, Δs) # quadratic part weight
qua = [Δci, Δd, Δp, Δidc]

lin_labels = ["max_aci_", "max_ad_", "max_ap_", "max_aidc_", "max_aidc_upp_"] # max lin part
qua_labels = ["greedy_aci_", "greedy_ad_", "greedy_ap_", "greedy_aidc_"]

# The initial of different measures
_, init_aci, init_ad, init_ap, init_aidc = Approx(G, s)
init_result_greedy = [init_aci, init_ad, init_ap, init_aidc, init_aidc]
#
init_labels = "init_".*name_greedy
for i in 1:length(name_greedy)
   merge!(jsdict, Dict(init_labels[i]=>init_result_greedy[i]))
end

#println("init aci is ", init_aci)

######################### For the pairs

ris_me, T_me, seeds_me, R_me, R_sample = [], [], [], [], []

for i in 1:length(qua)
    (ris, seeds_qua, R_size, T_qua, R) = ris_measure_greedy(G2, k, lin_qua[i], qua[i], delta) ## RR to select seed nodes
    push!(ris_me, ris)
    push!(T_me, T_qua)
    push!(R_me, R_size)
    push!(R_sample, R)
    push!(seeds_me, seeds_qua)
end

for i in 1:length(qua)
    nn = qua_labels[i].*name_greedy
    lin_nn = qua_labels[i]*"rr" ## init rr_algorithm
    lin_t = qua_labels[i]*"time"
    lin_size = qua_labels[i]*"r_size"

    merge!(jsdict, Dict(nn[i]=>(ris_me[i] + init_result_greedy[i])))
    merge!(jsdict, Dict(lin_t=>T_me[i]))
    merge!(jsdict, Dict(lin_size=>R_me[i]))
end

################### For maximize linear part

ris_me, seeds_me, T_me, R_me = [], [], [], []
# Calculation
for i = 1:length(lin)
    flag = 1
    (ris, seeds_lin, R_size, T_lin) = ris_measure_linear(G2, k, lin[i], flag, delta) ## RR to select seed nodes
    push!(ris_me, ris)
    push!(T_me, T_lin)
    push!(R_me, R_size)
    push!(seeds_me, seeds_lin)
end

for i = 1:length(lin)
    nn = lin_labels[i].*name
    lin_nn = lin_labels[i]*"lin" # only the linear part
    lin_t = lin_labels[i]*"time"
    lin_size = lin_labels[i]*"r_size"

    if i < 5
        result = evaluation_ris(R_sample[i], n, lin_qua[i], qua[i], seeds_me[i])
    else
        result = evaluation_ris(R_sample[4], n, lin_qua[4], qua[4], seeds_me[i])
    end

    merge!(jsdict, Dict(nn[i]=>(result + init_result_greedy[i])))
    merge!(jsdict, Dict(lin_t=>T_me[i]))
    merge!(jsdict, Dict(lin_size=>R_me[i]))
    merge!(jsdict, Dict(lin_nn=>ris_me[i] + init_result_greedy[i]))
end

####################### IO

js = JSON.json(jsdict)

println(jsdict)

if ARGS[5] == "1"
    open("out/greedy_marketing.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
elseif ARGS[5] == "2"
    open("out/greedy_backfire.json", "a+") do ff
        write(ff, js)
        write(ff, "\n")
    end
end


text:
"This is just a wrapper"
struct FunctionalPolicy{Ta <: Function} <: AbstractPolicy
    action_generator::Ta
end

(p::FunctionalPolicy)(s) = p.action_generator(s)


text:
# arithmetics.jl
# Routines having to do with computations involving domains

+(domain::Domain) = domain

function +(domain::Domain{T}, x::S) where {T,S}
    c = convert_space(spacetype(T), x)
    c isa S && error("Cannot add $x to $domain")
    (+)(domain, c)
end

*(map::AbstractMap, domain::Domain) = forwardmap_domain(map, domain)

*(a::Number, domain::Domain{T}) where {T} = LinearMap{T}(a) * domain
*(domain::Domain, a::Number) = a*domain

/(domain::Domain{T}, a::Number) where {T} = LinearMap{T}(inv(a)) * domain

+(d::Domain, x::SVector{N,T}) where {N,T} = Translation(x) * d

# Assume commutativity
+(x::AbstractVector, d::Domain) = d+x

# Rotation around the origin
rotate(d::EuclideanDomain{2}, θ) = rotation_map(θ) * d

rotate(d::EuclideanDomain{3}, phi, theta, psi) = rotation_map(phi,theta,psi) * d
# Rotation around a fixed center.
rotate(d::EuclideanDomain{2}, θ, center::SVector{T}) where {T} = Translation(center) * (rotation_map(θ) * (Translation(-center) * d))

rotate(d::EuclideanDomain{3}, phi, theta, psi, center::SVector{T}) where {T} = Translation(center) * (rotation_map(phi,theta,psi) * (Translation(-center) * d))


text:
export CyclicRange

#= make CyclicRange For ScaLAPACK =#

struct CyclicRange{T <: Integer} <: AbstractRange{T}
    start::T # full array's first
    stop::T # full array's last
    blocksize::T # blocksize
    n_split::T # number of process
    split_index::T # what process

    CyclicRange(;start::T, stop::T, blocksize::T, n_split::T, split_index::T) where T <: Integer = new{T}(start, stop, blocksize, n_split, split_index)
end

number_element(r::CyclicRange) = r.stop - r.start + 1

function number_element_local(r::CyclicRange)
    n_element = number_element(r)
    n_block_local_min = number_block_local_min(r)
    n_extra_block = number_block_global_extra(r)

    n_element_local = r.blocksize*n_block_local_min
    n_extra_element = n_element % r.blocksize

    # add extra block
    if r.split_index <= n_extra_block
        n_element_local += r.blocksize
    end
    # add extra element,
    if r.split_index == n_extra_block + 1 && n_extra_element != 0
        n_element_local += n_extra_element
    end
    return n_element_local
end

function number_block(r::CyclicRange)
    n_element = number_element(r)
    n_block = n_element ÷ r.blocksize
    return n_block
end

function number_block_local_min(r::CyclicRange{T}) where T <: Integer
    n_block = number_block(r)
    n_block_local = n_block ÷ r.n_split
    return n_block_local
end

function number_block_global_extra(r::CyclicRange{T}) where T <: Integer
    n_block = number_block(r)
    n_extra_block = n_block % r.n_split
    return n_extra_block
end

function local_index(r::CyclicRange, i::Integer)
    i_temp = i - 1
    block_position = i_temp ÷ r.blocksize
    local_index = (block_position ÷ r.n_split)*r.blocksize + (i_temp % r.blocksize) + 1
    split_index = block_position % r.n_split + 1
    # mpi rank = split_index - 1, in 0-based local index = extra_element - 1
    return split_index, local_index
end

Base.firstindex(::CyclicRange) = 1

Base.first(r::CyclicRange) = getindex(r, 1)

Base.last(r::CyclicRange) = getindex(r, length(r))

Base.length(r::CyclicRange{T}) where T <: Integer = number_element_local(r)

Base.getindex(r::CyclicRange, i::Integer) = r.blocksize*(r.split_index-1) + r.n_split*r.blocksize*div(i-1, r.blocksize) + mod(i-1, r.blocksize) + 1

Base.show(io::IO, r::CyclicRange) = print(io, "$(r.start):[$(r.split_index) $(r.blocksize) $(r.n_split)]:$(r.stop)")


text:
@discord_object struct UserNickChange
    nick::String
end


text:
function ground(args...)
    background("white")
    sethue("black")
end

@testset "Javis Viewer" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    render(vid; pathname = "")

    action_list = [back, star_obj]

    viewer_win, frame_dims, r_slide, tbox, canvas, actions, total_frames, video =
        Javis._javis_viewer(vid, 100, action_list, false)
    visible(viewer_win, false)

    @test get_gtk_property(viewer_win, :title, String) == "Javis Viewer"

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    second_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 2

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 1

    @test first_frame != second_frame

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    last_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == total_frames

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == 1

    @test last_frame != first_frame
end


text:
using LogisticOptTools
using JSOSolvers, NLPModels

const LOT = LogisticOptTools

mutable struct LogReg <: AbstractNLPModel
    meta::NLPModelMeta
    counters::Counters
    logreg::LOT.LogisticRegressor
end

function LogReg(logreg)
    p = LOT.nfeatures(logreg)
    meta = NLPModelMeta(p, x0=zeros(p),
                        name="LogisticRegression")
    return LogReg(meta, Counters(), logreg)
end

function NLPModels.obj(nlp :: LogReg, x :: AbstractVector)
    return LOT.loss(x, nlp.logreg)
end

function NLPModels.grad!(nlp :: LogReg, x :: AbstractVector, gx :: AbstractVector)
    LOT.gradient!(gx, x, nlp.logreg)
    return gx
end

function NLPModels.hprod!(nlp :: LogReg, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight=1.0)
    LOT.hessvec!(Hv, x, v, nlp.logreg)
    return Hv
end

DATASET = "covtype.libsvm.binary.bz2"
if !isdefined(Main, :dataset)
    #= dataset = LOT.SparseLogitData(DATASET, scale_data=true) =#
    # Non-scaling data yields more difficult optimization problems
    dataset = LOT.LogitData(DATASET, scale_data=false)
end

penalty = LOT.L2Penalty(1e-2)
logreg = LOT.LogisticRegressor(dataset.X, dataset.y, penalty=penalty)

nlp = LogReg(logreg)
results = @time tron(nlp, atol=1e-6, max_time=300.0)


text:
# Integration tests on some simple test cases.

using Test
using StaticArrays, FiniteVolumes
using FiniteVolumes: UnidirectionalFlux
using FiniteVolumes.CartesianMeshes: Half

# TOOLS
getfield(w, var) = [wi[var] for wi in w]

const ϵ = 1e-15
boundedness(w, lb, ub) = all(lb - ϵ .≤ w .≤ ub + ϵ)
boundedness(w, lb, ub, var) = all(lb - ϵ .≤ getfield(w, var) .≤ ub + ϵ)
maximum_principle(w, w₀) = boundedness(w, minimum(w₀), maximum(w₀))
maximum_principle(w, w₀, var) = maximum_principle(getfield(w, var), getfield(w₀, var))

riemann_problem(mesh, w₁, w₂, step_position=0.5) = [x[1] < step_position ? w₁ : w₂ for x in cell_centers(mesh)]

@testset "Cases" begin

@testset "Scalar problems" begin
    @testset "Trivial cases" begin
        mesh = CartesianMesh(10)

        # No velocity => no change
        w₀ = rand(FiniteVolumes.nb_cells(mesh))
        flux = LinearAdvectionFlux(0.0)
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=0.001, nb_time_steps=3, verbose=false)
        @test all(w .== w₀)

        # Uniform initial condition => no change
        w₀ = ones(FiniteVolumes.nb_cells(mesh))
        flux = LinearAdvectionFlux(1.0)
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=0.001, nb_time_steps=3, verbose=false)
        @test all(w .== w₀)
    end

    @testset "1D linear advection" begin
        mesh = CartesianMesh(100)

        sine_w₀ = map(x -> sin(x[1]), cell_centers(mesh))
        falling_step_w₀ = riemann_problem(mesh, 1.0, 0.0)
        rising_step_w₀ = riemann_problem(mesh, 0.0, 1.0)
        initial_conditions = [sine_w₀, falling_step_w₀, rising_step_w₀]

        forward_flux = LinearAdvectionFlux(1.0)
        backward_flux = LinearAdvectionFlux(-1.0)

        settings = (time_step=FixedCourant(0.1), nb_time_steps=5, verbose=false, numerical_flux=Upwind())

        for w₀ in initial_conditions
            # Test stability
            t, wf = FiniteVolumes.run(forward_flux, mesh, w₀; settings...)
            @test maximum_principle(wf, w₀)

            # Test left-right symmetry
            t, wb = FiniteVolumes.run(backward_flux, mesh, reverse(w₀); settings...)
            @test all(wf .== reverse(wb))

        end

        t, w_unstable = FiniteVolumes.run(forward_flux, mesh, sine_w₀; time_step=FixedCourant(1.5), nb_time_steps=5, verbose=false)
        @test !maximum_principle(w_unstable, sine_w₀)
    end

    @testset "1D Wave equation" begin
        mesh = CartesianMesh(10)
        w₀ = [SVector{2}(sin(2π*x), cos(2π*x)) for x in cell_centers(mesh)]

        # Anonymous model
        f(α, n) = n * SVector(α[2], α[1])
        t, w₁ = FiniteVolumes.run(f, mesh, w₀, time_step=FixedCourant(0.2), nb_time_steps=10, verbose=false)

        # Fancy model
        struct Wave1DFlux{T}
            velocity::T
        end
        (wf::Wave1DFlux)(w, n) = wf.velocity*SVector{2}(w[2], w[1])*n
        jacobian(f::Wave1DFlux, w, n) = f.velocity*SMatrix{2, 2}(0.0, 1.0, 1.0, 0.0)*n
        FiniteVolumes.eigvals(f::Wave1DFlux, w, n) = SVector{2}(-1.0, 1.0)*n
        FiniteVolumes.eigen(f::Wave1DFlux, w, n) = (SVector{2}(-1.0, 1.0)*n, SMatrix{2, 2}(-0.707107, 0.707107, 0.707107, 0.707107)*n)
        t, w₂ = FiniteVolumes.run(Wave1DFlux(1.0), mesh, w₀, time_step=FixedCourant(0.2), nb_time_steps=10, verbose=false)

        @test all(isapprox.(w₁, w₂, atol=1e-6))
    end

    @testset "1D Burger" begin
        mesh = CartesianMesh(100)
        flux(u, n) = 0.5*u.^2 * n
        step(x) = 0.33 < x[1] < 0.66 ? 1.0 : 0.0
        w₀ = map(step, cell_centers(mesh))
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=FixedCourant(0.2), nb_time_steps=10, verbose=false)
        @test maximum_principle(w, w₀)
    end

    @testset "Time dependant advection" begin
        mesh = CartesianMesh(100)
        flux = FiniteVolumes.TimeDependantFlux(t -> LinearAdvectionFlux(sin(t)))
        step(x) = 0.33 < x[1] < 0.66 ? 1.0 : 0.0
        w₀ = map(step, cell_centers(mesh))
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=0.005, nb_time_steps=10, verbose=false)
        @test maximum_principle(w, w₀)
    end

    @testset "Splitted 2D linear advection" begin
        mesh = PeriodicCartesianMesh(20, 20)
        splitted_advection(u₀) = FiniteVolumes.directional_splitting(LinearAdvectionFlux(u₀))

        @test splitted_advection(SVector(1.0, 1.0)) == (UnidirectionalFlux{1}(LinearAdvectionFlux(SVector(1.0, 0.0))), UnidirectionalFlux{2}(LinearAdvectionFlux(SVector(0.0, 1.0))))

        settings = (time_step=FixedCourant(0.1), nb_time_steps=5, numerical_flux=Upwind(), verbose=false)

        # Scalar
        square(x, side=0.5) = all(0.5-side/2 .<= x .<= 0.5+side/2) ? 1.0 : 0.0
        w₀ = map(square, cell_centers(mesh))
        t, w = FiniteVolumes.run(splitted_advection(SVector(1.0, 1.0)), mesh, w₀; settings...)
        @test maximum_principle(w, w₀)

        # Multi-scalar
        function triple_point(x)
            if x[1] < 0.5
                return SVector(1.0, 0.0, 0.0)
            elseif x[2] < 0.5
                return SVector(0.0, 1.0, 0.0)
            else
                return SVector(0.0, 0.0, 1.0)
            end
        end
        w₀ = map(triple_point, cell_centers(mesh))
        t, w = FiniteVolumes.run(splitted_advection(SVector(1.0, 1.0)), mesh, w₀; settings...)
        @test maximum_principle(w, w₀, 1)
        @test maximum_principle(w, w₀, 2)
        @test maximum_principle(w, w₀, 3)
        @test all(sum.(w) .≈ 1.0)
    end

    @testset "2D rotation" begin
        mesh = PeriodicCartesianMesh(40, 40)
        flux = RotationFlux(SVector(0.5, 0.5))

        square(x, side=0.5) = all(0.5-side/2 .<= x .<= 0.5+side/2) ? 1.0 : 0.0
        w₀ = map(square, cell_centers(mesh))

        settings = (time_step=FixedCourant(0.1), nb_time_steps=10, numerical_flux=Upwind(), verbose=false)
        t, w = FiniteVolumes.run(flux, mesh, w₀; settings...)
        @test maximum_principle(w, w₀)

        t, w = FiniteVolumes.run(directional_splitting(flux), mesh, w₀; settings...)
        @test maximum_principle(w, w₀)
    end

    @testset "2D diagonal Burger" begin
        #   ∂_t α + u₀ ⋅ ∇ α^2/2 = 0
        # ⟺ ∂_t α + div(α^2/2 * u₀) = 0
        # ≡ rotated 1D Burger
        mesh = PeriodicCartesianMesh(20, 20)
        flux(α, n) = 0.5*α.^2 * n' * [1.0, 1.0]
        vertical_band(x) = 0.33 < x[1] < 0.66 ? 1.0 : 0.0
        w₀ = map(vertical_band, cell_centers(mesh))
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=FixedCourant(0.2), nb_time_steps=10, verbose=false)
        @test maximum_principle(w, w₀)
    end

    @testset "Splitted 3D linear advection" begin
        mesh = PeriodicCartesianMesh(20, 20, 20)
        splitted_advection(u₀) = FiniteVolumes.directional_splitting(LinearAdvectionFlux(u₀))

        settings = (time_step=FixedCourant(0.1), nb_time_steps=5, numerical_flux=Upwind(), verbose=false)

        # Scalar
        cube(x, side=0.5) = all(0.5-side/2 .<= x .<= 0.5+side/2) ? 1.0 : 0.0
        w₀ = map(cube, cell_centers(mesh))
        t, w = FiniteVolumes.run(splitted_advection([1.0, 1.0, 1.0]), mesh, w₀; settings...)
        @test maximum_principle(w, w₀)
    end
end

@testset "Shallow water" begin

    total_mass(mesh, v) = sum(v[i][1] * FiniteVolumes.cell_volume(mesh, i) for i in FiniteVolumes.all_cells(mesh))
    total_momentum(mesh, v) = sum(v[i][2] * FiniteVolumes.cell_volume(mesh, i) for i in FiniteVolumes.all_cells(mesh))

    @testset "1D" begin
        flux = ShallowWater(9.8)
        @test Upwind()(flux, CartesianMesh(2), [SVector(1.0, 0.0), SVector(1.0, 0.0)], (Half(3,))) == SVector(0.0, 9.8/2)
        @test Upwind()(flux, CartesianMesh(2), [SVector(2.0, 0.0), SVector(1.0, 0.0)], (Half(3,))) ≈ SVector(1.9170289512680811, 12.25)

        mesh = CartesianMesh(100)
        w₀ = riemann_problem(mesh, SVector(2.0, 0.0), SVector(1.0, 0.0))
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=FixedCourant(0.4),
                                 nb_time_steps=50, verbose=false)

        w₀ = map(x -> SVector(1.0 + exp(-500*(x[1]-0.5)^2), 0.0), cell_centers(mesh))
        t, w = FiniteVolumes.run(flux, mesh, w₀, time_step=FixedCourant(0.4),
                                 nb_time_steps=50, verbose=false)

        @test total_mass(mesh, w₀) .≈ total_mass(mesh, w)
        @test total_momentum(mesh, w₀) .≈ total_momentum(mesh, w) atol=1e-15
        # plot(mesh, w, 1); plot!(mesh, w, 2)
    end

    @testset "2D" begin
        mesh = PeriodicCartesianMesh(40, 40)
        flux = ShallowWater(9.8)

        v₀ = map(x -> SVector(1.0 + exp(-500*((x[1]-0.5)^2 + (x[2]-0.5)^2)), 0.0, 0.0), cell_centers(mesh))
        FiniteVolumes.div(flux, mesh, v₀, Upwind())
        t, v = FiniteVolumes.run(flux, mesh, v₀, time_step=FixedCourant(0.3),
                                 nb_time_steps=40, verbose=false)

        @test total_mass(mesh, v₀) .≈ total_mass(mesh, v)
        @test total_momentum(mesh, v₀) .≈ total_momentum(mesh, v) atol=1e-15
        # plot(mesh, v, 1)
    end

end

end


text:
include("test_CG.jl")
include("test_CG_memory_vs_disk.jl")
include("test_FFTOp.jl")
include("test_MWNI.jl")
include("test_POCS.jl")
include("test_SeisLinearEvents.jl")
#include("test_AcousticWave2D.jl")


text:
# HTML 作成
module HtmlTable
using DataFrames
using CSV
using Tables

export fromDataFrame, fromMatrix, fromCSV

# HTML 文書のテンプレート
const HTMLFMT = """<!DOCTYPE html>

  <meta name="viewport" content="width=device-width,initial-scale=1" />
  <meta charset="utf-8" />
  <title>%title</title>

%title

    %table

"""

# スタイル
css = """\n
body {
    margin-left:3%;
    margin-right:3%;
}

h1 {
    text-align:center;
    padding:10px;
}

hr {
    margin-bottom:20px;
}

table {
    width:100%;
    border: solid thin gray;
    color: dimgray;
    border-collapse:collapse;
    margin-top:30px;
    margin-bottom:40px;
}

th {
    border: solid thin gray;
    background-color: azure;
    padding:3px;
    border-collapse:collapse;
}

td {
    border: solid thin gray;
    font-size:small;
    padding:3px;
    border-collapse:collapse;
}\n
"""

# HTML の 特殊文字をエスケープするかどうかのフラグ
flagEscape = true

# HTML の 特殊文字をエスケープする。
function htmlEscape(text::String)
    return replace(text, "&"=>"&", "<"=>"<", ">"=>">")
end

# HTML タグを作る。
function tag(tagname::String, text::Any; escape=true, style="")
    tagn = tagname
    if style != ""
        tagn *= " "
        tagn *= style
    end
    if escape
        return "<" * tagn * ">" * htmlEscape(string(text)) * "</" * tagname * ">"
    else
        return "<" * tagn * ">" * text * "</" * tagname * ">"
    end
end

# テーブル行を作る。
function tableRow(data::Vector; th=false, rown=-1)
    td = th ? "th" : "td"
    html = "
"
    # 行番号を付けるかどうか
    if rown == 0
        html *= tag(td, "Row", style="style=\"width:100px;\"")  # 表題行の場合
    elseif rown > 0
        html *= tag(td, rown)  # 表題行でない場合
    else
        # rown < 0 のときは行番号なし。
    end
    for d = data
        html *= tag(td, d, escape=flagEscape)
    end
    html *= "
\n"
    return html
end

# Matrix を元に作成する。
function fromMatrix(mat::Matrix, columns::Vector; title="HTML Table", rownum=false)
    table = "
\n
\n"
    table *= tableRow(columns, th=true, rown=rownum ? 0 : -1)
    table *= "
\n
\n"
    nrow, ncol = size(mat)
    for i = 1:nrow
        table *= tableRow(mat[i, :], rown=rownum ? i : -1)
    end
    table *= "
\n
\n"
    replace(HTMLFMT, "%title" => title, "%table" => table, "%style" =>css)
end

# DataFrame を元に作成する。
function fromDataFrame(df::DataFrame; title="HTML Table", rownum=false)
    columns = names(df)
    mat = Tables.matrix(df)
    fromMatrix(mat, columns, title=title, rownum=rownum)
end

# CSV ファイルを元に作成する。
function fromCSV(path::String; header=true, delimtter=",", stripQuote=true, title="HTML Table", rownum=false)
    df = DataFrame(CSV.File(path, header=header, quoted=stripQuote))
    fromDataFrame(df, title=title, rownum=rownum)
end

end # HtmlTable

max_stars_repo_name: makandat/julia

max_issues_repo_name: makandat/julia

max_forks_repo_name: makandat/julia


text:
immutable Remove end

# Recursively walk an Expr, applying a function to each subexpr and leaf to
# build up a new Expr. If the function returns an Expr, that Expr will be
# walked as well. If the function returns the singleton value
# ExpressionUtils.Remove, that node will be removed.
#
walk(f, leaf) = f(leaf)
function walk(f, ex::Expr)
    ex = copy(ex)

    function reducer(args, e)
        e = f(e)
        is(e, Remove) ? args :
        isa(e, Expr)    ? push!(args, walk(f, e)) :
        push!(args, e)
    end

    ex.args = reduce(reducer, Any[], ex.args)
    ex
end

getesc(d::Dict, k, default) = haskey(d, k) ? esc(d[k]) : default

Base.map(f, ex::Expr) =
    Expr(ex.head, [isa(e, Expr) ? map(f, e) : f(e) for e in ex.args]...)

remove_quote_block(val) = val
remove_quote_block(ex::Expr) =
    ex.head == :block && length(ex.args) == 2 ? ex.args[2] : ex

function remove_line_ann(ex::Expr)
    remove_line_nodes(node::LineNumberNode) = Remove
    remove_line_nodes(e) = e
    walk(remove_line_nodes, ex)
end

symbeginswith(sym, s) = startswith(string(sym), s)
issplat(sym)   = isa(sym, Symbol) && symbeginswith(sym, "_SPLAT_")
isunsplat(sym) = isa(sym, Symbol) && symbeginswith(sym, "_UNSPLAT_")
isesc(sym)     = isa(sym, Symbol) && symbeginswith(sym, "_ESC_")
function istemplate(sym::Symbol)
    s = string(sym)
    s[1] == '_' && s[end] == '_'
end

removeprefix(sym) =
    if issplat(sym)
        Symbol(string(sym)[7:end])
    elseif isunsplat(sym)
        Symbol(string(sym)[9:end])
    elseif isesc(sym)
        Symbol(string(sym)[5:end])
    else
        sym
    end

function expr_bindings(ex, template::Symbol, collected::Dict)
    if isesc(template)
        collected[removeprefix(template)] = esc(ex)
    elseif istemplate(template)
        collected[template] = ex
    end
    collected
end
function expr_bindings(ex::Expr, template::Expr, collected::Dict)
    ex.head == template.head || error("Cannot bind dissimilar syntax:\n$ex\nto\n$template")
    for (i, (exsub, templatesub)) in enumerate(zip(ex.args, template.args))
        if issplat(templatesub)
            collected[removeprefix(templatesub)] = ex.args[i:end]
            break
        end
        expr_bindings(exsub, templatesub, collected)
    end
    collected
end
expr_bindings(ex, template, collected) = collected
expr_bindings(ex, template) =
    expr_bindings(ex, template, Dict{Symbol, Any}())

function unsplat(ex::Expr, bindings::Dict)
    r = (args, e) -> isunsplat(e) ?
                     vcat(args, get(bindings, removeprefix(e), [e])) :
                     push!(args, get(bindings, removeprefix(e), e))

    ex = copy(ex)
    ex.args = reduce(r, [], ex.args)
    ex
end

function expr_bind(ex::Expr, bindings::Dict)
    walker(e::Expr) = unsplat(e, bindings)
    walker(s::Symbol) = get(bindings, removeprefix(s), s)
    walker(e) = e

    ex = copy(ex)

    walk(walker, ex)
end

function expr_replace(ex, template, out)
    ex, template, out = map((e) -> remove_line_ann(remove_quote_block(e)),
                            [ex, template, out])

    bindings = expr_bindings(ex, template)
    expr_bind(remove_quote_block(out), bindings)
end


text:
function inner_blas2(x::Array{Float32,1}, y::Array{Float32,1}, dist::Array{Float32,1}, len::UInt64)
    @fastmath @inbounds for i = 1:len
        #println("left = ", (triang(i-1)+1), "\tright = ", triang(i))
        fragment!(view(dist, (triang(i-1)+1):triang(i)), view(x, 1:i), view(y, 1:i), x[i], y[i])
    end
end

function blas_stat_st2(indx, step::Float64, maxrng)
    F = Tuple.(indx)
    x, y= Float32.(first.(F)), Float32.(last.(F))
    dist = Vector{Float32}(undef,triang(length(x)))
    len::UInt64 = length(x)
    inner_blas2(x, y, dist, len)
    #print(length(fit(Histogram, dist, 0:step:(maxrng)).weights))
    return fit(Histogram, dist, 0:step:(maxrng)).weights
end

function fragment!(dist, arrx, arry, xi, yi)
    dist .= sqrt.((arrx .- xi) .^2 .+ (arry .- yi) .^2)
end

q = rand(Bool, 200,100)

v = findall(x->x==true, q)

function cart2arr(v)
    x, y = Float32.(first.(Tuple.(v))), Float32.(last.(Tuple.(v)))
end
function test(v)
    for i in 1:10
        x, y = Float32.(first.(Tuple.(v))), Float32.(last.(Tuple.(v)))
    end
end

#@btime ReconTPCF.blas_stat_st2(v, 1.0, 100)
x,y = cart2arr(v)

m = 3
n = 8
nx = 10000
ny = 8

a = rand(Float64, m, nx)
b = rand(Float64, m, ny)

#@btime fit(Histogram, vec(pairwise(Euclidean(), a, dims=2)))

using NearestNeighbors

using ReconTPCF

data = rand(1:1.0:100, 2, 10^4)
kdtree = KDTree(data; leafsize = 4)
NearestNeighbors.show(kdtree)

kdtree.hyper_rec


text:
#=
Given a stream of elements too large to store in memory, pick a random element from the stream with uniform probability.
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module Akagera

using JSON

include("animecontainers.jl")

# I need exportall
export Animator,
       init!,
       set_msg!,
       start!,
       update_msg!,
       finish!,

       AnimeContainer,
       dots,
       dots2,
       dots3,
       dots4,
       dots5,
       dots6,
       dots7,
       dots8,
       dots9,
       dots10,
       dots11,
       dots12,
       line,
       line2,
       pipe,
       simpleDots,
       simpleDotsScrolling,
       star,
       star2,
       flip,
       hamburger,
       growVertical,
       growHorizontal,
       balloon,
       balloon2,
       noise,
       bounce,
       boxBounce,
       boxBounce2,
       triangle,
       arc,
       circle,
       squareCorners,
       circleQuarters,
       circleHalves,
       squish,
       toggle,
       toggle2,
       toggle3,
       toggle4,
       toggle5,
       toggle6,
       toggle7,
       toggle8,
       toggle9,
       toggle10,
       toggle11,
       toggle12,
       toggle13,
       arrow,
       arrow2,
       arrow3,
       bouncingBar,
       bouncingBall,
       smiley,
       monkey,
       hearts,
       clock,
       earth,
       moon,
       runner,
       pong,
       shark,
       dqpb,
       weather,
       christmas

"""
       Animator(interval::Int, frames::Vector{String}, animate_type::String)
       Animator(json::String; animate_type="")
       Animator(fc::AnimeContainer)

Construct an `Animator` with animation loaded from given arguments.
"""
mutable struct Animator
    flag::Bool
    interval::Int
    msg::String
    frames::Vector{String}
    animate_type::String

    function Animator(interval::Int, frames::Vector{String}, animate_type::String)
        self = new()
        self.flag = false
        self.interval = interval
        self.frames = frames
        self.animate_type = animate_type

        return self
    end

    function Animator(json::String; animate_type="") # JSON
        json_dict = JSON.parse(json)

        self = new()
        self.flag = false
        self.interval = begin
            if haskey(json_dict, "interval")
                json_dict["interval"]
            else
                throw(ArgumentError("interval is not specified. You should set interval in json."))
            end
        end

        self.frames = begin
            if haskey(json_dict, "frames")
                json_dict["frames"]
            else
                throw(ArgumentError("frames is not specified. You should set frames in json."))
            end
        end

        self.animate_type = begin
            if animate_type != ""
                animate_type
            else
                warn("animate_type is not specified. Default value `linear` will be set.")
                "linear"
            end
        end

        return self
    end

    function Animator(fc::AnimeContainer)
        self = new()
        self.flag = false

        self.interval = fc.interval
        self.frames = fc.frames
        self.animate_type = fc.animate_type

        return self
    end
end

"""
    init!(a::Animator)

Set `a.flag = true` to be ready to start the animation. This should not be called manually.
"""
function init!(a::Animator)
    a.flag = true
end

"""
    set_msg!(a::Animator, msg::String)

Set `msg` field of `Animator`.

```julia-repl
julia> a = Animator(arrow3)
julia> set_msg!(a, "Hello, world!!!")
julia> a.msg
Hello, world!!!
```
"""
function set_msg!(a::Animator, msg::String)
    a.msg = msg
end

"""
    start!(a::Animator)

Start the animation.
"""
function start!(a::Animator; io=STDERR)
    init!(a)
    is_first = true

    function render(i::Int)
        if !is_first
            print(io, "\r\033[K")
        else
            println(io, " ") # Put empty line not to hurt old prints
            is_first = false
        end
        print(io, a.frames[i] * " " * a.msg)
        sleep(a.interval/1000)
    end

    @schedule while a.flag
        if a.animate_type == "linear"
            for i = 1:length(a.frames)
                # render i th element
                render(i)
            end
        elseif a.animate_type == "swing"
            for i = 1:2*length(a.frames) - 2
                # render ( i<=length(a.frames) ? i : (2*length(a.frames) - i) ) th element
                render(
                    begin
                        if i <= length(a.frames)
                            i
                        else
                            2*length(a.frames) - i
                        end
                    end)
            end
        end
    end
end

"""
    update_msg!(a::Animator, new_msg::String)

Change `msg` field of `Animator`. This works even while animation is activated.
"""
function update_msg!(a::Animator, new_msg::String)
    set_msg!(a, new_msg)
end

"""
    finish!(a::Animator)

Finish the animation.
"""
function finish!(a::Animator)
    a.flag = false
end

end # module


text:
using Koma

obj = read_phantom_jemris("koma/column1d_phantom.h5")
seq = read_seq("koma/epi_100x100_TE100_FOV230.seq")
sys = Scanner()
#Time
warmup = simulate(obj,seq[1:8],sys) #warmup, to precompile
signal = simulate(obj,seq,sys)


text:
export shaping_rnn,make_sequential_batch,get_rate_of_change,standardization!,decode_standardization!,normalization!,decode_normalization

#================
データ整形
=================#
function shaping_rnn(data,N)

    X,D = size(data,1),size(data,2)
    re_data = []

    #データの補充
    # if size(data,1)%N != 0
    #     #(n,D)分，paddingで埋める
    #     n = N - size(data,1)%N
    #     padding_arr = fill(padding,n,D)
    #     data = vcat(data,padding_arr) #補充
    #     X,D = size(data)
    # end

    # re_data = reshape(data,Int(X/N),N,D)

    #補充じゃなく，はみ出たデータは切り捨て
    re_data = reshape(data[1:X-X%N],Int(X÷N),N,D)
    re_data = permutedims(re_data,(2,1,3)) #軸の入れ替え

    return re_data
end

function make_sequential_batch(data,N,window_size)
    #T = window * 任意の数 + 1
    T_a = size(data,1) ÷ N
    n = (T_a - 1) ÷ window_size
    T = window_size * n + 1
    D = size(data,2)
    ed = N * T #はみ出したデータは削る

    re_data = reshape(data[1:ed,:],(T,N,D))
    return permutedims(re_data,(2,1,3))  #軸の入れ替え
end

function get_rate_of_change(array)
    rate_arr = zeros(length(array)-1)
    for i in 1:length(array)-1
        rate = (array[i+1]-array[i])/array[i]
        rate_arr[i] = rate
    end
    #行列化
    rate_arr = reshape(rate_arr,(length(rate_arr),1))
    return rate_arr
end

function standardization!(model::Models,data,name)
    if haskey(model.common.std_params,name)
        avg,variance,std = model.common.std_params[name]

        return (data.-avg)/std
    end

    avg = sum(data)/length(data)
    variance = sum((data.-avg).^2)/length(data)
    std = sqrt(variance)

    std_data = (data.-avg)/std

    model.common.std_params[name] = [avg,variance,std]

    return std_data
end
function decode_standardization(model::Models,data,name)
    avg,variance,std = model.common.std_params[name]

    return data .* std .+ avg
end

function normalization!(model::Models,data,name)
    if haskey(model.common.norm_params,name)
        max,min = model.common.norm_params[name]

        return (data.-min)./(max-min)
    end

    max = findmax(data)[1]
    min = findmin(data)[1]

    model.common.norm_params[name] = [max,min]

    if max == min return zeros(size(data)) end

    norm_data = (data.-min)./(max-min)

    return norm_data
end
function decode_normalization(model::Models,data,name)
    max,min = model.common.norm_params[name]

    return data .* (max-min) .+ min
end

max_stars_repo_name: catmanium/Spike

max_issues_repo_name: catmanium/Spike

max_forks_repo_name: catmanium/Spike


text:
#####
##### Compilation context
#####

# User-extended context allows parametrization of the pipeline through
# our subtype of AbstractInterpreter
abstract type CompilationContext end

struct NoContext <: CompilationContext end

@doc(
"""
    abstract type CompilationContext end

Parametrize the Mixtape pipeline by inheriting from `CompilationContext`. Similar to the context objects in [Cassette.jl](https://julia.mit.edu/Cassette.jl/stable/contextualdispatch.html). By using the interface methods [`transform`](@ref) and [`optimize!`](@ref) -- the user can control different parts of the compilation pipeline.
""", CompilationContext)

transform(ctx::CompilationContext, b) = b
transform(ctx::CompilationContext, b, sig) = transform(ctx, b)

@doc(
"""
    transform(ctx::CompilationContext, b::Core.CodeInfo)::Core.CodeInfo
    transform(ctx::CompilationContext, b::Core.CodeInfo, sig::Tuple)::Core.CodeInfo

User-defined transform which operates on lowered `Core.CodeInfo`. There's two versions: (1) ignores the signature of the current method body under consideration and (2) provides the signature as `sig`.

Transforms might typically follow a simple "swap" format using `CodeInfoTools.Builder`:

```julia
function transform(::MyCtx, src)
    b = CodeInfoTools.Builder(b)
    for (k, st) in b
        b[k] = swap(st))
    end
    return CodeInfoTools.finish(b)
end
```

but more advanced formats are possible. For further utilities, please see [CodeInfoTools.jl](https://github.com/JuliaCompilerPlugins/CodeInfoTools.jl).
""", transform)

optimize!(ctx::CompilationContext, b) = julia_passes!(b)

@doc(
"""
    optimize!(ctx::CompilationContext, b::OptimizationBundle)::Core.Compiler.IRCode

User-defined transform which operates on inferred IR provided by an [`OptimizationBundle`](@ref) instance.

The fallback implementation is:
```julia
optimize!(ctx::CompilationContext, b::OptimizationBundle) = julia_passes!(b)
```
which runs a set of standard (and required) Julia passes to the lowered and inferred `Core.Compiler.IRCode`.

!!! warning
    If you overload this method, _you are responsible for the optimization pass_! This means that you should know what you're doing (in general), and you will likely also want to call `julia_passes!` yourself. Be aware of this -- or else you'll receive verification errors on `Core.Compiler.IRCode`.
""", optimize!)

allow(f::C, args...) where {C <: CompilationContext} = false
function allow(ctx::CompilationContext, mod::Module, fn, args...)
    return allow(ctx, mod) || allow(ctx, fn, args...)
end

@doc(
"""
    allow(f::CompilationContext, args...)::Bool

Determines whether the user-defined [`transform`](@ref) and [`optimize!`](@ref) are allowed to look at a lowered `Core.CodeInfo` or `Core.Compiler.IRCode` instance.

The user is allowed to greenlight modules:

```julia
allow(::MyCtx, m::Module) == m == SomeModule
```

or even specific signatures

```julia
allow(::MyCtx, fn::typeof(rand), args...) = true
```
""", allow)


text:
# simulations.jl
# Convenience functions which run and plot the simulation

function overlap(psia,psib)
    overlaps=zeros(eltype(psib), size(psia,1))
    for i in 1:size(psia,1)
        overlap=dot(psia[:,i], psib)
        overlaps[i]=overlap
    end
    return overlaps
end

" Decompose Hamiltonian into Diagonal/S/PE and Off-diag/J/KE elements"
function decompose_H(H)
    S=eye(N).*H # elementwise to select for just diagonal terms
    J=H-S
    return S,J
end

" Plot spectrum of (H)amiltonian, other useful info."
function plot_H(H)
    # display -> use Julia's prettyprinting, a la the REPL
    display(H)
    println("Eigenvalues: ")
    display(eigvals(H))
    println("Min Eigenvalue: ",eigmin(H))

    myvecs=eigvecs(H)
    myvals=eigvals(H)

    println("Eig Vecs: ")
    display(myvecs)

    # TODO: UPDATE THIS FOR PLOTS.jl
    siteenergies=[H[i,i] for i in 1:size(H,2)]
    myplot=lineplot(siteenergies,name="Site Energies",color=:red,width=80)

    for i in 1:size(myvecs,2)
        vec=myvecs[i,:] + myvals
        println("Vec: ",vec)
        lineplot!(myplot,vec,name="Eigenvectors",color=:green)
    end
    print(myplot)
end

framecounter=0 # variable to keep track of which frame / plot for later movie we are in

function outputpng()
    global framecounter

    Plots.png(@sprintf("%05D.png",framecounter)) # Save plot to PNG file; with XXXXX.png filename
    framecounter=framecounter+1
    #println("Just output frame: $framecounter")
end

function SCFthenUnitary(dampening, SCFcycles, Unitarycycles,eigenergy::Bool=false; PNG::Bool=false)

    S,E,H,psi,dipoles=prepare_model()
    density = psi.^2
    if eigenergy
        Energies = zeros((N,SCFcycles+Unitarycycles))
    end
    # Self consistent field loop; Adiabatic response of lattice + polaron
    # Sets up distorted lattice with polaron, before time-based propagation (should you want it)
    for i in 1:SCFcycles
        @printf("\tSCF loop: %d\n",i)
        S,H,psi,density,dipoles = AdiabaticPropagation(S,dipoles,density,E,dampening)
        if eigenergy
            Energies[(i-1)*N+1:i*N]=eigvals(H)
        else
            plot_model(S,psi,density,dipoles,title="Dampening: $dampening SCF (Adiabatic): Cycle $i / $SCFcycles")
        end
        if PNG outputpng() end
    end

    H=diagm(E,-1)+diagm(S)+diagm(E,1) #build full matrix from diagonal elements; for comparison
    psi=eigvecs(H)[:,2] # 1=gnd state, 2=1st excited state, etc.

    # Setup wavefunction for time-based propagation
    psi=psi+1.5*nondispersive_wavepacket(35,4.0)

    #        psi=nondispersive_wavepacket(20,8.0) # Centered on 20, with Width (speed?) 8.0
    #        psi=psi+nondispersive_wavepacket(40,-20.0) # Fight of the wavepackets!
    #        psi=planewave(8.0) # Plane wave, lambda=8.0 lattice units

    dt=1000 # Time step
    # As energy is in eV; hbar=1; we believe unit is ħ/q = ~0.658 fs

    println("Psi: ",psi)
    #myplot=lineplot(psi,name="Psi",color=:red,width=80,ylim=[-1,1])
    for i in 1:Unitarycycles
        @printf("\n\tUnitary Propagation Loop: %d\n",i)
        #        psi=eigvecs(H)[:,1] # gnd state
        S,H,psi,density,dipoles = UnitaryPropagation(dipoles,density,S,E,psi,dt,dampening,slices=1)
        if eigenergy
            j = i+SCFcycles
            Energies[(j-1)*N+1:j*N]=eigvals(H)
        else
            plot_model(S,psi,density,dipoles,title="Dampening: $dampening TDSE: Cycle $i / $Unitarycycles")
        end
        if PNG outputpng() end

        H=diagm(E,-1)+diagm(S)+diagm(E,1) #build full matrix from diagonal elements; for comparison
#        psi=eigvecs(H)[:,2] # 1=gnd state, 2=1st excited state, etc.

        ## Prototype surface hopping implementation
        overlaps=abs(overlap(eigvecs(H),psi)) # Overlap of current time-dep ψ with the ADIABATIC set of Φ for this H.
        closestAdiabaticState=indmax(overlaps) # index of maximum overlap vector; i.e. this is the adiabatic state with max overlap

        println("Wavefunction overlap; polaron c.f. complete set of (adiabatic) states:")
        for j in overlaps
            @printf("%.2f ",j)
        end
        println("\nSum overlaps: ",sum(overlaps)," Mean: ",sum(overlaps)/N)

        println("Maximum overlap with state $closestAdiabaticState =",overlaps[closestAdiabaticState])

        # Should be minimum switching algorithm / partition function sampling
        # For now we just randomly jump
        if i%50==0
            println(" BORED! JUMPING STATE!")
            psi=eigvecs(H)[:,closestAdiabaticState] # TODO: +1 is a lie; but otherwise it just jumps to the ground state again and again

            if eigenergy
                j = i+SCFcycles
                Energies[(j-1)*N+1:j*N]=eigvals(H)
            else
                plot_model(S,psi,density,dipoles,title="JUMPING JACK FLASH TO: $closestAdiabaticState")
            end
            if PNG outputpng() end
        end
    end
    if eigenergy
        return Energies
    end
end

"""
    UnitarySim(dampening,Unitarycycles; slices=5, dt=1.0, PNG::Bool=false)

Simple Unitary (time evolution) simulation.

slices - number of slices in Trotter decomposition of Hamiltonian for matrix exponentiation
dt - Time step. As energy is in eV; hbar=1; we believe unit is ħ/q = ~0.658 fs
"""
function UnitarySim(dampening,Unitarycycles; slices=5, dt=1.0, PNG::Bool=false)
    S,E,H,psi,dipoles=prepare_model()

#    dipoles=[ -(x-N/2)^2*0.04 for x in 1:N ] # quadratic set of dipoles, to give energy slope across simulation
    dipoles=[ (x-N/2)^3*0.0005 for x in 1:N ]

    # Setup wavefunction for time-based propagation
    psi=nondispersive_wavepacket(10,16.0)

    #        psi=nondispersive_wavepacket(20,8.0) # Centered on 20, with Width (speed?) 8.0
    #        psi=psi+nondispersive_wavepacket(40,-20.0) # Fight of the wavepackets!
    #        psi=planewave(8.0) # Plane wave, lambda=8.0 lattice units
    println("Psi: ",psi)

    for i in 1:Unitarycycles
        @printf("\n\tUnitary Propagation Loop: %d\n",i)
        #        psi=eigvecs(H)[:,1] # gnd state
        S,psi,density,dipoles = UnitaryPropagation(dipoles,E,psi,dt,dampening,slices=slices)
        plot_model(S,psi,density,dipoles,title="Dampening: $dampening TDSE: Cycle $i / $Unitarycycles")
        if PNG outputpng() end
    end
end


text:
using LinearAlgebra

#=
SortedOffsets is a vector of vectors (directional arrows) that span a tensor.
!!! Not indexes ( which are usually one based and positive).
Offsets are zero based and signed.

The zero offset
===============

For the Resynthesizer algorithm, for most uses:

!!! Does not include the zero length vector, e.g. for 2D, (0,0)
IOW, a patch is a donut, with a hole in the middle.

When parameters.withReplacement (which is most use cases)
we are synthesizing the center of the patch,
and don't want to compare the in-progress synthesized value with the corpus,
so we exclude the zero offset.
In other words, the being-synthesized point it not its own neighbor.

Otherwise, when we are just finding the best match,
we might as well compare the centers of the patches, target and corpus.

Order
=====

Ordered by distance from the origin.
The ordering makes patches that are compact around the distinguished point.
A compact patch is crucial.

E.G. for 2D the set of shortest distance offsets is [(-1,0), (1,0), (0,-1), (0,1)]

Spanning
========

Includes enough offsets to span an Array (from any corner to the opposite corner.)

Consider this worst case diagram of a target image and its synth region

CCC
CSS
CSZ

Legend:
C = context
S = synth region
Z = distinguished point of a patch (also in the synth region)

To reach from the distinguished point in a corner, to the opposite corner
requires an offset the size of the dimensions of the target.
Also, requires a negative offset.
Consequently, offsets that will reach from any point in the target
to any other point are all the vectors (of all signedness)
whose x,y are the dimensions of the target.
IOW, the length of the vector of sorted offsets is four times the length of the target.

Use of the entire vector
========================
For the Resynthesizer algorithm:
Only in rare use cases are all the sorted offsets actually used.
That is, the above case is a worst case.

Only on the first pass does the Resynthesizer algorithm need offsets that span the entire target.
On subsequent passes, the algorithm only needs offsets that span a patch,
since the synth region has values, so any patch will be contiguous instead of scattered.

The sorted offsets for the first pass is large and expensive.
We only create it once, and pass it.

We don't bother to create a different, smaller instance of SortedOffsets for subsequent passes.

TODO Faster in most use cases to generate lazily.

=#

#=
Not sure why I don't just define a type alias: (maybe it affects performance?)
SortedOffsets = Vector{CartesianIndex{DimensionCount}} where {DimensionCount}
=#
struct SortedOffsets{DimensionCount}
    offsets::Vector{CartesianIndex{DimensionCount}}
end

#=
Outer constructor from a tensor.

Create from SortedOffets from the tensor the offsets should span.

Minimizing allocations is not too important because this is only called once.
=#
function SortedOffsets(
    tensor::Array{ValueType, DimensionCount}
    ) where {ValueType, DimensionCount}
    # mock
    # return createMockSortedOffsets()

    offsets = createOffsets(tensor)

    #=
    sort() alone would sort on y, then x
    Instead sort by manhattan AKA cityblock distance

    sort! sorts in place (but with no fewer allocations?)
    =#
    sortedOffsets = sort!(offsets, by=cityBlockDistance)

    #=
    if parameters.withReplacement
        sortedOffsets = excludeFirstElement(sortedOffsets)
    end
    =#

    # call default constructor to encapsulate vector in a struct
    result = SortedOffsets(sortedOffsets)
    @debug "SortedOffsets of type, size" typeof(result) size(result.offsets)
    @debug "SortedOffsets" result.offsets
    return result
end

#=
I struggled to get this performant.
It still does many allocations that could probably be eliminated,
but it is only done once.

Strategies:

pop() or deleteat() does 25M allocations for medium data
deleteat!(sortedOffsets, 1)

Slice is performant.  It creates a view, or a copy?

This still does many allocations??? TODO
sortedOffsets[2:end]

FAIL, no signature for SubArray
@view sortedOffsets[2:end]

# Assert size is one less than before
=#

# Exclude the first element, the zero offset
function excludeFirstElement(vector)
    result = vector[2:length(vector)]
end

#=
Create a vector of all the arrows to span a tensor.

The strategy is to generate a large array (four times larger than tensor)
of CartesianIndex of suitable ranges (negative:positve), then convert
to a Vector.

Another strategy would be to generate all words of (-1,1) of length k,
then multiply by CartesianIndices of the tensor, shifted by substracting 1,
to get zeros.
=#
function createOffsets(
                tensor::Array{ValueType, DimensionCount}
            )::Vector{CartesianIndex{DimensionCount}}  where{ValueType, DimensionCount}

    sizeVec = size(tensor)
    # sizeVec is a tuple, one element per dimension, giving size on an axis

    rawArray = createArrayofSpanningIndex(sizeVec)

    rawVector = vec(rawArray)   # reshape from array to vector
    # !!! Contains (0,0) still
    return rawVector
end

#=
Create array of CartesianIndex{DimensionCount}
where elements are indices in range -size:size

Since in all generality, there are infinite sizeVec's (for DimensionCount in range 1:Inf)
we could use @generated function,
Or use runtime determined
=#
function createArrayofSpanningIndex(sizeVec)
    #=
    Runtime dispatch.
    Performance is not an issue, we only do this once.
    But it is not fully general for all inputs to algorithm beyond dimension coded here.
    =#
    # Using comprehension
    countDimensions = length(sizeVec)
    if countDimensions == 1
        # TODO 1D indexes are Int64, not CartesianIndex{1}
        return [ CartesianIndex(x) for x = -sizeVec[1]:sizeVec[1] ]
    elseif countDimensions == 2
        return [ CartesianIndex(x,y) for x = -sizeVec[1]:sizeVec[1], y = -sizeVec[2]:sizeVec[2]  ]
    elseif countDimensions == 3
        return [ CartesianIndex(x,y,z) for  x = -sizeVec[1]:sizeVec[1],
                                            y = -sizeVec[2]:sizeVec[2],
                                            z = -sizeVec[3]:sizeVec[3]  ]
    elseif countDimensions == 4
        return [ CartesianIndex(x,y,z,u) for x = -sizeVec[1]:sizeVec[1],
                                            y = -sizeVec[2]:sizeVec[2],
                                            z = -sizeVec[3]:sizeVec[3],
                                            u = -sizeVec[4]:sizeVec[4]  ]
    else
        throw(ErrorException("Unhandled dimension > 4"))
    end
end

#=
Manhattan/taxicab/L1 distance of a vector is norm(vector,1) using LinearAlgebra pkg.
The so-called 1-norm.  Here the literal 1 means the 1-norm.
That is, the literal 1 is the value of the argument with keyword "p",
i.e. this is from a family of norms called p-norms.

Result is a scalar.

Alternatively use Distances pkg?
=#
function cityBlockDistance(coordinates) # ::CartesianIndex{DimensionCount})
    # convert to Tuple because norm requires an iterable and CartsianIndex is not iterable?
    return norm(Tuple(coordinates), 1)
end

#=
Create a short SortedOffsets, for testing.
=#
function createMockSortedOffsets()
    return [
        # !!! Exclude CartesianIndex(0,0),
        CartesianIndex(0,1),
        CartesianIndex(0,-1)
    ]
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module TableWidgets

using Formatting
using IJulia.CommManager

export TextWidget, writetable, update, TableWidget

include("textwidget.jl")
include("writetable.jl")
include("tablewidget.jl")

end # module


text:
using LinearAlgebra: qr, diagm, diag
using Flux: glorot_uniform, gate, zeros, @functor
import Flux: RNNCell, LSTMCell, GRUCell, hidden, Recur

export orthogonal

RNNCell(in::Integer, out::Integer, σ = tanh; init = glorot_uniform, kernel_init = glorot_uniform) =
  RNNCell(σ, init(out, in), glorot_uniform(out, out), init(out), zeros(out))

function LSTMCell(in::Integer, out::Integer; init = glorot_uniform, kernel_init = glorot_uniform)
    cell = LSTMCell(init(out * 4, in), kernel_init(out * 4, out), init(out * 4), zeros(out), zeros(out))
    cell.b[gate(out, 2)] .= 1
    return cell
end

GRUCell(in, out; init = glorot_uniform, kernel_init = glorot_uniform) =
    GRUCell(init(out * 3, in), kernel_init(out * 3, out), init(out * 3), zeros(out))

"""
    orthogonal(dim)
Return a random orthogonal maxtrix of size `(dim, dim)`. This is commonly used for kernel initialization of recurrent layers.
# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> Flux.orthogonal(2)
2×2 Array{Float32,2}:
 -0.633973  -0.773356
 -0.773356   0.633973
```
See [Exact solutions to the nonlinear dynamics of learning in deep linear neural networks](http://arxiv.org/abs/1312.6120)
"""
orthogonal(dim) = orthogonal_matrix(dim, dim)

function orthogonal(nrow, ncol)
    @assert nrow >= ncol && nrow % ncol == 0
    vcat([orthogonal(ncol) for n in 1:(nrow ÷ ncol)]...)
end

"""
    orthogonal_matrix(nrow, ncol)
If the shape of the matrix to initialize is two-dimensional, it is initialized
with an orthogonal matrix obtained from the QR decomposition of a matrix of
random numbers drawn from a normal distribution.
If the matrix has fewer rows than columns then the output will have orthogonal
rows. Otherwise, the output will have orthogonal columns.
"""
function orthogonal_matrix(nrow, ncol)
    shape = reverse(minmax(nrow, ncol))
    a = randn(Float32, shape)
    q, r = qr(a)
    q = Matrix(q) * diagm(sign.(diag(r)))
    nrow < ncol ? permutedims(q) : q
end

mutable struct GGRUCell{G, A, V}
    g::G
    Wi::A
    Wh::A
    b::V
    h::V
end

GGRUCell(in, out, g = 1f0; init = glorot_uniform, kernel_init = glorot_uniform) =
    GGRUCell(g, init(out * 3, in), kernel_init(out * 3, out), init(out * 3), zeros(out))

function (m::GGRUCell)(h, x)
    b, o = m.b, size(h, 1)
    gx, gh = m.Wi * x, m.Wh * h
    r = σ.((r_ = gate(gx, o, 1) .+ gate(gh, o, 1) .+ gate(b, o, 1);))
    z = σ.((z_ = gate(gx, o, 2) .+ gate(gh, o, 2) .+ gate(b, o, 2);))
    h̃ = tanh.(m.g * (h̃_ = gate(gx, o, 3) .+ r .* gate(gh, o, 3) .+ gate(b, o, 3);))
    h′ = (1 .- z) .* h̃ .+ z .* h
    return h′, h′
end

hidden(m::GGRUCell) = m.h

@functor GGRUCell

Base.show(io::IO, l::GGRUCell) =
    print(io, "GGRUCell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷3, ")")

GGRU(a...; ka...) = Recur(GGRUCell(a...; ka...))


text:
__precompile__(true)

module TenhouLogs

include("playdata.jl")
include("parseutils.jl")
include("parserdict.jl")

using SQLite
using Tables
using Dates
using Downloads
using ZipFile
using CodecZlib
using CodecLz4

include("dbcreate.jl")
include("analyse.jl")

end


text:
#
# GLVisualize
# - Write a C program to convert FBX with animations
# -

max_stars_repo_name: laxtools/article

max_issues_repo_name: laxtools/article

max_forks_repo_name: laxtools/article


text:
using RenyiEntropyOfMagic
RenyiEntropyOfMagic.command_main()


text:
using DashBootstrapComponents, DashHtmlComponents

cards = html_div([
    dbc_card(dbc_cardbody("This is some text within a card body"), className = "mb-3"),
    dbc_card("This is also within a body", body = true),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using StaticTools
using Test
using ManualMemory: MemoryBuffer

const GROUP = get(ENV, "GROUP", "All")

@static if GROUP == "Core" || GROUP == "All"
    @testset "IO" begin include("testllvmio.jl") end
    @testset "libc" begin include("testllvmlibc.jl") end
    @testset "Parsing" begin include("testparse.jl") end
    @testset "StaticString" begin include("teststaticstring.jl") end
    @testset "MallocString" begin include("testmallocstring.jl") end
    @testset "MallocArray" begin include("testmallocarray.jl") end
    @testset "StaticRNG" begin include("teststaticrng.jl") end
    @testset "MemoryBuffer" begin

        # Test direct buffer constructor
        buf = mm"Hello, world! 🌍"
        @test isa(buf, MemoryBuffer{18, UInt8})

        # Test ascii escaping
        a = mm"\0\a\b\f\n\r\t\v'\"\\\0"::MemoryBuffer
        b = MemoryBuffer((codeunits("\0\a\b\f\n\r\t\v'\"\\\0")...,))
        @test a.data == b.data

    end
end

using LoopVectorization
@static if GROUP == "Integration" || GROUP == "All"
    @testset "StaticCompiler" begin include("teststaticcompiler.jl") end
end


text:
setup = CylinderSetup(
    name = "cylinders/somecylinders",
    ncell = 2,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    bend = 0.0,
    twist = 0.0,
    height = 20.0,
    include_in = true,
    in_ratio = 0.6,
    ecs_shape = "convex_hull",
    ecs_ratio = 0.5,
    refinement = 1,
    D_in = 0.002 * I(3),
    D_out = 0.002 * I(3),
    D_ecs = 0.002 * I(3),
    T₂_in = Inf,
    T₂_out = Inf,
    T₂_ecs = Inf,
    ρ_in = 1.0,
    ρ_out = 1.0,
    ρ_ecs = 1.0,
    κ_in_out = 1e-3,
    κ_out_ecs = 1e-3,
    κ_in = 0.0,
    κ_out = 0.0,
    κ_ecs = 0.0,
)

experiment = Experiment(
    gradient = (
        directions = create_directions([0.0; 0.0; 1.0]),
        sequences = [PGSE(1000.0, 5000.0)],
        values = [2000.0],
        values_type = "b",
    ),
    btpde = (odesolver = QNDF(), reltol = 1e-4, abstol = 1e-6, nsave = 1),
    btpde_midpoint = (θ = 0.5, timestep = 5),
    # mf = (length_scale = 3, neig_max = 400, ninterval = 500),
    hadc = (odesolver = QNDF(), reltol = 1e-4, abstol = 1e-6),
    karger = (odesolver = MagnusGL6(), timestep = 5),
)


text:
function optimize(f::Function, lower::T, upper::T;
        alg::UnivariateAlgorithm=GoldenSection(),
        reltol=sqrt(eps(T)),
        abstol=eps(T),
        maxiter::Integer=1_000) where {T<:AbstractFloat}
    lower < upper || error("Variable lower bound has to be smaller than upper bound.")
    _optimize(f, lower, upper, alg; reltol=T(reltol), abstol=T(abstol), maxiter)
end

"""
    optimize(f::Function, lower::Real, upper::Real; kwargs...)

Compute the minimizer of the function `f` using the bracket (`lower`, `upper`).

# Keywords
- `alg::UnivariateAlgorithm=GoldenSection()`: the algorithm choice
- `reltol=sqrt(eps(T))`: relative tolerance
- `abstol=eps(T)`: absolute tolerance
- `maxiter::Integer=1_000`: maximum numer of iterations
"""
function optimize(f::Function, lower::Real, upper::Real; kwargs...)
    T = promote_type(typeof(lower/1), typeof(upper/1))
    optimize(f, T(lower), T(upper); kwargs...)
end


text:
### FlatS2 types
# (spin-2 fields are just special FieldTuples)
const FlatQUMap{P,T,M}     = FieldTuple{QUMap,     NamedTuple{(:Q,:U),NTuple{2,FlatMap{P,T,M}}}, T}
const FlatQUFourier{P,T,M} = FieldTuple{QUFourier, NamedTuple{(:Q,:U),NTuple{2,FlatFourier{P,T,M}}}, Complex{T}}
const FlatEBMap{P,T,M}     = FieldTuple{EBMap,     NamedTuple{(:E,:B),NTuple{2,FlatMap{P,T,M}}}, T}
const FlatEBFourier{P,T,M} = FieldTuple{EBFourier, NamedTuple{(:E,:B),NTuple{2,FlatFourier{P,T,M}}}, Complex{T}}
# some handy Unions
const FlatS2{P,T,M}=Union{FlatEBMap{P,T,M},FlatEBFourier{P,T,M},FlatQUMap{P,T,M},FlatQUFourier{P,T,M}}
const FlatQU{P,T,M}=Union{FlatQUMap{P,T,M},FlatQUFourier{P,T,M}}
const FlatEB{P,T,M}=Union{FlatEBMap{P,T,M},FlatEBFourier{P,T,M}}
const FlatS2Map{P,T,M}=Union{FlatQUMap{P,T,M},FlatEBMap{P,T,M}}
const FlatS2Fourier{P,T,M}=Union{FlatQUFourier{P,T,M},FlatEBFourier{P,T,M}}

### convenience constructors
for (F,F0,(X,Y),T) in [
        (:FlatQUMap,     :FlatMap,     (:Qx,:Ux), :T),
        (:FlatQUFourier, :FlatFourier, (:Ql,:Ul), :(Complex{T})),
        (:FlatEBMap,     :FlatMap,     (:Ex,:Bx), :T),
        (:FlatEBFourier, :FlatFourier, (:El,:Bl), :(Complex{T}))
    ]
    doc = """
        # main constructor:
        $F($X::AbstractMatrix, $Y::AbstractMatrix[, θpix={resolution in arcmin}, ∂mode={fourier∂ or map∂})

        # more low-level:
        $F{P}($X::AbstractMatrix, $Y::AbstractMatrix) # specify pixelization P explicilty
        $F{P,T}($X::AbstractMatrix, $Y::AbstractMatrix) # additionally, convert elements to type $T
        $F{P,T,M<:AbstractMatrix{$T}}($X::M, $Y::M) # specify everything explicilty

    Construct a `$F` object. The top form of the constructor is most convenient
    for interactive work, while the others may be more useful for low-level code.
    """
    @eval begin
        @doc $doc $F
        $F($X::AbstractMatrix, $Y::AbstractMatrix; kwargs...) = $F{Flat(Nside=size($X,2);kwargs...)}($X, $Y)
        $F{P}($X::AbstractMatrix, $Y::AbstractMatrix) where {P} = $F{P}($F0{P}($X), $F0{P}($Y))
        $F{P,T}($X::AbstractMatrix, $Y::AbstractMatrix) where {P,T} = $F{P,T}($F0{P,T}($X), $F0{P,T}($Y))
    end
end

### properties
function propertynames(f::FlatS2)
    (:fs, propertynames(f.fs)...,
     (Symbol(string(k,(f isa FlatMap ? "x" : "l"))) for (k,f) in pairs(f.fs))...)
end
getproperty(f::FlatQUMap,     ::Val{:Qx}) = getfield(f,:fs).Q.Ix
getproperty(f::FlatQUMap,     ::Val{:Ux}) = getfield(f,:fs).U.Ix
getproperty(f::FlatQUFourier, ::Val{:Ql}) = getfield(f,:fs).Q.Il
getproperty(f::FlatQUFourier, ::Val{:Ul}) = getfield(f,:fs).U.Il
getproperty(f::FlatEBMap,     ::Val{:Ex}) = getfield(f,:fs).E.Ix
getproperty(f::FlatEBMap,     ::Val{:Bx}) = getfield(f,:fs).B.Ix
getproperty(f::FlatEBFourier, ::Val{:El}) = getfield(f,:fs).E.Il
getproperty(f::FlatEBFourier, ::Val{:Bl}) = getfield(f,:fs).B.Il
getproperty(f::FlatS2,        ::Val{:P})  = f
function getindex(f::FlatS2, k::Symbol)
    B = @match k begin
        (:P)         => Basis
        (:E  || :B)  => f isa FlatQUMap ? EBMap : f isa FlatQUFourier ? EBFourier : Basis
        (:Q  || :U)  => f isa FlatEBMap ? QUMap : f isa FlatEBFourier ? QUFourier : Basis
        (:Ex || :Bx) => EBMap
        (:El || :Bl) => EBFourier
        (:Qx || :Ux) => QUMap
        (:Ql || :Ul) => QUFourier
        _ => throw(ArgumentError("Invalid FlatS2 index: $k"))
    end
    getproperty(B(f),k)
end

### basis conversion

QUFourier(f::FlatQUMap) = FlatQUFourier(map(Fourier,f.fs))
QUFourier(f::FlatEBMap) = f |> EBFourier |> QUFourier
QUFourier(f::FlatEBFourier{P,T}) where {P,T} = begin
    @unpack sin2ϕ, cos2ϕ = fieldinfo(f)
    Ql = @. - f.El * cos2ϕ + f.Bl * sin2ϕ
    Ul = @. - f.El * sin2ϕ - f.Bl * cos2ϕ
    FlatQUFourier(Q=FlatFourier{P}(Ql), U=FlatFourier{P}(Ul))
end

QUMap(f::FlatQUFourier)  = FlatQUMap(map(Map,f.fs))
QUMap(f::FlatEBMap)      = f |> EBFourier |> QUFourier |> QUMap
QUMap(f::FlatEBFourier)  = f |> QUFourier |> QUMap

EBFourier(f::FlatEBMap) = FlatEBFourier(map(Fourier,f.fs))
EBFourier(f::FlatQUMap) = f |> QUFourier |> EBFourier
EBFourier(f::FlatQUFourier{P,T}) where {P,T} = begin
    @unpack sin2ϕ, cos2ϕ = fieldinfo(f)
    El = @. - f.Ql * cos2ϕ - f.Ul * sin2ϕ
    Bl = @.   f.Ql * sin2ϕ - f.Ul * cos2ϕ
    FlatEBFourier(E=FlatFourier{P}(El), B=FlatFourier{P}(Bl))
end

EBMap(f::FlatEBFourier) = FlatEBMap(map(Map,f.fs))
EBMap(f::FlatQUMap)     = f |> QUFourier |> EBFourier |> EBMap
EBMap(f::FlatQUFourier) = f |> EBFourier |> EBMap

QUFourier(f′::FlatQUFourier, f::FlatQUMap) = (tmap(Fourier,f′.fs,f.fs); f′)
QUMap(f′::FlatQUMap, f::FlatQUFourier) = (tmap(Map,f′.fs,f.fs); f′)

Map(f::FlatQUFourier) = QUMap(f)
Map(f::FlatEBFourier) = EBMap(f)
Map(f::FlatS2Map) = f
Fourier(f::FlatQUMap) = QUFourier(f)
Fourier(f::FlatEBMap) = EBFourier(f)
Fourier(f::FlatS2Fourier) = f

### simulation and power spectra
function white_noise(::Type{F2}) where {F2<:FlatS2}
    F = (((::Type{<:FieldTuple{B,NamedTuple{Names,NTuple{2,F}}}}) where {B,Names,F}) -> F)(F2)
    FlatEBMap(E=white_noise(F), B=white_noise(F))
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S2}, CℓEE::InterpolatedCℓs, CℓBB::InterpolatedCℓs) where {P,T,M}
    Diagonal(FlatEBFourier(E=Cℓ_to_Cov(P,T,S0,CℓEE).diag, B=Cℓ_to_Cov(P,T,S0,CℓBB).diag))
end

function get_Cℓ(f1::FlatS2, f2::FlatS2=f1; which=(:EE,:BB), kwargs...)
    Cℓ = (;[Symbol(x1*x2) => get_Cℓ(getindex(f1,Symbol(x1)),getindex(f2,Symbol(x2)); kwargs...) for (x1,x2) in split.(string.(ensure1d(which)),"")]...)
    which isa Symbol ? Cℓ[1] : Cℓ
end

function ud_grade(f::FlatS2{P}, args...; kwargs...) where {P}
    f′ = FieldTuple(map(f->ud_grade(f, args...; kwargs...), f.fs))
    B′ = (f′[1] isa FlatMap) ? (f isa FlatQU ? QUMap : EBMap) : (f isa FlatQU ? QUFourier : EBFourier)
    FieldTuple{B′}(f′)
end


text:
using SpecialFunctions

"""
    mc_volume(i, xs, nmc, searcher)

Estimate the area and volume of the `i`-th Voronoi cell from the Voronoi Diagram generated by `xs`
by a Monte Carlo estimate from random rays
"""
function mc_volume(i, xs, nmc=1000, searcher = Raycast(xs))
    y, δy, V, A = mc_integrate(x->1, i, xs, nmc, 0, searcher)
    A, V
end

"""
    mc_volumes(xs::Points, nmc=1000)

Estimate the areas and volumes of the Voronoi Cells generated by `xs` using `nmc` Monte Carlo samples.

# Returns
- `SparseMatrix`: the areas of the common boundaries of two cells
- `Vector`: the volumes of each cell
"""
function mc_volumes(xs::Points, nmc=1000)
    V = zeros(length(xs))
    A = spzeros(length(xs), length(xs))

    searcher = Raycast(xs)

    for i in 1:length(xs)
        a, v = mc_volume(i, xs, nmc, searcher)
        V[i] = v
        A[:, i] = a
    end
    #A = (A + A') / 2
    A, V
end

"""
    mc_volumes(sig::Vertices, xs::Points, nmc=1000)

In the case when the simplicial complex is already known this information can be used
to speed up the Monte-Carlo sampling by restricting the search space
"""
function mc_volumes(σ::Vertices, xs::Points, nmc=1000)
    V = zeros(length(xs))
    A = spzeros(length(xs), length(xs))

    searcher = Raycast(xs)
    neigh = neighbors(σ)

    for i in 1:length(xs)
        ix = [i; neigh[i]]
        a, v = mc_volume(1, xs[ix], nmc, RaycastBruteforce())
        V[i] = v
        A[ix, i] = a
    end
    #A = (A + A') / 2
    A, V
end

"""
    neighbors(sig::Vertices) --> Dict{Int, Vector{Int}}
Compute the neighbors of all generators.
"""
neighbors(σ::Vertices) = neighbors(keys(σ))
function neighbors(σ)
    d = Dict{Int, Vector{Int}}()
    for s in σ
        for i in s
            v = get!(d, i, Vector{Int}())
            append!(v, s)
        end
    end
    for i in keys(d)
        unique!(d[i])
        filter!(x->x != i, d[i])
    end
    return d
end

"""
    mc_integrate(f::Function, i::Int, xs::Points, nmc=1000, nmc2=1000, searcher=Raycast(xs))
Integrate function `f` over cell `i` and its boundary using `nmc` rays per cell and `nmc2` points per ray for the volume integral.

# Returns:
- `Vf::Real`: the volume integral of `f`
- `Af::SparseVector`: Af[j] is the surface integral of `f` over the intersection between cells i and j.
- `V::Real`: the volume of cell `i`
- `A::SparseVector`: A[j] is the surface area of the intersection between cells i and j.
"""
function mc_integrate(f::Function, i::Int, xs::Vector, nmc=1000, nmc2=1000, searcher = Raycast(xs))
    x = xs[i]
    d = length(x)

    # volume/area computations are basically for free, leaving them commented for now
    V = 0.
    A = spzeros(length(xs))

    y = 0.
    δy = spzeros(length(xs))

    for _ in 1:nmc
        u = normalize(randn(d))
        (j, t) = raycast([i], x, u, xs, searcher)

        V += t^d
        for _ in 1:nmc2
            r = t * rand()
            x′ = x + u * r
            y += f(x′) * r^(d-1) * t
        end

        if t < Inf
            j = j[1] == i ? j[2] : j[1]
            normal = normalize(xs[j] - xs[i])
            dA = t ^ (d-1) / abs(dot(normal, u))
            A[j] += dA
            δy[j] += dA * f(x + t*u)
        end
    end

    c_vol = pi^(d/2) / gamma(d/2 + 1)
    c_area = d * c_vol

    V *= c_vol / nmc
    A *= c_area / nmc

    y *= c_area / nmc / nmc2
    δy *= c_area / nmc

    y, δy, V, A
end

function mc_integrate(f, σ::Vertices, xs::Points, nmc, nmc2)
    return mc_integrate_i((x,i)->f(x), σ, xs, nmc, nmc2)
end

# MC integration using f(x, i) function which is also passed the index of the cell
function mc_integrate_i(f, σ::Vertices, xs::Points, nmc, nmc2)
    n = length(xs)
    ys = zeros(n)
    ∂ys = zeros(n)

    neigh = neighbors(σ)

    for i in 1:n
        ix = [i; neigh[i]]
        y, ∂y, V, A = mc_integrate(x->f(x, i), 1, xs[ix], nmc, nmc2, RaycastBruteforce())
        ys[i] = y
        ∂ys[i] = sum(∂y)
    end
    return ys, ∂ys
end


text:
using Test, SafeTestsets

@safetestset "Political Democracy" begin include("political_democracy/political_democracy.jl") end
@safetestset "Recover Parameters" begin include("recover_parameters/recover_parameters_twofact.jl") end
@safetestset "Multigroup" begin include("multigroup/multigroup.jl") end


text:
using Printf
using Random
using LinearAlgebra
using Distributed
using MAT
@everywhere using DistributedArrays
@everywhere using RCAM
@everywhere using Random
@everywhere Random.seed!(123)

#Load the data
fid = matopen("../data/Data.mat")
d = read(fid)
X = d["D"]

function pCholtest(X)
    # Choltest is just a wrapper for this script

    # Params
    r = 50 #Rank
    alt = 5 #Iterations
    perc = 0.8 #Percent of missing values
    s1,s2 = size(X)

    # Subsample
    ind = randperm(s1*s2)
    inds = ind[1:Int(round(perc*s1*s2))]
    b = copy(X)
    b[inds] .= 0

    # Interpolate
    @time ownerL, ownerR = dclr(b, r, alt)

    # Gather L and R results
    L = fetch(@spawnat ownerL Main.L)
    R = fetch(@spawnat ownerR Main.R)

    return L,R
end

snr(raw,interp) = -20*log10(norm(interp-raw)/norm(raw))

L,R = pCholtest(X)
println(snr(X,L*R'))


text:
mutable struct State{A<:AbstractVector{<:Complex}}
  ψ::A
end

State(ψ::AbstractVector{<:Real}) = State(float(complex(ψ)))

function State(n::Integer, vs::NTuple{N,Integer}) where N
  ψ = zeros(ComplexF64, 2^n)
  foreach(v -> ψ[v+1] = 1/√length(vs), vs)
  return State(ψ)
end

State(n::Integer) = State(n, (0,))

nbits(s::State) = nbits(s.ψ)
Base.length(s::State) = nbits(s)

labels(s::State) = map(n -> Base.bitstring(n-1)[end-nbits(s)+1:end], 1:length(s.ψ))

probabilities(s::State) = abs2.(s.ψ)

function Base.show(io::IO, s::State; n = 5)
  basis = reverse(collect(Iterators.filter(x -> x[2] ≉ 0, zip(labels(s), s.ψ))))
  length(basis) == 1 && basis[1][2] ≈ 1 && return print(io, "|$(basis[1][1])⟩")
  join(io, ["|$(l)⟩ * ($p)" for (l, p) in take(basis, n)], " +\n")
  length(basis) > n && print(io, "\n    ⋮")
end

function apply!(U::AbstractMatrix, s::State)
  s.ψ = U * s.ψ
  return s
end

function apply!(U::AbstractMatrix, s::State, is...)
  @assert nbits(U) == length(is)
  U = pad(U, 0, nbits(s)-nbits(U))
  S = headfirst(nbits(s), is...)
  apply!(S'*U*S, s)
end

function probabilities(s::State, is...)
  n = length(is)
  M = measure(n, nbits(s)-length(is))*headfirst(nbits(s), is...)
  M*probabilities(s)
end

function sample(ps)
  x = rand()
  s = zero(eltype(ps))
  for i = 1:length(ps)
    s += ps[i]
    s > x && return i
  end
end

function measure!(s::State, is...)
  n = length(is)
  M = measure(n, nbits(s)-length(is))*headfirst(nbits(s), is...)
  ps = M*probabilities(s)
  i = sample(ps)
  s.ψ .*= M[i,:]./sqrt(ps[i])
  bits(i-1, Val{length(is)})
end


text:
##############################################################################
##
## TODO
##
##############################################################################

# - Convert to c function. Hard to read so many ccall's.

function init()
  ccall((:OpenSSL_add_all_digests, "libcrypto"), Void, ())
end

function cleanup()
  ccall((:EVP_cleanup, "libcrypto"), Void, ())
end

function digest(name::AbstractString, data::ByteString; is_hex=false)
  if is_hex
    L = round(Int,length(data)/2)
    data = [@compat(parse(@compat(UInt8),SubString(data,2*i-1,2*i), 16)) for i in 1:L]
  else
    data = data.data
  end
  digest(name, data)
end

function digest(name::AbstractString, data::Array{@compat(UInt8)})
  ctx = ccall((:EVP_MD_CTX_create, "libcrypto"), Ptr{Void}, ())
  try
    # Get the message digest struct
    md = ccall((:EVP_get_digestbyname, "libcrypto"), Ptr{Void}, (Ptr{Uint8},), bytestring(name))
    if(md == C_NULL)
      error("Unknown message digest $name")
    end

    # Add the digest struct to the context
    ccall((:EVP_DigestInit_ex, "libcrypto"), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, md, C_NULL)

    # Update the context with the input data
    ccall((:EVP_DigestUpdate, "libcrypto"), Void, (Ptr{Void}, Ptr{Uint8}, Uint), ctx, data, length(data))

    # Figure out the size of the output string for the digest
    size = ccall((:EVP_MD_size, "libcrypto"), Uint, (Ptr{Void},), md)
    uval = Array(Uint8, size)

    # Calculate the digest and store it in the uval array
    ccall((:EVP_DigestFinal_ex, "libcrypto"), Void, (Ptr{Void}, Ptr{Uint8}, Ptr{Uint}), ctx, uval, C_NULL)

    return uval
  finally
    ccall((:EVP_MD_CTX_destroy, "libcrypto"), Void, (Ptr{Void},), ctx)
  end
end

digest(name::AbstractString, data::IOBuffer) = digest(name, takebuf_array(data))


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for algorithms.
############################################################################################################################

# Traversals
@testset "Algorithms" begin
   g = loadgraph("testgraph.txt")

   @test bfs(g, 1) == [0,1,1,3,4,4,4,4,4,4]
   @test bfs(g, 1, 2) == [0,1,1,3,-1,-1,-1,-1,-1,-1]
   @test bfs(g, 1, 3) == [0,1,1,3,4,4,4,4,4,4]

   @test bfs_list(g, 1) == [2,3,4,5,6,7,8,9,10]
   @test bfs_list(g, 2) == [1,3,4,5,6,7,8,9,10]
   @test bfs_list(g, 1, 1, 2) == [2,3,4]
   @test bfs_list(g, 1, 2, 3) == [4,5,6,7,8,9,10]

   g1 = bfs_tree(g, 1, 2)
   @test nv(g1) == 4
   @test ne(g1) == 3

   g2 = bfs_tree(g, 1, 4)
   @test nv(g2) == 10
   @test ne(g2) == 9

   g3 = bfs_subgraph(g, 1, 2)
   @test nv(g3) == 4
   @test ne(g3) == 8

   g4 = bfs_subgraph(g, 1, 4)
   @test nv(g4) == 10
   @test ne(g4) == 28

   ###
   # External interface
   ###
   setlabel!(g, map(string, 1:10))
   @test hoplist(g, "1", 1, 2) == ["2","3","4"]

   g5 = hoptree(g, "1", 2)
   @test nv(g5) == 4
   @test ne(g5) == 3

   g6 = hopgraph(g, "1", 2)
   @test nv(g6) == 4
   @test ne(g6) == 8

   ###
   # Mutual friends
   ###
   g = completegraph(10)
   setlabel!(g, map(string, 1:10))
   @test mutualcount(g, "1", "2") == 8
   @test mutual(g, "1", "2") == map(string, [3,4,5,6,7,8,9,10])
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

#Pkg.add("ForwardDiff")
#Pkg.update()
#Pkg.checkout("ForwardDiff")
#Pkg.status()
using ForwardDiff

include("common.jl")

# matrix version which works with ForwardDiff
function logsumexp_AD_mat(X)
  n = size(X,2)
  mX = [maximum(X[:,i]) for i in 1:n]
  log(sum(exp(X .- mX'),1)) + mX'
end
# input should be 1 dimensional
function logsumexp_AD_vec(x)
  mx = maximum(x)
  log(sum(exp(x - mx))) + mx
end

@inbounds function gmm_objective_AD(alphas,means,icf,x,wishart::Wishart)
  d = size(x,1)
  n = size(x,2)
  CONSTANT = -n*d*0.5*log(2 * pi)

  sum_qs = sum(icf[1:d,:],1)
  slse = sum(sum_qs)
  Qs = [get_Q(d,icf[:,ik]) for ik in 1:k]
  main_term = zeros(eltype(alphas),k,n)

  for ik=1:k
    Qxcentered = Qs[ik] * (x .- means[:,ik])
    main_term[ik,:] = alphas[ik] + sum_qs[ik] - 0.5*sumabs2(Qxcentered,1)
  end

  CONSTANT + sum(logsumexp_AD_mat(main_term)) - n*logsumexp_AD_vec(alphas) + log_wishart_prior(wishart, sum_qs, Qs, icf)
end

function logsumexp(X,axis=2)
  mX = maximum(X,axis)
  log(sum(exp(X .- mX),axis)) + mX
end
typealias Mat Matrix{Float64}

@inbounds function gmm_objective(alphas::Mat,means::Mat,icf::Mat,x::Mat,wishart::Wishart)
  d = size(x,1)
  n = size(x,2)
  CONSTANT = -n*d*0.5*log(2 * pi)

  sum_qs = sum(icf[1:d,:],1)
  slse = sum(sum_qs)
  Qs = [get_Q(d,icf[:,ik]) for ik in 1:k]
  main_term = zeros(eltype(alphas),k,n)

  for ik=1:k
    Qxcentered = Qs[ik] * (x .- means[:,ik])
    main_term[ik,:] = alphas[ik] + sum_qs[ik] - 0.5*sumabs2(Qxcentered,1)
  end

  lse_alphas = n*logsumexp(alphas);
  CONSTANT + sum(logsumexp(main_term,1)) - lse_alphas[1] + log_wishart_prior(wishart, sum_qs, Qs, icf)
end

# Read instance
#dir_in = "C:\\Users\\t-filsra\\Workspace\\autodiff\\"
#dir_out = "C:\\Users\\t-filsra\\Workspace\\autodiff\\"
#fn = "gmm"
dir_in = ARGS[1]
dir_out = ARGS[2]
fn = ARGS[3]
nruns_f = parse(Int,ARGS[4])
nruns_J = parse(Int,ARGS[5])
replicate_point = size(ARGS,1) >= 6 && ARGS[6] == "-rep"

fn_in = string(dir_in, fn)
fn_out = string(dir_out, fn)

alphas,means,icf,x,wishart = read_gmm_instance(string(fn_in,".txt"),replicate_point)
d = size(means,1)
k = size(means,2)
n = size(x,2)

# Objective
precompile(gmm_objective,(typeof(alphas),typeof(means),typeof(icf),typeof(x),Wishart))
err = 0.
tic()
for i in 1:nruns_f
  err = gmm_objective(alphas,means,icf,x,wishart)
end
tf = toq()/nruns_f
@printf "tf: %g\n" tf
#println(err)

# Gradient helper
function wrapper_gmm_objective_AD(packed)
  alphas,means,icf = unpack(d,k,packed)
  gmm_objective_AD(alphas,means,icf,x,wishart)
end

# Gradient
g = ForwardDiff.gradient(wrapper_gmm_objective_AD)
precompile(pack,(typeof(alphas),typeof(means),typeof(icf)))
precompile(g,(typeof(alphas),typeof(means),typeof(icf)))
J = zeros(1,length(alphas)+length(means)+length(icf))
tic()
for i in 1:nruns_J
  J = g(pack(alphas,means,icf))
end
tJ = toq()/nruns_J
@printf "tJ: %g\n" tJ
#println("J:")
#println(J)

name = "Julia_F_vector"

write_J(string(fn_out,"_J_",name,".txt"),J)
write_times(string(fn_out,"_times_",name,".txt"),tf,tJ)

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench