text:
#test time evolution (DifferentialEquations 4.0)
using ProjectedGPE

function timeevolution2()
sinfo = Params()

#== start template ==#
#fundamental constants/units
  ħ = 1.0545718e-34
  kB = 1.38064852e-23
  amu = 1.660339040e-27
  Bohr = 5.29e-11
#Rb87 mass and scattering length
  m = 86.909180527*amu
  as = 100*Bohr
#trap frequencies
  ωx = 2π
  ωy = 4ωx
  ωz = 0.
#choice of time, length, energy units
  t0 = 1.0/ωy
  x0 = sqrt(ħ*t0/m)
  E0 = ħ/t0
#interactions
  #g  = (4*pi*ħ^2*as/m)*x0^3/E0 #dimensionless 3D
  g = 0.1 #test 2D
#damping parameters (dimensionless)
  γ = 0.05
  ℳ  = 0.0
#chemical potential (dimensionless)
  μ  = 12.0
#time evolution parameters
  ti = 0.0
  tf = 3.0/γ  #evolve for 2 damping times
  Nt = 50
  t = collect(linspace(ti,tf,Nt))
  dt = 0.01π/μ #integrate step size [ - should have dt ≪ 2π/μ]
#== end template ==#

@pack sinfo = ωx,ωy,ωz,γ,ℳ,g,t0,x0,E0,μ,ti,tf,Nt,t,dt

#Initialize CField (dimensionless)
  basis = "Hermite"
  ecut = 30*ħ*ωy/E0
  Ω = [ωx; ωy]*t0
  cinfo = makecinfo(ecut,Ω,basis)
  tinfo = maketinfo(cinfo,4)
  @unpack en,P,M = cinfo;
  #x,wx,Tx,y,wy,Ty = makealltrans(M,4,Ω)
  #W = wx.*wy'
#test transform
  c0  = randn(M...)+im*randn(M...); c0=P.*c0
  ψ0  = c2x(c0,tinfo) #initial condition
  ψ = similar(ψ0)  #a field to write to in place

#equation of motion
  function nlin!(dc,c,tinfo)
      c2x!(ψ,c,tinfo)
      x2c!(dc,abs2.(ψ).*ψ,tinfo)
  end

  #in place
  function Lgp!(dc,c,p,t)
      nlin!(dc,c,tinfo)
      dc .= P.*(-im*(1-im*γ)*((en - μ).*c .+ g*dc))
  end

tspan = (t[1],t[end])
prob = ODEProblem(Lgp!,c0,tspan)
alg =    DP5() #abstol = 1e-3, reltol = 1e-3

println("Started evolution ...")
#@time sol = solve(prob,alg,dt=dt,saveat=t,alg_hints=[:stiff],save_everystep=false,dense=false);
@time sol = solve(prob,alg,dt=dt,saveat=t);
println("... Finished.")
return sinfo,cinfo,sol
end

sinfo,cinfo,sol = timeevolution2()

#plot solution for 2D
## Transform to cartesian grid
@unpack ħ,m,ωx,ωy,ωz,γ,ℳ,g,x0,t0,E0,μ,ti,tf,Nt,t,dt = sinfo
#@unpack M,Ω,ecut,P,en = cinfo

#examine solution
using PyPlot

#lets be explicit about units:
Rx = sqrt(2μ*E0/m/ωx^2)
Ry = sqrt(2μ*E0/m/ωy^2)
yMax=1.5Ry
xMax=1.5Rx

Nx = 400
Ny = Nx
x = collect(linspace(-xMax,xMax,Nx))
y = collect(linspace(-yMax,yMax,Ny))
tinfop = maketinfoplot(cinfo,x/x0,y/x0)

#Plot
#using PyPlot
f=figure(figsize=(12,3))
#@manipulate for i=1:length(t) withfig(f,clear=true) do
    ψ = c2x(sol[end],tinfop);
    pcolormesh(x/x0,y/x0,g*abs2.(ψ'))
    xlabel("x/x0")
    ylabel("y/x0")
    colorbar()
    #end
#end


text:
@testset "compute_fvm_tpfa_matrix" begin
    m = 20
    n = 10
    h = 0.1
    kvalue = rand(m*n)
    bc = bcedge("all", m, n, h)
    pval = rand(size(bc, 1))
    K, rhs = compute_fvm_tpfa_matrix(kvalue,bc,pval,m,n,h)
    K1, rhs1 = compute_fvm_tpfa_matrix(constant(kvalue),bc,pval,m,n,h)
    K2 = run(sess, K1)
    rhs2 = run(sess, rhs1)
    @test rhs≈rhs2
    @test K ≈ K2

    m = 20
    n = 10
    h = 0.1
    kvalue = rand(m*n)
    K = compute_fvm_tpfa_matrix(kvalue,m,n,h)
    K1 = compute_fvm_tpfa_matrix(constant(kvalue),m,n,h)
    K2 = run(sess, K1)
    @test K ≈ K2
end

function advection_jl(m, n, h, v, u)
    # v: velocity, (m*n) * 2
    # u: temperature, (m*n)
    a = zeros(m*n)
    for j=1:n
        for i=1:m
            vx = v[(j-1)*m+i, 1]
            vy = v[(j-1)*m+i, 2]
            if vx >= 0
                if i==1
                    a[(j-1)*m+i] += vx * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vx * h * (u[(j-1)*m+i] - u[(j-1)*m+i-1])
                end
            else
                if i==m
                    a[(j-1)*m+i] += - vx * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vx * h * (u[(j-1)*m+i+1] - u[(j-1)*m+i])
                end
            end

            if vy >= 0
                if j==1
                    a[(j-1)*m+i] += vy * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vy * h * (u[(j-1)*m+i] - u[(j-2)*m+i])
                end
            else
                if j==n
                    a[(j-1)*m+i] += - vy * 2 * h * u[(j-1)*m+i]

                else
                    a[(j-1)*m+i] += vy * h * (u[j*m+i] - u[(j-1)*m+i])
                end
            end
        end
    end
    return a
end

@testset "compute_fvm_advection_term" begin
    m = 10
    n = 10
    h = 0.1
    v = rand(m*n,2)
    u = rand(m*n)
    jl_u = advection_jl(m, n, h, v, u)
    u = compute_fvm_advection_term(v,u,m,n,h)
    @show run(sess, u)≈jl_u
end

@testset "grad on gauss points" begin
    m = 10
    n = 10
    h = 0.1
    u = rand(2*(m+1)*(n+1))
    ref = eval_grad_on_gauss_pts(u, m, n, h)
    out = eval_grad_on_gauss_pts(constant(u), m, n, h)
    @test run(sess, out)≈ref
end

@testset "fem_to_gauss_points2" begin
    m = 10
    n = 10
    h = 0.1
    xy = fem_nodes(m, n, h)
    x, y = xy[:,1], xy[:,2]
    u0 = @. x*y*(1-x^2)
    v0 = fem_to_gauss_points(u0, m, n, h)
    @test run(sess, fem_to_gauss_points(constant(u0), m, n, h)) ≈ v0

    figure(figsize=(10,4))
    subplot(121)
    visualize_scalar_on_fem_points(u0, m, n, h)
    subplot(122)
    visualize_scalar_on_gauss_points(v0, m, n, h)
end

@testset "eval_grad_on_gauss_pts" begin
    m = 10
    n = 10
    h = 0.1
    xy = fem_nodes(m, n, h)
    x, y = xy[:,1], xy[:,2]
    u0 = @. x*y*(1-x^2)
    u1 = @. sin(x) * cos(y)
    G = eval_grad_on_gauss_pts([u0;u1], m, n, h)

    xy = gauss_nodes(m, n, h)
    x, y = xy[:,1], xy[:,2]
    g11 = @. y-3x^2*y
    g12 = @. x-x^3
    g21 = @. cos(x) * cos(y)
    g22 = @. -sin(x) * sin(y)

    figure(figsize=(10,10))
    subplot(221)
    visualize_scalar_on_gauss_points(G[:,1,1], m, n, h)
    subplot(222)
    visualize_scalar_on_gauss_points(G[:,1,2], m, n, h)
    subplot(223)
    visualize_scalar_on_gauss_points(G[:,2,1], m, n, h)
    subplot(224)
    visualize_scalar_on_gauss_points(G[:,2,2], m, n, h)

    figure(figsize=(10,10))
    subplot(221)
    visualize_scalar_on_gauss_points(g11, m, n, h)
    subplot(222)
    visualize_scalar_on_gauss_points(g12, m, n, h)
    subplot(223)
    visualize_scalar_on_gauss_points(g21, m, n, h)
    subplot(224)
    visualize_scalar_on_gauss_points(g22, m, n, h)
end

@testset "compute_fem_laplace_term1" begin
    m = 10
    n = 20
    h = 0.1
    u = rand((m+1)*(n+1))
    kappa = rand(4*m*n)
    out = compute_fem_laplace_term1(constant(u), kappa, m, n, h)
    M = compute_fem_laplace_matrix1(kappa, m, n, h)
    out_ = M * u
    @test run(sess, out)≈out_
end

@testset "compute_interaction_term" begin
    m = 10
    n = 20
    h = 0.1
    p = rand(m*n)
    t1 = compute_interaction_term(p, m, n, h)
    t2 = compute_interaction_term(constant(p), m, n, h)
    @test run(sess, t2)≈t1
end

ext: jl


text:
@testset "AIM selection" begin
    # download data
    vcffile = "chr22.1kg.phase3.v5a.vcf.gz"
    isfile(vcffile) ||
        download("http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.vcf/chr22.1kg.phase3.v5a.vcf.gz",
        joinpath(pwd(), vcffile))

    # cd to test folder and read population origin
    joinpath(pathof(VCFTools), "test")
    df = CSV.read("1000genomes.population.txt", DataFrame)

    # create dictionary with key = ID, value = population
    sampleID_to_population = Dict{String, String}()
    for (id, population) in eachrow(df)
        sampleID_to_population[id] = population
    end
    sampleID_to_population

    pvals = VCFTools.aim_select(vcffile, sampleID_to_population)
    @test all(0.0 .≤ pvals .≤ 1.0)
    @test length(pvals) == 424147
    @test pvals[1] == 1.0
    @test pvals[11] ≈ 5.684678238417218e-120
end


text:
# Delving into parsed DXF files.

export walk, dxffind

"""
    walk(action, object::DXFObject; path=[])
Recursively invoke action on `object and each of its components.
Path is a vector of the `DXFObjects` that contain `object` starting
with the subject of the outermost call to `walk`, and not including
object`itself.
Action is called with both `object` and `path` as arguments.  It's
return value is ignored.
"""
function walk end

WalkPathType = Vector{DXFObject}

function walk(action, o::Nothing; path=WalkPathType())
end

function walk(action, o::DXFGroup; path=WalkPathType())
    action(o, path)
end

function walk(action, o::DXFContentsObject; path=WalkPathType())
    action(o, path)
    for child in o
        walk(action, child; path=[path..., o])
    end
end

function walk(action, o::DocumentStart; path=WalkPathType())
    action(o, path)
end

function walk(action, o::HeaderVariable; path=WalkPathType())
    action(o, path)
    path = (path..., o)
    walk(action, o.name; path=path)
    walk(action, o.value; path=path)
end

function walk(action, o::DXFPoint; path=WalkPathType())
    action(o, path)
    path = (path..., o)
    walk(action, o.pointX; path=path)
    walk(action, o.pointY; path=path)
    walk(action, o.pointZ; path=path)
end

function dxffind(predicate, o::DXFObject)
    found = []
    walk(o) do o, path
        if predicate(o)
            push!(found, (path..., o))
        end
    end
    found
end


text:
#以一个部件在所有时层的状况为一个子图，约束主要是子图间的空间关系
#也就是部件与部件的关系
using Plasmo, Ipopt
using PlasmoTest

G = OptiGraph()         #建立主图
O = Ipopt.Optimizer     #求解器为Ipopt

sub1 = OptiGraph()      #建立子图
sub2 = OptiGraph()
add_subgraph!(G,sub1)   #把子图加到主图中
add_subgraph!(G,sub2)

@optinode(sub1,ng1[1:3])    #子图有1-3三个结点
@optinode(sub2,ng2[1:3])

for node in getnodes(sub1)      #批量添加结点，这些结点内部数学表达是相同的
    @variable(node,0 <= x <= 2)
    @variable(node,2 <= y <= 4)
    @constraint(node,x + y >= 5)
end

for node in getnodes(sub2)
    @variable(node,x)
    @NLnodeconstraint(node,exp(x) >= 2)
    @variable(node,y >= 0)
    @constraint(node,x + y >= 3)
end

#添加约束，“ng?[x]”表示 “第?个子图的第x个结点”
@linkconstraint(G,ng1[1][:x]==ng2[1][:y])
@linkconstraint(G,ng2[1][:x]==ng1[2][:y]+1)
@linkconstraint(G,ng2[1][:y]==ng2[2][:y])
@linkconstraint(G,ng2[2][:x]==ng1[3][:y])
@linkconstraint(G,ng2[2][:y]==ng2[3][:y])

#添加目标函数
@objective(G, Min, ng2[1][:x] + ng2[2][:x] + ng2[3][:x])
optimize!(G,O)          #G和O都在最开始有定义
println("Objective = ", digit5(objective_value(G)))

#把图的关系保存成图片
#plt_G = plt_G = Plots.plot(G,node_labels = true, markersize=60, labelsize=30, linewidth=4, layout_options = Dict(:tol => 0.01, :iterations => 2))
#savefig("006.png")

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
"""
    NamedColor(name::String,color::RGB)

A named `ColorTypes.RGB`

# Examples

```julia-repl
julia> NamedColor("DarkGreen")
NamedColor("DarkGreen",(0.0, 0.39, 0.0))
```
"""
struct NamedColor
    name::String
    color::ColorTypes.RGB
end

NamedColor(;name="",color=ColorTypes.RGB(0,0,0)) =
                                NamedColor(name,color)

NamedColor(S::String) = NamedColor(name=S,color=parse(Colors.Colorant,S))
NamedColor(S::Symbol) = NamedColor(string(S))
NamedColor(C::ColorTypes.RGB) = NamedColor(color=C)
function NamedColor(S::String,C::Tuple{<:Real,<:Real,<:Real})
    NamedColor(S,ColorTypes.RGB(C...))
end

function NamedColor(C::Tuple{<:Real,<:Real,<:Real})
    NamedColor("",ColorTypes.RGB(C...))
end

NamedColor(r::Real,g::Real,b::Real) = NamedColor((r,g,b))

softblue = NamedColor(0.92,0.95,0.99);
softred = NamedColor(0.99, 0.92, 0.91);
softyellow = NamedColor(0.98, 0.98, 0.9);
softgreen = NamedColor(0.96, 0.995, 0.98);

function *(k::Real,C::NamedColor)
    NamedColor("",ColorTypes.RGB(k*C.color.r,k*C.color.g,k*C.color.b))
end
dot(k::Real,C::NamedColor) = k*C
zero(::Type{NamedColor}) = NamedColor(0,0,0)
/(C::NamedColor,r::Real) = 1/r*C
function +(C::NamedColor,D::NamedColor)
    NamedColor("",ColorTypes.RGB(C.color.r + D.color.r,
                                 C.color.g + D.color.g,
                                 C.color.b + D.color.b))
end

function ==(C::NamedColor,D::NamedColor)
    C.name == D.name && C.color == D.color
end

function Base.show(io::IO,NC::NamedColor)
    t = map(x->round(Float64(x);digits=2),(NC.color.r,NC.color.g,NC.color.b))
    print(io,"""NamedColor("$(NC.name)",$t)""")
end

defaultcolor() = NamedColor(name="default")

"""
    Pen(color::NamedColor,
        opacity::Real,
        linewidth::Real,
        other::String)

Store drawing properties

# Examples

```julia-repl
julia> Pen(color="DarkGreen",opacity=0.5)
Pen(color=DarkGreen,opacity=0.5)
```
"""
struct Pen
    color::NamedColor
    opacity::Real
    linewidth::Real
    fontsize::Real
    other::String
end

const _DEFAULT_PEN_ARGS =
    OrderedDict(
        :color => defaultcolor(),
        :opacity => 1,
        :linewidth => 0,
        :fontsize => 12,
        :other => ""
    )

Pen(;kwargs...) = Pen(updatedvals(_DEFAULT_PEN_ARGS,kwargs)...)

NoPen() = Pen(other="NoPen")

==(P::Pen,Q::Pen) = all(getfield(P,k)==getfield(Q,k)
                                for k in keys(_DEFAULT_PEN_ARGS))

convert(::Type{NamedColor},name::Union{Symbol,AbstractString}) = NamedColor(name)
convert(::Type{NamedColor},color::ColorTypes.RGB) = NamedColor(;color=color)
convert(::Type{Pen},string::AbstractString) = Pen(;other=string)
convert(::Type{Pen},color::NamedColor) = Pen(color=color)

function Base.show(io::IO,P::Pen)
    kwargs = kwargstring(P,_DEFAULT_PEN_ARGS,semicolon=false)
    print(io,"Pen($kwargs)")
end

isdefault(P::Pen) = (P == Pen())
isdefault(S::AbstractString) = (S == "")

"""
    Arrow(name::AbstractString,
          size::Real,
          position::Real)

Store instructions for drawing an arrow

`position` is an element of [0,1] which indicates
how far along the path the arrow should be drawn

NoArrow() returns a no-arrow instruction, while
Arrow3() gives an arrow suitable for 3D paths
"""
struct Arrow
    name::AbstractString
    size::Real
    position::Real
end

const _DEFAULT_ARROW_KWARGS =
    OrderedDict(
        :name => "Arrow",
        :size => 7.5,
        :position => 1
    )

function Arrow(;kwargs...)
    return Arrow(updatedvals(_DEFAULT_ARROW_KWARGS,kwargs)...)
end

function Arrow3(;kwargs...)
    D = deepcopy(_DEFAULT_ARROW_KWARGS)
    merge!(D,Dict(kwargs))
    merge!(D,Dict(:name => "Arrow3"))
    return Arrow(values(D)...)
end

NoArrow() = Arrow(;name="None")
Arrow(n) = Arrow(size=n)
Arrow3(n) = Arrow3(size=n)
Arrows(;kwargs...) = Arrow(name="Arrows";kwargs...)

==(A::Arrow,B::Arrow) = all(getfield(A,k) == getfield(B,k)
                                    for k in fieldnames(Arrow))

function AsyString(A::Arrow)
    position = A.position == 1 ? "" : "position=$(A.position)"
    if A.name == "None"
        return "None"
    elseif A.size ≠ 7.5
        return "$(A.name)($(filterjoin(A.size,position)))"
    else
        return "$(A.name)($position)"
    end
end

function string(A::Arrow)
    if A.name == "None"
        return "NoArrow()"
    else
        return AsyString(A)
    end
end
print(io::IO,A::Arrow) = print(io,string(A))
show(io::IO,A::Arrow) = print(io,string(A))

function string(C::NamedColor)
    if C.name ≠ "default" && C.name ≠ ""
        return C.name
    elseif C.name == "default"
        return ""
    else
        colorvals = [Float64(getfield(C.color,s)) for s in [:r,:g,:b]]
        return """rgb($(join(map(string,colorvals),",")))"""
    end
end
print(io::IO,C::NamedColor) = print(io,string(C))

function string(P::Pen)
    color = string(P.color)
    opacity = P.opacity == 1 ? "" : "opacity($(P.opacity))"
    linewidth = P.linewidth == 0 ? "" : "linewidth($(P.linewidth))"
    fontsize = P.fontsize == 12 ? "" : "fontsize($(P.fontsize))"
    other = P.other
    join(filter(x->x≠"",[color,opacity,linewidth,fontsize,other]),"+")
end
print(io::IO,P::Pen) = print(io,string(P))

function Ticks(major::Union{Vector,AbstractRange},minor::Union{Vector,AbstractRange,Nothing};beginlabel=true,endlabel=true)
    majorticks = "Ticks=new real[] {$(join(major,","))}"
    minorticks = isa(minor,Nothing) ? "" : "ticks=new real[] {$(join(minor,","))}"
    bl = beginlabel ? "" : "beginlabel=false"
    el = endlabel ? "" : "endlabel=false"
    args = filterjoin(majorticks,minorticks,bl,el)
    "Ticks($args)"
end

function Ticks(major::Union{Vector,AbstractRange};kwargs...)
    Ticks(major,nothing;kwargs...)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4a40583e-5993-11eb-3121-0319695416d3
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
end

# ╔═╡ 98ff088e-5994-11eb-1b32-d928c49e8466
using Flux

# ╔═╡ a0b5cd0e-5994-11eb-2572-8334d710e797
using Plots

# ╔═╡ cd235348-5994-11eb-09f1-ad1bd636ccc4
using Statistics

# ╔═╡ 52a3fc2a-5992-11eb-27b4-0da25f26e08d
md"""
# Chapter 8.2 Dyna: Integrated Planning, Acting, and Learning

To demonstrate the flexibility of `ReinforcementLearning.jl`, the `DynaAgent` is also included and we'll explore its performance in this notebook.
"""

# ╔═╡ 50dea4e0-5993-11eb-1cd3-c3e6a79800a4
md"""
## The Maze Environment

In this chapter, the authors introduced a specific maze environment. So let's define it by implementing the interfaces in `ReinforcementLearning.jl`.
"""

# ╔═╡ a75c916a-5993-11eb-1928-ef4651464d8e
const LRUD = [
    CartesianIndex(0, -1),  # left
    CartesianIndex(0, 1),   # right
    CartesianIndex(-1, 0),  # up
    CartesianIndex(1, 0),   # down
]

# ╔═╡ b8c9bd9a-5993-11eb-3acc-6532b72f9bb2
begin
	mutable struct MazeEnv <: AbstractEnv
		walls::Set{CartesianIndex{2}}
		position::CartesianIndex{2}
		start::CartesianIndex{2}
		goal::CartesianIndex{2}
		NX::Int
		NY::Int
	end

	function MazeEnv()
		walls = Set([
			[CartesianIndex(i, 3) for i = 2:4]
			CartesianIndex(5, 6)
			[CartesianIndex(j, 8) for j = 1:3]
		])
		start = CartesianIndex(3, 1)
		goal = CartesianIndex(1, 9)
		MazeEnv(walls, start, start, goal, 6, 9)
	end

	function (env::MazeEnv)(a::Int)
		p = env.position + LRUD[a]
		if p == env.goal
			env.position = env.goal
		elseif !(p ∈ env.walls)
			env.position = CartesianIndex(min(max(p[1], 1), env.NX), min(max(p[2], 1), env.NY))
		end
		nothing
	end
end

# ╔═╡ c289a410-5993-11eb-2601-cdb3102ee32e
RLBase.action_space(env::MazeEnv) = Base.OneTo(length(LRUD))

# ╔═╡ 35054552-5994-11eb-0541-47fc7dba1724
RLBase.reward(env::MazeEnv) = Float64(env.position == env.goal)

# ╔═╡ 46a51362-5994-11eb-3645-87be0a146644
RLBase.is_terminated(env::MazeEnv) = env.position == env.goal

# ╔═╡ 494d0848-5994-11eb-0845-857771b8ec2d
RLBase.reset!(env::MazeEnv) = env.position = env.start

# ╔═╡ 0b6f59b2-5994-11eb-2255-b96582fb8d2e
begin
	import Base: *

	function extend(p::CartesianIndex{2}, n::Int)
		x, y = Tuple(p)
		[CartesianIndex(n * (x - 1) + i, n * (y - 1) + j) for i = 1:n for j = 1:n]
	end

	function remap(p::CartesianIndex{2}, n::Int)
		x, y = Tuple(p)
		CartesianIndex((x - 1) * n + 1, (y - 1) * n + 1)
	end

	function *(env::MazeEnv, n::Int)
		walls = Set{CartesianIndex{2}}(ww for w in env.walls for ww in extend(w, n))
		start, position, goal = remap(env.start, n), remap(env.position, n), remap(env.goal, n)
		NX, NY = env.NX * n, env.NY * n
		MazeEnv(walls, position, start, goal, NX, NY)
	end
end

# ╔═╡ bdabb1bc-5993-11eb-2b28-2d49b0e98b68
RLBase.state_space(env::MazeEnv) = Base.OneTo(env.NX * env.NY)

# ╔═╡ 4422f2e6-5994-11eb-1e9e-33ffb5cbcbdb
RLBase.state(env::MazeEnv) = (env.position[2] - 1) * env.NX + env.position[1]

# ╔═╡ 524b4b96-5994-11eb-07c9-212ded5607d7
x = MazeEnv()

# ╔═╡ 5d4f1c16-5994-11eb-1206-0bb9f1b95d70
md"""
## Figure 8.2
"""

# ╔═╡ 77c8c1e6-5994-11eb-03d8-85cd80e3f4b0
function plan_step(n)
	env = MazeEnv()
	ns = length(state_space(env))
	na = length(action_space(env))

	agent = DynaAgent(
		policy=QBasedPolicy(
			learner=TDLearner(
				approximator=TabularQApproximator(
					n_state=ns,
					n_action=na,
					opt=Descent(0.1)
					),
				γ=0.1,
				method=:SARS
				),
			explorer=EpsilonGreedyExplorer(0.1;is_break_tie=true)
			),
		model=ExperienceBasedSamplingModel(),
		trajectory=VectorSARTTrajectory(),
		plan_step=n
	)

	hook = StepsPerEpisode()
	run(agent, env, StopAfterEpisode(50),hook)
	hook.steps
end

# ╔═╡ c85451a2-5994-11eb-1fa3-c1f9fc128a3f
begin
	fig_8_2 = plot(legend=:topright, xlabel="Episodes", ylabel="Steps per episode")
	for n in [0, 5, 50]
		plot!(fig_8_2, mean(plan_step(n) for _ in 1:30), label="plan_step = $n")
	end
	fig_8_2
end

# ╔═╡ e5b3dbfe-5995-11eb-3c89-659c945e3d79
md"""
## Figure 8.4
"""

# ╔═╡ 19c4c11c-5996-11eb-1d02-4b044c55ed30
function cumulative_dyna_reward(model, walls, nstep1, change, nstep2)
    env = MazeEnv(
		walls,
		CartesianIndex(6, 4),
		CartesianIndex(6, 4),
		CartesianIndex(1, 9),
		6,
		9
	)
	ns = length(state_space(env))
	na = length(action_space(env))
    agent = DynaAgent(
        policy=QBasedPolicy(
            learner=TDLearner(
                approximator=TabularQApproximator(
					n_state=ns,
					n_action=na,
					opt=Descent(1.)
					),
                γ=0.95,
                method=:SARS
                ),
            explorer=EpsilonGreedyExplorer(0.1;is_break_tie=true)
            ),
        model=model,
        trajectory=VectorSARTTrajectory(),
        plan_step=10
	)

    hook = StepsPerEpisode()
    run(agent, env, StopAfterStep(nstep1;is_show_progress=false),hook)
    change(env.walls)
    run(agent, env, StopAfterStep(nstep2;is_show_progress=false),hook)

	cumulative_reward = []
	for (i, n) in enumerate(hook.steps)
		for _ in 1:n
			push!(cumulative_reward, i)
		end
	end
	for _ in (nstep1+nstep2):-1:length(cumulative_reward)
		push!(cumulative_reward, length(hook.steps))
	end
	cumulative_reward
end

# ╔═╡ 3e68d0a0-5997-11eb-28b0-9fe9f480212e
walls() = Set([CartesianIndex(4, j) for j in 1:8])

# ╔═╡ 5040c99a-5997-11eb-3b1f-ff6a7e7fc122
function change_walls(walls)
    pop!(walls, CartesianIndex(4,1))
    push!(walls, CartesianIndex(4,9))
end

# ╔═╡ d3463894-59a0-11eb-2bcf-0dc30ae600dc
begin
	fig_8_4 = plot(legend=:topleft, xlabel="Time steps", ylabel="Cumulative reward")
	plot!(fig_8_4, mean(cumulative_dyna_reward(ExperienceBasedSamplingModel(), walls(), 1000, change_walls, 2000) for _ in 1:30), label="Dyna-Q")
	plot!(fig_8_4, mean(cumulative_dyna_reward(TimeBasedSamplingModel(;n_actions=4), walls(), 1000, change_walls, 2000) for _ in 1:30), label="Dyna-Q+")
	fig_8_4
end

# ╔═╡ af1b3f58-59a2-11eb-2cdc-073dd3f3fbc4
md"""
## Figure 8.5
"""

# ╔═╡ 3b4c6844-59a1-11eb-3e95-8d9ebfbe393b
new_walls() = Set([CartesianIndex(4, j) for j in 2:9])

# ╔═╡ 4fb4ac62-59a1-11eb-3fe1-ebeea8006edb
function new_change_walls(walls)
    pop!(walls, CartesianIndex(4,9))
end

# ╔═╡ b1d8fc48-59a1-11eb-27f0-cdc20954dbb3
begin
	fig_8_5 = plot(legend=:topleft, ylabel="Cumulative reward", xlabel="Time steps")
	plot!(fig_8_5, mean(cumulative_dyna_reward(ExperienceBasedSamplingModel(), new_walls(), 3000, new_change_walls, 3000) for _ in 1:50), label="dyna-Q")
	plot!(fig_8_5, mean(cumulative_dyna_reward(TimeBasedSamplingModel(n_actions=4, κ = 1e-3), new_walls(), 3000, new_change_walls, 3000) for _ in 1:50), label="dyna-Q+")
	fig_8_5
end

# ╔═╡ bb2fa194-59a2-11eb-1b5b-39fe96f94dd9
md"""
## Example 8.4
"""

# ╔═╡ 23b5f902-59a3-11eb-1af3-09163b251737
function run_once(model, ratio=1)
    env = MazeEnv() * ratio
    ns = length(state_space(env))
    na = length(action_space(env))
    agent = DynaAgent(
        policy=QBasedPolicy(
            learner=TDLearner(
                approximator=TabularQApproximator(
					n_state=ns,
					n_action=na,
					opt=Descent(0.5)),
                γ=0.95,
                method=:SARS
                ),
            explorer=EpsilonGreedyExplorer(0.1;is_break_tie=true)
            ),
        model=model,
        trajectory=VectorSARTTrajectory(),
        plan_step=5
    )
    hook = StepsPerEpisode()
    run(agent, env, (args...) -> length(hook.steps) > 0 && hook.steps[end] <= 14 * ratio * 1.2,hook)
    model.sample_count
end

# ╔═╡ 721bf66c-59a3-11eb-3dc2-b965593cf427
begin
	p = plot(legend=:topleft)
	plot!(mean([run_once(ExperienceBasedSamplingModel(), ratio) for ratio in 1:6] for _ in 1:5), label="Dyna", yscale=:log10)
	plot!(mean([run_once(PrioritizedSweepingSamplingModel(), ratio) for ratio in 1:6] for _ in 1:5), label="Prioritized", yscale=:log10)
	p
end

# ╔═╡ Cell order:
# ╟─52a3fc2a-5992-11eb-27b4-0da25f26e08d
# ╠═4a40583e-5993-11eb-3121-0319695416d3
# ╟─50dea4e0-5993-11eb-1cd3-c3e6a79800a4
# ╠═a75c916a-5993-11eb-1928-ef4651464d8e
# ╠═b8c9bd9a-5993-11eb-3acc-6532b72f9bb2
# ╠═bdabb1bc-5993-11eb-2b28-2d49b0e98b68
# ╠═c289a410-5993-11eb-2601-cdb3102ee32e
# ╠═35054552-5994-11eb-0541-47fc7dba1724
# ╠═4422f2e6-5994-11eb-1e9e-33ffb5cbcbdb
# ╠═46a51362-5994-11eb-3645-87be0a146644
# ╠═494d0848-5994-11eb-0845-857771b8ec2d
# ╠═0b6f59b2-5994-11eb-2255-b96582fb8d2e
# ╠═524b4b96-5994-11eb-07c9-212ded5607d7
# ╟─5d4f1c16-5994-11eb-1206-0bb9f1b95d70
# ╠═98ff088e-5994-11eb-1b32-d928c49e8466
# ╠═77c8c1e6-5994-11eb-03d8-85cd80e3f4b0
# ╠═a0b5cd0e-5994-11eb-2572-8334d710e797
# ╠═cd235348-5994-11eb-09f1-ad1bd636ccc4
# ╠═c85451a2-5994-11eb-1fa3-c1f9fc128a3f
# ╟─e5b3dbfe-5995-11eb-3c89-659c945e3d79
# ╠═19c4c11c-5996-11eb-1d02-4b044c55ed30
# ╠═3e68d0a0-5997-11eb-28b0-9fe9f480212e
# ╠═5040c99a-5997-11eb-3b1f-ff6a7e7fc122
# ╠═d3463894-59a0-11eb-2bcf-0dc30ae600dc
# ╟─af1b3f58-59a2-11eb-2cdc-073dd3f3fbc4
# ╠═3b4c6844-59a1-11eb-3e95-8d9ebfbe393b
# ╠═4fb4ac62-59a1-11eb-3fe1-ebeea8006edb
# ╠═b1d8fc48-59a1-11eb-27f0-cdc20954dbb3
# ╟─bb2fa194-59a2-11eb-1b5b-39fe96f94dd9
# ╠═23b5f902-59a3-11eb-1af3-09163b251737
# ╠═721bf66c-59a3-11eb-3dc2-b965593cf427


text:
"""
    LatentGP(f<:AbstractGP, lik, Σy)

 - `f` is a `AbstractGP`.
 - `lik` is the likelihood function which maps samples from `f` to the
   corresponding conditional likelihood distributions (i.e., `lik` must return a
   `Distribution` compatible with the observations).
 - `Σy` is the noise under which the latent GP is "observed"; this represents
   the jitter used to avoid numeric instability and should generally be small.
"""
struct LatentGP{Tf<:AbstractGP,Tlik,TΣy}
    f::Tf
    lik::Tlik
    Σy::TΣy
end

"""
    LatentFiniteGP(fx<:FiniteGP, lik)

 - `fx` is a `FiniteGP`.
 - `lik` is the likelihood function which maps samples from `f` to the
   corresponding conditional likelihood distributions (i.e., `lik` must return a
   `Distribution` compatible with the observations).
"""
struct LatentFiniteGP{Tfx<:FiniteGP,Tlik}
    fx::Tfx
    lik::Tlik
end

(lgp::LatentGP)(x) = LatentFiniteGP(lgp.f(x, lgp.Σy), lgp.lik)

Base.length(lgpx::LatentFiniteGP) = length(lgpx.fx)

function Distributions.rand(rng::AbstractRNG, lfgp::LatentFiniteGP)
    f = rand(rng, lfgp.fx)
    y = rand(rng, lfgp.lik(f))
    return (f=f, y=y)
end

"""
    logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})

```math
    log p(y, f; x)
```
The joint log density of the Gaussian process output `f` and observation `y`.
"""
function Distributions.logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})
    return logpdf(lfgp.fx, y.f) + logpdf(lfgp.lik(y.f), y.y)
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using LinearAlgebra

T = SymTridiagonal(2ones(5), -ones(4))

T * randn(5)     # Matrix-vector multiplication
T^3              # Matrix cube

T + 2I     # OK, since T is 3-by-3

I4 = Matrix{Float64}(I, 4, 4)


text:
module SSLTestIteration
using SimpleSDMLayers
using Test

S = SimpleSDMPredictor([1 2 3; 4 5 6; 7 nothing 9], 0.0, 1.0, 0.0, 1.0)

val = sort([v for v in S])
@test all(val .== vec([1 2 3 4 5 6 7 9]))
@test all(val .== sort(collect(S)))

end


text:
if basename(pwd()) == "aoc"
    cd("2018/11")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using GraphMatFun, LsqFit, LinearAlgebra, Random
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=8
target_n=400 # Measure error wrt this discr
n=100;        # Optimize wrt this discr
rho=13.5;  # Increased SID rho m=7

its=5;
droptol0=1e-20;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its,
                :errtype => :relerr)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;
base.params[:kickit_mode]=1;
base.params[:cref_count]=0;

oldrho="6_4";
oldrho2="3_59";
## First the optimization free
#filename="simulations/newgraphs/exp_m$(m)_sid_$oldrho.cgr";
#sid_org=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)
#
##filename="simulations/newgraphs/exp_m$(m)_bbc_$oldrho.cgr";
##bbc_org=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)
#
#filename="simulations/newgraphs/exp_m$(m)_sastre_$oldrho.cgr";
#sastre_org=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true)
#
#
##sastre_org.cref=get_degopt_crefs(sastre_org.graph);
#
#
#filename="simulations/newgraphs/exp_m$(m)_ps_$oldrho.cgr";
#ps_org=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)
#
#
#filename="simulations/newgraphs/exp_m$(m)_mono_$oldrho.cgr";
#mono_org=init_state_file!(deepcopy(base),:mono,filename,showmeta=true)
#

#filename="simulations/newgraphs/exp_m$(m-1)_mono+GN_$oldrho.cgr";
##filename="simulations/graphs/exp_m8_mono_taylor_13_5.cgr";
#mono0=init_state_file!(deepcopy(base),:mono,filename,showmeta=true,scale_and_square=true)
#mono0.cref=mono0.cref[8:70]
#(mono,simlist,commandlist)=
#     run_sequence(mono0,"ksssssdddsssssddddssssssddddddsssssssdddsssssssssssdddssssGGsssddddggggggssssssssGGGsssGGGssssdddddggggggggsssssssGGGsGGGGGssssdddddggggggggssssssssGGGsGGGGGsssssddddgggggggsssGGGGGGGssssssdddggggggggggssssssggssssssGsGsGsGsGGsGGsGGGsGGsss");
#

datadir=joinpath("..","..","data","exp");
#filename="simulations/graphs/exp_m8_mono_taylor_13_5.cgr";
filename=joinpath(datadir,"exp_mono_m$(m-1)_opt_rho$(oldrho).cgr");
mono0=init_state_file!(deepcopy(base),:mono,filename,showmeta=true,
                       scale_and_square=true)

degopt_crefs=get_degopt_crefs(mono0.graph);
#mono0.cref=mono0.cref[8:end-8]
mono0.cref=[#get_degopt_crefs(mono0.graph)[1][end-1][1];
            #get_degopt_crefs(mono0.graph)[1][end-1][2];
            degopt_crefs[1][end][1];
            degopt_crefs[1][end][2];
            degopt_crefs[2] ];
(mono,simlist,commandlist)=
run_sequence(deepcopy(mono0),"sSSNSkSSGGGSS22SdddddSSSSSzS0pSZpSSSSzSzq");

filename=joinpath(datadir,"exp_ps_m$(m-1)_opt_rho$(oldrho).cgr");
ps0=init_state_file!(deepcopy(base),:ps,filename,showmeta=true,
                         scale_and_square=true)
ps0.cref=[degopt_crefs[1][end][1];
            degopt_crefs[1][end][2];
            degopt_crefs[2]];

(ps,simlist,commandlist)=
     run_sequence(ps0,
                  "SSNSSSSGGGSSSSSDDDSSSSZSSSZSSSddSddSzZSSSSddZSSSSZpSSSzq");

filename=joinpath(datadir,"exp_sastre_m$(m-1)_opt_rho$(oldrho).cgr");
sastre0=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true,
                         scale_and_square=true)
sastre0.cref=[degopt_crefs[1][end][1];
            degopt_crefs[1][end][2];
            degopt_crefs[2]];

(sastre,simlist,commandlist)=
     run_sequence(sastre0,
                  "SSSSNSSSSkSSGGGSSddddSdSdSzSkZSSSSSSSZsSSzZSSSSpZpSSSSddSddddSkSSS2zq");

#filename=joinpath(datadir,"exp_sid_m$(m-1)_opt_rho$(oldrho).cgr");
#sid0=init_state_file!(deepcopy(base),:sid,filename,showmeta=true,
#                         scale_and_square=true)
#sid0.cref=[degopt_crefs[1][end][1];
#            degopt_crefs[1][end][2];
#            degopt_crefs[2]];
#
#(sid,simlist,commandlist)=
#     run_sequence(sid0,
#                  "SSSSNSSSSkSSGGGSSddddSdSdSzSkZSSSSSSSZsSSzZSSSSpZpSSSSddSddddSkSSS2zq");
#
#asd
#sSSNSkSSGGGSS22Sq");
##sSkSSSddddSSSSSSS22GGG22ddd2SS22n222tk22222222223tdddtSS333k333333333dddd3N333k3333333k3333333k3333333dddds333333333333333333333sSSStztkk3333N3nn3333n33NN333n33tz33NNN3333N3tS333ts333t33Skkk3333333tTz");
#
#graph1=mono1.graph;
##degopt=Degopt(mono1.graph);
###degopt.x[1][1] += sqrt(eps());;
##degopt.x[1][1][1] += 0.000001;
##normalize!(degopt);
##degopt.x[1][1][1] -= 0.000001;
##graph1=graph_degopt(degopt)[1];
#export_compgraph(graph1,"/tmp/mono1_m8_tmp.cgr");
#
#
#mono0.cref=[#get_degopt_crefs(mono0.graph)[1][end-1][1];
#            #get_degopt_crefs(mono0.graph)[1][end-1][2];
#            get_degopt_crefs(mono0.graph)[1][end][1];
#            get_degopt_crefs(mono0.graph)[1][end][2];
#            get_degopt_crefs(mono0.graph)[2] ];
#filename_tmp="/tmp/mono1_m8_tmp.cgr";
#mono2=init_state_file!(deepcopy(base),:mono,filename_tmp,showmeta=true);
#setdiff!(mono2.cref,mono0.cref);
#(mono,simlist,commandlist)=
#run_sequence(deepcopy(mono2),"dddddddddddddddddd22GGG22n2S");
#
#
#
#

#(mono,simlist,commandlist)=
#run_sequence(deepcopy(mono0),"ssssGGGssgddssssssddddggggggggggsssGGGsss");

#             "kssssGGGssgddssssssddddggggggggggsssGGGssssGGGssGGGGsssssksssdddsssssddssssssdddsdgggggggsssssGGGGsGGGssssssNssGssddssssddssssdggggggggsssssGGGGssssGGGGsssddsdsdgggggggggggssssGGGGGGGGGGsGssssggggggggdsdsdskssGGGGGGGGssNsszq");
#sGsGsGGGssssssssddsssGGGGGGsssksssssNsq
#             "sssGGGsssdddsddsssssdddsssdddsssssq");

#filename="simulations/newgraphs/exp_m$(m-1)_mono+GN_$(oldrho).cgr";
#mono0=init_state_file!(deepcopy(base),:mono,filename,showmeta=true,scale_and_square=true)
#(mono,simlist,commandlist)=
#     run_sequence(mono0,"ddddddddddddddddddddddddddsssGGGsssdddsddsssssdddsssdddsssssq");
#
#filename="simulations/newgraphs/exp_m$(m)_ps+GN_$oldrho.cgr";
#ps0=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)
#(ps,simlist,commandlist)=
#     run_sequence(ps0,
#     "sq");
#
#filename="simulations/newgraphs/exp_m$(m-1)_sastre+GN_$oldrho.cgr";
#sastre0=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true,scale_and_square=true)
#(sastre,simlist,commandlist)=
#     run_sequence(sastre0,
#                  "kssGGGssddddsdddddsdsssssdddssssddssssddsgggsssssGGGssq");
#
#filename="simulations/newgraphs/exp_m$(m)_sid+GN_$oldrho.cgr";
#sid0=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)
#(sid,simlist,commandlist)=
#     run_sequence(sid_org,
#     "nssNdddddddddddddddddddssGGGsssddsddsddsddsd")
#
##filename="simulations/newgraphs/exp_m$(m)_bbc+GN_$oldrho.cgr";
##bbc0=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)
##(bbc,simlist,commandlist)=
##     run_sequence(bbc0,
##     "sdddddsdddddssssssddddsssssssddddsssssssssssssdddssssssssssddddssssssssdsggggggsdssssssssssssssssGGGGGsssGGGsssq")
##
#
include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using Bukdu

struct WelcomeController <: ApplicationController
    conn::Conn
end

function index(c::WelcomeController)
    render(JSON, "Hello World")
end

if PROGRAM_FILE == basename(@__FILE__)

routes() do
    get("/", WelcomeController, index)
end

Bukdu.start(8080)

# Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module FiniteVolumeMethod

    import FullyThreadedTree
    import FullyThreadedTree: Tree, Face, centroid, at_boundary, at_refinement, regular, active, level, parent_of_active, initialized, cells, faces, levels, face_area, area, volume, cell_volume

    export Tree,
           Face,
           at_boundary,
           at_refinement,
           regular,
           active,
           level,
           parent_of_active,
           initialized,
           centroid

    include("abstract.jl")

    export CellVar,
           FaceVar,
           size

    include("treegrid.jl")

    export TreeGrid,
           index,
           refine!,
           coarsen!

    include("cartesiangrid.jl")

    export CartesianGrid,
           area,
           volume,
           cell_volume

    include("interface.jl")

    export show,
           faces,
           cells

    include("operators.jl")

    export gradient,
           divergence

    include("functionals.jl")

    export integral,
           inner_product,
           norm

    include("plotting.jl")

    export plot

end # module


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# This file contains the routines that interfaces with SpMP library.

const LIB_PATH = libcsr

@doc """
Create a knob for a matrix with the given properties of it, while this matrix
might not have been created yet.

Properties:
constant_valued:        The matrix is a constant in value. It implies
                        constant_structured below.
constant_structured:    The matrix has always the same structure, even if its
                        value may be changed.
is_symmetric:           The matrix is symmetric in value.
                        It implies is_structure_symmetric below.
is_structure_symmetric: The matrix is symmetric in structure.
is_structure_only :     Only the structure of matrix should be used.
is_single_def:          The matrix is statically defined only once.
"""
function new_matrix_knob(
    A                      :: Symbol, # Useless. Only for better readability of the generated code
    constant_valued        = false,
    constant_structured    = false,
    is_symmetric           = false,
    is_structure_symmetric = false,
    is_structure_only      = false,
    is_single_def          = false
 )
    # Note: is_structure_symmetric is only for CoSP2. In general, the matrix
    # must be constant valued or structured.
    # dump(A)
    # @printf("%d %d %d %d\n", constant_valued, constant_structured, is_structure_symmetric, collective_structure_prediction_enabled)
    assert(constant_valued || constant_structured || is_structure_symmetric || !collective_structure_prediction_enabled)
    assert(!constant_valued || constant_structured || !collective_structure_prediction_enabled)
    assert(!is_symmetric || is_structure_symmetric || !collective_structure_prediction_enabled)
    mknob = ccall((:NewMatrixKnob, LIB_PATH), Ptr{Void},
                   (Bool, Bool, Bool, Bool, Bool, Bool),
                   constant_valued, constant_structured, is_symmetric,
                   is_structure_symmetric, is_structure_only, is_single_def)
end

@doc """ Fill matrix info into the knob """
function fill_matrix_info_once(
    mknob :: Ptr{Void},
    A     :: SparseMatrixCSC
 )
    ccall((:FillMatrixInfoOnce, LIB_PATH), Void,
            (Ptr{Void}, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
             mknob, A.m, A.n, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval))
end

@doc """ Set the matrix as constant_structured. """
function set_constant_valued(
    mknob :: Ptr{Void}
)
    ccall((:SetConstantValued, LIB_PATH), Void, (Ptr{Void},), mknob)
end

@doc """ Set the matrix as constant_structured. """
function set_constant_structured(
    mknob :: Ptr{Void}
)
    ccall((:SetConstantStructured, LIB_PATH), Void, (Ptr{Void},), mknob)
end

@doc """ Set the matrix as value symmetric. """
function set_value_symmetric(
    mknob :: Ptr{Void}
)
    ccall((:SetValueSymmetric, LIB_PATH), Void, (Ptr{Void},), mknob)
end

@doc """ Set the matrix as structurally symmetric. """
function set_structure_symmetric(
    mknob :: Ptr{Void}
)
    ccall((:SetStructureSymmetric, LIB_PATH), Void, (Ptr{Void},), mknob)
end

@doc """ set matrix derivative """
function set_derivative(
    mknob           :: Ptr{Void},
    derivative_type :: Int,
    derivative      :: Ptr{Void}
)
    ccall((:SetDerivative, LIB_PATH), Void, (Ptr{Void}, Cint, Ptr{Void}), mknob, derivative_type, derivative)
end

@doc """
Derivative types between two matrix knobs.
NOTE: when adding any new type, or changing any type' integer value,
make sure to change int2derivative_map as well.
"""
const DERIVATIVE_TYPE_TRANSPOSE = 0
const DERIVATIVE_TYPE_SYMMETRIC = 1
const DERIVATIVE_TYPE_LOWER_TRIANGULAR = 2
const DERIVATIVE_TYPE_UPPER_TRIANGULAR = 3

const int2derivative_map = Dict(
    0 => GlobalRef(Sparso, :DERIVATIVE_TYPE_TRANSPOSE),
    1 => GlobalRef(Sparso, :DERIVATIVE_TYPE_SYMMETRIC),
    2 => GlobalRef(Sparso, :DERIVATIVE_TYPE_LOWER_TRIANGULAR),
    3 => GlobalRef(Sparso, :DERIVATIVE_TYPE_UPPER_TRIANGULAR)
)

@doc """ Delete a matrix knob. """
function delete_matrix_knob(
    mknob :: Ptr{Void}
)
    ccall((:DeleteMatrixKnob, LIB_PATH), Void, (Ptr{Void},), mknob)
end

@doc """ Propagate a matrix knob's information to another. """
function propagate_matrix_info(
    to_mknob   :: Ptr{Void},
    from_mknob :: Ptr{Void}
)
    ccall((:PropagateMatrixInfo, LIB_PATH), Void, (Ptr{Void}, Ptr{Void}), to_mknob, from_mknob)
end

@doc """ Associate a matrix knob with a function knob. """
function add_mknob_to_fknob(
    mknob :: Ptr{Void},
    fknob :: Ptr{Void}
)
    ccall((:AddMatrixKnob, LIB_PATH), Void, (Ptr{Void}, Ptr{Void}),
        fknob, mknob)
end

@doc """ Create a function knob. """
function new_function_knob(
    fknob_creator :: AbstractString
)
    assert(fknob_creator != "")

    # Simulate ccall((fknob_creator, LIB_PATH), Ptr{Void}, ()). We cannot directly
    # use this ccall here because (fknob_creator, LIB_PATH) is treated as a tuple,
    # instead of a pointer or expression.
    expr = Expr(:call,
        TopNode(:ccall),
        Expr(:call, TopNode(:tuple), QuoteNode(fknob_creator), LIB_PATH),
        Ptr{Void},
        Expr(:call, TopNode(:svec))
    )
    eval(expr)
end
new_function_knob() = new_function_knob("NewFunctionKnob")

@doc """ Delete a function knob. """
function delete_function_knob(
    fknob_deletor :: AbstractString,
    fknob         :: Ptr{Void}
)
    assert(fknob_deletor != "")

    # Simulate ccall((fknob_deletor, LIB_PATH), Void, (Ptr{Void},), fknob). We cannot directly
    # use this ccall here because (fknob_deletor, LIB_PATH) is treated as a tuple,
    # instead of a pointer or expression.
    expr = Expr(:call,
        TopNode(:ccall),
        Expr(:call, TopNode(:tuple), QuoteNode(fknob_deletor), LIB_PATH),
        GlobalRef(Main, :Void),
        Expr(:call,
                TopNode(:svec),
                Expr(:call,
                      TopNode(:apply_type),
                      GlobalRef(Main, :Ptr),
                      GlobalRef(Main, :Void)
                )
        ),
        Expr(:call,
                TopNode(:unsafe_convert),
                Expr(:call,
                      TopNode(:apply_type),
                      GlobalRef(Main, :Ptr),
                      GlobalRef(Main, :Void)
                ),
                fknob
        ),
        fknob
    )
    eval(expr)
end
delete_function_knob(fknob :: Ptr{Void}) = delete_function_knob("DeleteFunctionKnob", fknob)

function replace_lower_with_UT(
    A     :: SparseMatrixCSC,
    U     :: SparseMatrixCSC
)
    ccall((:ReplaceLowerWithUT, LIB_PATH), Void,
           (
            Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
            Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}
           ),
            A.m, A.n, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval),
            U.m, U.n, pointer(U.colptr), pointer(U.rowval), pointer(U.nzval))
end

@doc """
Context-sensitive forward triangular solver, equivalent to Base.SparseMatrix.
fwdTriSolve!(L, b).
"""
function fwdTriSolve!(
    L     :: SparseMatrixCSC{Float64, Int32},
    b     :: Vector,
    fknob :: Ptr{Void}
 )
    if fknob != C_NULL
        mknob_L = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
        assert(mknob_L != C_NULL)
        fill_matrix_info_once(mknob_L, L)
    end

    ccall((:ForwardTriangularSolve, LIB_PATH), Void,
           (
            Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void},
           ),
            L.m, L.n, pointer(L.colptr), pointer(L.rowval), pointer(L.nzval),
            pointer(b), pointer(b), fknob)
    b
end

@doc """
Context-sensitive forward triangular solver, equivalent to x = Base.SparseMatrix.
fwdTriSolve!(L, b).
"""
function fwdTriSolve!(
    x     :: Vector,
    L     :: SparseMatrixCSC{Float64, Int32},
    b     :: Vector,
    fknob :: Ptr{Void}
 )
    # ISSUE: if x is not allocated space yet in the caller,can we really copy b to x in this way?
    x[:] = copy(b)
    fwdTriSolve!(L, x, fknob)
    x
end

@doc """
Context-sensitive backward triangular solver.
"""
function bwdTriSolve!(
    U     :: SparseMatrixCSC{Float64, Int32},
    b     :: Vector,
    fknob :: Ptr{Void}
 )
    if fknob != C_NULL
        mknob_U = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
        assert(mknob_U != C_NULL)
        fill_matrix_info_once(mknob_U, U)
    end

    ccall((:BackwardTriangularSolve, LIB_PATH), Void,
              (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}),
               U.m, U.n, pointer(U.colptr), pointer(U.rowval), pointer(U.nzval),
               pointer(b), pointer(b), fknob)
    b
end

@doc """
Context-sensitive backward triangular solver, equivalent to x = Base.SparseMatrix.
bwdTriSolve!(L, b).
"""
function bwdTriSolve!(
    x     :: Vector,
    U     :: SparseMatrixCSC{Float64, Int32},
    b     :: Vector,
    fknob :: Ptr{Void}
 )
    # ISSUE: if x is not allocated space yet in the caller,can we really copy b to x in this way?
    x[:] = copy(b)
    bwdTriSolve!(U, x, fknob)
    x
end

function speye_int32(m::Integer)
  rowval = [Int32(x) for x in [1:m;]]
  colptr = [Int32(x) for x in [rowval; m + 1]]
  nzval = ones(Float64, m)
  return SparseMatrixCSC(m, m, colptr, rowval, nzval)
end

function ilu(
    A     :: SparseMatrixCSC{Float64, Int32},
    fknob :: Ptr{Void} = C_NULL
 )
    LU_nzval = zeros(nnz(A))

    ccall((:ILU, LIB_PATH), Void,
          (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}),
          size(A, 1), size(A, 2),
          A.colptr, A.rowval, A.nzval,
          LU_nzval,
          fknob)

    LU = SparseMatrixCSC{Cdouble, Cint}(size(A, 1), size(A, 2), copy(A.colptr), copy(A.rowval), LU_nzval)
    LU = LU'

    L = tril(LU, -1) + speye_int32(size(LU, 1))
    U = triu(LU)

    L, U
end

function ichol(
    A     :: SparseMatrixCSC{Float64, Int32}
 )
    L_nzval = zeros(nnz(A))

    ccall((:IChol, LIB_PATH), Void,
          (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}),
          size(A, 1), size(A, 2),
          A.colptr, A.rowval, A.nzval,
          L_nzval,
          C_NULL)

    L = SparseMatrixCSC{Cdouble, Cint}(size(A, 1), size(A, 2), copy(A.colptr), copy(A.rowval), L_nzval)
    L = L'
    L
end

@doc """
Reorder a sparse matrix A and store the result in new_A. A itself is not changed.
If get_permutation is true, then compute the permutation and inverse permutation
vector P and inverse_P. Otherwise, P and inverse_P are given inputs.
One_based_input/output tells if A and new_A are 1-based.
ISSUE: We hard code a sparse matrix format to be SparseMatrixCSC{Cdouble, Cint},
and a permutation and inverse permutation vector to be Vector{Cint}.
Should make it general in future
"""
function reorder_matrix(
    A                :: SparseMatrixCSC{Float64, Int32},
    new_A            :: SparseMatrixCSC{Float64, Int32},
    P                :: Vector,
    inverse_P        :: Vector
)
    ccall((:ReorderMatrix, LIB_PATH), Void,
              (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cint}, Ptr{Cint}),
               A.m, A.n, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval),
               pointer(new_A.colptr), pointer(new_A.rowval), pointer(new_A.nzval),
               pointer(P), pointer(inverse_P))
end

function reorder_matrix!(
    A                :: SparseMatrixCSC{Float64, Int32},
    P                :: Vector,
    inverse_P        :: Vector,
)
    ccall((:ReorderMatrixInplace, LIB_PATH), Void,
              (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cint}, Ptr{Cint}),
               A.m, A.n, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval),
               pointer(P), pointer(inverse_P))
    A
end

@doc """ Inversely reorder vector V with permutation vector P """
function inverse_reorder_vector!(
    V     :: Vector,
    P     :: Vector
)
    ccall((:reorderVectorWithInversePermInplace, LIB_PATH), Void,
         (Ptr{Cdouble}, Ptr{Cint}, Cint),
         pointer(V), pointer(P), length(V))
    V
end

function set_reordering_decision_maker(
    fknob            :: Ptr{Void},
    perm_restriction = PERM_VECTORS_ARE_INDEPENDENT
 )
    ccall((:SetReorderingDecisionMaker, LIB_PATH), Void,
         (Ptr{Void}, Cint),
         fknob, perm_restriction)
end

function get_reordering_vectors(
    fknob :: Ptr{Void}
 )
    m = Ref{Cint}(0)
    n = Ref{Cint}(0)
    row_perm = ccall((:GetRowReorderingVector, LIB_PATH), Ptr{Cint},
         (Ptr{Void}, Ref{Cint}),
         fknob, m)
    row_inv_perm = ccall((:GetRowInverseReorderingVector, LIB_PATH), Ptr{Cint},
         (Ptr{Void}, Ref{Cint}),
         fknob, m)
    col_perm = ccall((:GetColReorderingVector, LIB_PATH), Ptr{Cint},
         (Ptr{Void}, Ref{Cint}),
         fknob, n)
    col_inv_perm = ccall((:GetColInverseReorderingVector, LIB_PATH), Ptr{Cint},
         (Ptr{Void}, Ref{Cint}),
         fknob, n)
    pointer_to_array(row_perm, m[]), pointer_to_array(row_inv_perm, m[]), pointer_to_array(col_perm, n[]), pointer_to_array(col_inv_perm, n[])
end

@doc """
Return a new representation of the sparse matrix, which is in CSC
format, in CSR format.
"""
function create_CSR(A :: SparseMatrixCSC)
    ccall((:CSR_Create, LIB_PATH), Ptr{Void},
        (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
        A.n, A.m, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval))
end

@doc """ Destroy the CSR representation of the sparse matrix. """
function destroy_CSR(A :: Ptr{Void})
    ccall((:CSR_Destroy, LIB_PATH), Void, (Ptr{Void},), A)
end

@doc """ w = (alpha*A*x + beta*y + gamma).*z """
function SpMV!(
    w     :: Vector,
    alpha :: Number,
    A     :: SparseMatrixCSC,
    x     :: Vector,
    beta  :: Number,
    y     :: Vector,
    gamma :: Number,
    z     :: Vector,
    fknob :: Ptr{Void} = C_NULL # fknob == C_NULL => A is symmetric so we can directly call CSR_MultiplyWithVector with CSC matrix
)
    assert(length(w) == length(y) || length(y) == 0)
    assert(length(w) == length(z) || length(z) == 0)

    if use_SPMP
        if fknob == C_NULL
          A1 = create_CSR(A)
          ccall((:CSR_MultiplyWithVector, LIB_PATH), Void,
                (Ptr{Cdouble}, Cdouble, Ptr{Void}, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
                pointer(w), alpha, A1, pointer(x), beta, length(y) == 0 ? C_NULL : pointer(y), gamma, length(z) == 0 ? C_NULL : pointer(z))
          destroy_CSR(A1)
        else
          mknob_A = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
          assert(mknob_A != C_NULL)
          fill_matrix_info_once(mknob_A, A)

          ccall((:SpMV, LIB_PATH), Void,
                (Cint, Cint,
                 Ptr{Cdouble},
                 Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                 Ptr{Cdouble},
                 Cdouble, Ptr{Cdouble},
                 Cdouble,
                 Ptr{Cdouble},
                 Ptr{Void}),
                size(A, 1), size(A, 2),
                w,
                alpha, A.colptr, A.rowval, A.nzval,
                x,
                beta, length(y) == 0 ? C_NULL : y,
                gamma,
                length(z) == 0 ? C_NULL : z,
                fknob)
        end
        w
    else
        # use Julia implementation
        w[:] = alpha * A * x + beta * y + gamma
    end
end

@doc """ w = alpha*A*x + beta*y + gamma """
SpMV!(
    w     :: Vector,
    alpha :: Number,
    A     :: SparseMatrixCSC,
    x     :: Vector,
    beta  :: Number,
    y     :: Vector,
    gamma :: Number,
    fknob :: Ptr{Void} = C_NULL) =
    SpMV!(w, alpha, A, x, beta, y, gamma, Array{Float64,1}[], fknob)

@doc """ y = A*x """
SpMV!(y :: Vector, A :: SparseMatrixCSC, x :: Vector, fknob :: Ptr{Void} = C_NULL) =
    SpMV!(y, one(eltype(x)), A, x, zero(eltype(y)), y, zero(eltype(x)), fknob)

@doc """ w = alpha*A*x + beta*y """
function SpMV!(
    w     :: Vector{Complex128},
    alpha :: Number,
    A     :: SparseMatrixCSC,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128},
    fknob :: Ptr{Void} = C_NULL # fknob == C_NULL => A is symmetric so we can directly call CSR_MultiplyWithVector with CSC matrix
)
    assert(length(w) == length(y) || length(y) == 0)

    if use_SPMP
        if fknob == C_NULL
          A1 = create_CSR(A)
          ccall((:CSR_MultiplyWithComplexVector, LIB_PATH), Void,
                (Ptr{Complex128}, Cdouble, Ptr{Void}, Ptr{Complex128}, Cdouble, Ptr{Complex128}, Cdouble, Ptr{Complex128}),
                pointer(w), alpha, A1, pointer(x), beta, length(y) == 0 ? C_NULL : pointer(y), 0, C_NULL)
          destroy_CSR(A1)
        else
          mknob_A = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
          assert(mknob_A != C_NULL)
          fill_matrix_info_once(mknob_A, A)

          ccall((:SpMVComplex, LIB_PATH), Void,
                (Cint, Cint,
                 Ptr{Complex128},
                 Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                 Ptr{Complex128},
                 Cdouble, Ptr{Complex128},
                 Cdouble,
                 Ptr{Complex128},
                 Ptr{Void}),
                size(A, 1), size(A, 2),
                w,
                alpha, A.colptr, A.rowval, A.nzval,
                x,
                beta, length(y) == 0 ? C_NULL : y,
                0,
                C_NULL,
                fknob)
        end
    else
        # use Julia implementation
        w[:] = alpha * A * x + beta * y
    end
end

@doc """ w = alpha*A*x + beta*y + gamma """
# TODO: need a real 8 args SpMV! gamma is not handled yet.
SpMV!(
    w     :: Vector{Complex128},
    alpha :: Number,
    A     :: SparseMatrixCSC,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128},
    gamma :: Number,
    fknob :: Ptr{Void} = C_NULL) =
    SpMV!(w, alpha, A, x, beta, y, fknob)

@doc """ alpha*A*x + beta*y + gamma """
function SpMV(
    alpha :: Number,
    A     :: SparseMatrixCSC,
    x     :: Vector,
    beta  :: Number,
    y     :: Vector,
    gamma :: Number,
    fknob :: Ptr{Void} = C_NULL
)
   # When fknob is not null, we just use
   # the CSC input matrix as if it's CSR
   # assuming the input matrix is symmetric
   w = Array(Cdouble, size(A, 1))
   SpMV!(w, alpha, A, x, beta, y, gamma, fknob)
   w
end

@doc """ alpha*A*x + beta*y """
SpMV(alpha :: Number, A :: SparseMatrixCSC, x :: Vector, beta :: Number, y :: Vector, fknob :: Ptr{Void} = C_NULL) =
    SpMV(alpha, A, x, beta, y, zero(eltype(x)), fknob)

@doc """ alpha*A*x + y """
SpMV(alpha :: Number, A :: SparseMatrixCSC, x :: Vector, y :: Vector, fknob :: Ptr{Void} = C_NULL) =
    SpMV(alpha, A, x, one(eltype(y)), y, zero(eltype(x)), fknob)

@doc """ alpha*A*x """
SpMV(alpha :: Number, A :: SparseMatrixCSC, x :: Vector, fknob :: Ptr{Void} = C_NULL) =
    SpMV(alpha, A, x, zero(eltype(x)), Array{Float64,1}[], zero(eltype(x)), fknob)

@doc """ A*x """
SpMV(A :: SparseMatrixCSC, x :: Vector, fknob :: Ptr{Void} = C_NULL) =
    SpMV(one(eltype(x)), A, x, zero(eltype(x)), Array{Float64,1}[], zero(eltype(x)), fknob)

@doc """ Dot product of vector x and y """
function dot(
    x :: Vector,
    y :: Vector
)
  assert(length(x) == length(y))

  if use_SPMP #&& length(x) > 4096
    ccall((:dot, LIB_PATH), Cdouble,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), pointer(x), pointer(y))
  else
    Base.dot(x, y)
  end
end

@doc """ Dot product of vector x and y """
function dot(
    x :: Vector{Complex128},
    y :: Vector{Complex128}
)
  assert(length(x) == length(y))

  if use_SPMP #&& length(x) >= 3072
    ccall((:dot, LIB_PATH), Complex128,
          (Cint, Ptr{Complex128}, Ptr{Complex128}),
          length(x), pointer(x), pointer(y))
  else
    Base.dot(x, y)
  end
end

@doc """ Dot product of vector x and y """
function dot(
    x :: Vector{Complex128},
    y :: Vector{Complex128}
)
  assert(length(x) == length(y))

  if use_SPMP
    ccall((:dot, LIB_PATH), Complex128,
          (Cint, Ptr{Complex128}, Ptr{Complex128}),
          length(x), pointer(x), pointer(y))
  else
    dot(x, y)
  end
end

@doc """ 2-norm of vector x """
function norm(
    x :: Vector
)
  sqrt(dot(x, x))
end

@doc """ 2-norm of vector x """
function norm(
    x :: Vector{Complex128}
)

  if use_SPMP #&& length(x) >= 8192
    ccall((:norm_complex, LIB_PATH), Cdouble,
          (Cint, Ptr{Complex128}),
          length(x), pointer(x))
  else
    Base.norm(x)
  end
end

@doc """ sum of vector x """
function sum(
    x :: Vector
)
  if use_SPMP
    ccall((:sum, LIB_PATH), Cdouble,
          (Cint, Ptr{Cdouble}),
          length(x), x)
  else
    Base.sum(x)
  end
end

function minimum(
    x :: Vector
)
  if use_SPMP
    ccall((:minimum, LIB_PATH), Cdouble,
          (Cint, Ptr{Cdouble}),
          length(x), x)
  else
    Base.minimum(x)
  end
end

function abs!(
    w :: Vector,
    x :: Vector
)
  if use_SPMP
    ccall((:CSR_abs, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), w, x)
    w
  else
    w[:] = Base.abs(x)
  end
end

function abs!(
    w :: Vector{Float64},
    x :: Vector{Complex128}
)
  if use_SPMP
    ccall((:CSR_abs_complex, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Complex128}),
          length(x), w, x)
  else
    w[:] = Base.abs(x)
  end
end

function abs!(
    w :: Vector{Float64},
    x :: Vector{Complex128}
)
  if use_SPMP
    ccall((:CSR_abs_complex, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Complex128}),
          length(x), w, x)
  else
    w[:] = abs(x)
  end
end

function exp!(
    w :: Vector,
    x :: Vector
)
  if use_SPMP
    ccall((:CSR_exp, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), w, x)
    w
  else
    w[:] = Base.exp(x)
  end
end

function log1p!(
    w :: Vector,
    x :: Vector
)
  if use_SPMP
    ccall((:CSR_log1p, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), w, x)
    w
  else
    w[:] = Base.log1p(x)
  end
end

function min!(
    w     :: Vector,
    x     :: Vector,
    alpha :: Number
)
  assert(length(x) == length(w))

  if use_SPMP
    ccall((:min, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble),
          length(x), w, x, alpha)
    w
  else
    w[:] = Base.min(x, alpha)
  end
end

@doc """ copy x to y """
function copy!(
    y :: Vector,
    x :: Vector
)
  WAXPB!(y, 1, x, 0)
end

@doc """ w = x./y """
function element_wise_divide!(
    w :: Vector,
    x :: Vector,
    y :: Vector
)
  assert(length(x) == length(y))

  if use_SPMP
    ccall((:pointwiseDivide, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), pointer(w), pointer(x), pointer(y))
    w
  else
    w[:] = x./y
  end
end

@doc """" Alocate space for and return element-wise division of two vectors. """
function element_wise_divide(
    x :: Vector,
    y :: Vector
)
  w = Array(Cdouble, length(x))
  element_wise_divide!(w, x, y)
  w
end

@doc """" w = x.*y """
function element_wise_multiply!(
    w :: Vector,
    x :: Vector,
    y :: Vector
)
  assert(length(x) == length(y))

  if use_SPMP
    ccall((:pointwiseMultiply, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          length(x), pointer(w), pointer(x), pointer(y))
    w
  else
    w[:] = x.*y
  end
end

@doc """" w = x.*y """
function element_wise_multiply!(
    w :: Vector{Complex128},
    x :: Vector{Float64},
    y :: Vector{Complex128}
)
  assert(length(x) == length(y))

  if use_SPMP
    ccall((:pointwiseMultiplyRealWithComplex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Ptr{Float64}, Ptr{Complex128}),
          length(x), pointer(w), pointer(x), pointer(y))
  else
    w[:] = x.*y
  end
end

@doc """" w = x.*y """
function element_wise_multiply!(
    w :: Vector{Complex128},
    x :: Vector{Float64},
    y :: Vector{Complex128}
)
  assert(length(x) == length(y))

  if use_SPMP
    ccall((:pointwiseMultiplyRealWithComplex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Ptr{Float64}, Ptr{Complex128}),
          length(x), pointer(w), pointer(x), pointer(y))
  else
    w[:] = x.*y
  end
end

@doc """" Alocate space for and return element-wise multiplication of two vectors. """
function element_wise_multiply(
    x :: Vector{Float64},
    y :: Vector{Float64}
)
  w = Array(Cdouble, length(x))
  element_wise_multiply!(w, x, y)
  w
end

@doc """" Alocate space for and return element-wise multiplication of two vectors. """
function element_wise_multiply(
    x :: Vector{Float64},
    y :: Vector{Complex128}
)
  w = Array(Complex128, length(x))
  element_wise_multiply!(w, x, y)
  w
end

@doc """ w = alpha*x + beta*y """
function WAXPBY!(
    w     :: Vector,
    alpha :: Number,
    x     :: Vector,
    beta  :: Number,
    y     :: Vector
)
  assert(length(x) == length(y))
  assert(length(x) == length(w))

  if use_SPMP
    ccall((:waxpby, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
          length(x), pointer(w), alpha, pointer(x), beta, pointer(y))
    w
  else
    w[:] = alpha*x + beta*y
  end
end

@doc """ w = alpha*x + beta*y """
function WAXPBY!(
    w     :: Vector{Complex128},
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128}
)
  assert(length(x) == length(y))
  assert(length(x) == length(w))

  if use_SPMP
    ccall((:waxpby_complex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Complex128, Ptr{Complex128}, Complex128, Ptr{Complex128}),
          length(x), pointer(w), alpha, pointer(x), beta, pointer(y))
  else
    w[:] = alpha*x + beta*y
  end
end

@doc """ w = alpha*x + beta*y """
function WAXPBY!(
    w     :: Vector{Complex128},
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128}
)
  assert(length(x) == length(y))
  assert(length(x) == length(w))

  if use_SPMP
    ccall((:waxpby_complex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Complex128, Ptr{Complex128}, Complex128, Ptr{Complex128}),
          length(x), pointer(w), alpha, pointer(x), beta, pointer(y))
  else
    w[:] = alpha*x + beta*y
  end
end

@doc """ Allocate space for and return alpha*x + beta*y """
function WAXPBY(
    alpha :: Number,
    x     :: Vector{Float64},
    beta  :: Number,
    y     :: Vector{Float64}
)
  w = Array{Cdouble}(length(x))
  WAXPBY!(w, alpha, x, beta, y)
  w
end

@doc """ Allocate space for and return alpha*x + beta*y """
function WAXPBY(
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128}
)
  w = Array(Complex128, length(x))
  WAXPBY!(w, alpha, x, beta, y)
  w
end

@doc """ Allocate space for and return alpha*x + beta*y """
function WAXPBY(
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number,
    y     :: Vector{Complex128}
)
  w = Array(Complex128, length(x))
  WAXPBY!(w, alpha, x, beta, y)
  w
end

@doc """ w = alpha*x + beta """
function WAXPB!(
    w     :: Vector,
    alpha :: Number,
    x     :: Vector,
    beta  :: Number
)
  assert(length(w) == length(x))

  if use_SPMP
    ccall((:waxpb, LIB_PATH), Void,
          (Cint, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Cdouble),
          length(x), pointer(w), alpha, pointer(x), beta)
    w
  else
    w[:] = alpha*x + beta
  end
end

@doc """ w = alpha*x + beta """
function WAXPB!(
    w     :: Vector{Complex128},
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number
)
  assert(length(w) == length(x))

  if use_SPMP
    ccall((:waxpb_complex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Complex128, Ptr{Complex128}, Complex128),
          length(x), pointer(w), alpha, pointer(x), beta)
  else
    w[:] = alpha*x + beta
  end
end

@doc """ w = alpha*x + beta """
function WAXPB!(
    w     :: Vector{Complex128},
    alpha :: Number,
    x     :: Vector{Complex128},
    beta  :: Number
)
  assert(length(w) == length(x))

  if use_SPMP
    ccall((:waxpb_complex, LIB_PATH), Void,
          (Cint, Ptr{Complex128}, Complex128, Ptr{Complex128}, Complex128),
          length(x), pointer(w), alpha, pointer(x), beta)
  else
    w[:] = alpha*x + beta
  end
end

@doc """ Allocate space for and return alpha*x + beta """
function WAXPB(
    alpha :: Number,
    x     :: Vector,
    beta  :: Number
)
  w = Array(Cdouble, length(x))
  WAXPB!(w, alpha, x, beta)
  w
end

@doc """
Compute A * D * B, where A and B have constant structures, D is a diagonal
matrix.

The A and B are inputed as AT and BT (transposed A and B) in Julia CSC format,
which become A and B in CSR.

An additional fknob is passed in, which contains a matrix knob representing the
output of this function, and 3 other matrix knobs representing the input of this
function (A, D, B).

In addition to its original semantics, which is returning A * D *B, this
function also reuses context info (from the input matrix knobs), and generates
new context info (for the output matrix knob).

If the source statements are
    for
        X = A * D * B
they have been transformed into the following form before executing this function:
    C = Symbol("ouptut") # representing the ouptut of the call at this call site
    mknob_C = new_matrix_knob(properties of C)
    mknob_A = new_matrix_knob(properties of A)
    mknob_D = new_matrix_knob(properties of D)
    mknob_B = new_matrix_knob(properties of B)
    mknob_X = new_matrix_knob(properties of X)
    fknob_ADB = new_function_knob()
    add mknob_C, _A, _D, _B to fknob_ADB
    for
        X = ADB(AT, D, BT, fknob_ADB)
        knob_X = mknob_C  #propagate mknob_C information to mknob_X.
                          #This may or may not be one to one copy. For example,
                          # if C is constant valued, but X is only constant
                          # structured, then we may not pass C's memory
                          # to X directly; instead, we might pass a copy of C's
                          # memory.

In general, in context-sensitive optimizations, there are two problems:
(1) how a function generates and reuses context info?
(2) how to propagate context info?
Here we focus on the first problem.

ADB(A', D, B', fknob_ADB) functionality:
    if fknob_ADB == NULL
        return A * D * B

    get the output mknob_C, and input mknob_A, D, B from fknob_ADB
    assert(mknob_C != NULL)
    assert(mknob_A != NULL)
    assert(mknob_D != NULL)
    assert(mknob_B != NULL)

    # Check inputs
    assert(mknob_A->constant_structured)
    assert(mknob_D->diagonal)
    assert(mknob_B->constant_structured)

    # Compute output
    if mknob_C->matrix == NULL
        mknob_C->matrix = adb_inspect(A', B')
    mknob_C->matrix = CSR_ADB(mknob_C->matrix,  A', B', diag(D))

    # Return the result.
    return mknob_C->matrix in CSC format
"""
function ADB(
    A     :: SparseMatrixCSC,
    D     :: SparseMatrixCSC,
    B     :: SparseMatrixCSC,
    fknob :: Ptr{Void}
)
    if fknob == C_NULL
        # This case should never be reached. We leave it here only to show how
        # a library function can have the same interface with and without
        # context info.
        assert(false)
        return AT' * D * BT'
    end

    mknob_C = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
    mknob_A = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 1)
    mknob_D = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 2)
    mknob_B = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 3)
    assert(mknob_C != C_NULL)
    assert(mknob_A != C_NULL)
    assert(mknob_D != C_NULL)
    assert(mknob_B != C_NULL)

    # Check inputs and output
    assert(ccall((:IsConstantStructured, LIB_PATH), Bool, (Ptr{Void},), mknob_A))
    #TODO: assert(mknob_D->diagonal)
    assert(ccall((:IsConstantStructured, LIB_PATH), Bool, (Ptr{Void},), mknob_B))

    fill_matrix_info_once(mknob_A, A)
    fill_matrix_info_once(mknob_D, D)
    fill_matrix_info_once(mknob_B, B)

    m = size(A, 1)
    n = size(B, 2)
    k = size(A, 2)

    C_colptr_ref = Ref{Ptr{Cint}}(C_NULL)
    C_rowval_ref = Ref{Ptr{Cint}}(C_NULL)
    C_nzval_ref = Ref{Ptr{Cdouble}}(C_NULL)

    ccall((:ADB, LIB_PATH), Void,
         (Cint, Cint, Cint,
           Ref{Ptr{Cint}}, Ref{Ptr{Cint}}, Ref{Ptr{Cdouble}},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
           Ptr{Cdouble},
           Ptr{Void}),
         m, n, k,
         C_colptr_ref, C_rowval_ref, C_nzval_ref,
         A.colptr, A.rowval, A.nzval,
         B.colptr, B.rowval, B.nzval,
         D.nzval,
         fknob)

    C_colptr = pointer_to_array(C_colptr_ref[], n + 1)
    nnz = C_colptr[n + 1] - 1
    C_rowval = pointer_to_array(C_rowval_ref[], nnz)
    C_nzval = pointer_to_array(C_nzval_ref[], nnz)
    SparseMatrixCSC{Cdouble, Cint}(m, n, C_colptr, C_rowval, C_nzval)
end

# A*A'
function SpSquareWithEps(
    A     :: SparseMatrixCSC{Float64, Int32},
    eps   :: Number,
    fknob :: Ptr{Void} = C_NULL)

    if fknob != C_NULL
        mknob_A = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
        assert(mknob_A != C_NULL)
        fill_matrix_info_once(mknob_A, A)
    end

    m = size(A, 1)
    n = size(A, 2)

    C_colptr_ref = Ref{Ptr{Cint}}(C_NULL)
    C_rowval_ref = Ref{Ptr{Cint}}(C_NULL)
    C_nzval_ref = Ref{Ptr{Cdouble}}(C_NULL)

    ccall((:SpSquareWithEps, LIB_PATH), Ptr{Void},
          (Cint, Cint,
           Ref{Ptr{Cint}}, Ref{Ptr{Cint}}, Ref{Ptr{Cdouble}},
           Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
           Cdouble,
           Ptr{Void}),
          m, n,
          C_colptr_ref, C_rowval_ref, C_nzval_ref,
          A.colptr, A.rowval, A.nzval,
          eps,
          fknob)

    C_colptr = pointer_to_array(C_colptr_ref[], n + 1)
    nnz = C_colptr[n + 1] - 1
    C_rowval = pointer_to_array(C_rowval_ref[], nnz)
    C_nzval = pointer_to_array(C_nzval_ref[], nnz)
    SparseMatrixCSC{Cdouble, Cint}(m, n, C_colptr, C_rowval, C_nzval)
end

# alpha*A + beta*B
function SpAdd(
    alpha :: Number,
    A     :: SparseMatrixCSC{Float64, Int32},
    beta  :: Number,
    B     :: SparseMatrixCSC{Float64, Int32})

    m = size(A, 1)
    n = size(A, 2)
    assert(m == size(B, 1))
    assert(n == size(B, 2))

    C_colptr_ref = Ref{Ptr{Cint}}(C_NULL)
    C_rowval_ref = Ref{Ptr{Cint}}(C_NULL)
    C_nzval_ref = Ref{Ptr{Cdouble}}(C_NULL)

    ccall((:SpAdd, LIB_PATH), Ptr{Void},
          (Cint, Cint,
           Ref{Ptr{Cint}}, Ref{Ptr{Cint}}, Ref{Ptr{Cdouble}},
           Cdouble,
           Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
           Cdouble,
           Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
           Ptr{Void}),
          m, n,
          C_colptr_ref, C_rowval_ref, C_nzval_ref,
          alpha,
          A.colptr, A.rowval, A.nzval,
          beta,
          B.colptr, B.rowval, B.nzval,
          C_NULL)

    C_colptr = pointer_to_array(C_colptr_ref[], n + 1)
    nnz = C_colptr[n + 1] - 1
    C_rowval = pointer_to_array(C_rowval_ref[], nnz)
    C_nzval = pointer_to_array(C_nzval_ref[], nnz)
    SparseMatrixCSC{Cdouble, Cint}(m, n, C_colptr, C_rowval, C_nzval)
end

function trace(
    A   :: SparseMatrixCSC{Float64, Int32})

    ccall((:Trace, LIB_PATH), Cdouble,
          (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
          size(A, 1), A.colptr, A.rowval, A.nzval)
end

@doc """
Compute the sparse Cholesky factorization of a sparse matrix A, where A is
constant in structure.

Functionality:
    if fknob == NULL
        return cholfact_int32(A)

    get the output/input mknob_O and A from fknob

    # Check input
    assert(mknob_A->constant_structured)
    if mknob_A->matrix == NULL
        return cholfact_int32(A)

    if mknob_O->dss_handle == NULL
        mknob_O->dss_handle = dss_analyze(mknob_A->matrix)
    dss_factor(mknob_O->dss_handle, mknob_A->matrix)
"""
function cholfact_int32(
    A     :: SparseMatrixCSC,
    fknob :: Ptr{Void}
)
    if fknob == C_NULL
        # This case should never happen
        assert(false)
        return cholfact_int32(A)
    end

    mknob_A = ccall((:GetMatrixKnob, LIB_PATH), Ptr{Void}, (Ptr{Void}, Cint), fknob, 0)
    assert(mknob_A != C_NULL)

    fill_matrix_info_once(mknob_A, A)

    ccall((:CholFact, LIB_PATH), Ptr{Void},
         (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Void}),
         size(A, 1), size(A, 2), A.colptr, A.rowval, A.nzval, fknob)
end

@doc """
Compute y = R \ t, where R is constant in structure.

Functionality:
    if fknob == NULL
       return y = R \ t
   get the input mknob_R from fknob
   if mknob_R->dss_handle == NULL
       return y = R \ t2
   else
       opt = MKL_DSS_DEFAULTS
       dss_solve!(mknob_R->dss_handle, t, y)
"""
function cholfact_inverse_divide!(
    y     :: Vector,
    R     :: Any,
    b     :: Vector,
    fknob :: Ptr{Void}
)
    if fknob == C_NULL
        # This case should never happen
        assert(false)
        y[:] = R \ b
        return y
    end

    ccall((:CholFactInverseDivide, LIB_PATH), Void,
          (Ptr{Void}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}),
          R, y, b, fknob)
    y
end

@doc """
Compute y = R \ t, where R is constant in structure.

Functionality:
    if fknob == NULL
       return y = R \ t
   get the input mknob_R from fknob
   if mknob_R->dss_handle == NULL
       return y = R \ t2
   else
       opt = MKL_DSS_DEFAULTS
       dss_solve!(mknob_R->dss_handle, t, y)
"""
function cholfact_inverse_divide(
    R     :: Any,
    b     :: Vector,
    fknob :: Ptr{Void}
)
    y = zeros(length(b))
    cholfact_inverse_divide!(y, R, b, fknob)
    y
end

# Some symbolic names for each permutation vector.
const NOT_PERM_YET = 0 # Maybe permute, but the color is not decided yet.
const ROW_PERM     = 1
const ROW_INV_PERM = 2
const COL_PERM     = 3
const COL_INV_PERM = 4
const NEVER_PERM   = 5 # Never permute. This state can not change to any other.

# The row and column permutation vectors might have some restrictions from the
# source code. For example, from pcg source code: "z=L\r" and "dot(r,z)", we
# know "Lz=r" and "dot(r,z)", and therefore:
#    L's row permutation vector is the same as r's, which is the same as z's,
#    which is inverse to L's column permutation vector.
# Therefore,
#    L's row and column permutation vector must be inverse to each other.
const PERM_VECTORS_ARE_INDEPENDENT = 0
const PERM_VECTORS_ARE_INVERSE     = 1
const PERM_VECTORS_ARE_SAME        = 2

const LOG_REORDERING = false

@doc """
Reorder arrays based on the the given permutation and inverse permutation vectors.
If inverse_reorder is true, do inverse reordering instead.
"""
function reorder_arrays(
    inverse_reorder     :: Bool,
    permutation_vectors :: Vector,
    arrays...
)
    if LOG_REORDERING
      println("reorder_arrays")
    end

    matrices_done = false
    array_tuple   = arrays[1]
    i = 1
    while i <= length(array_tuple)
        A = array_tuple[i]
        if A == :__delimitor__
            # A delimiter that separate matrices from vectors
            # It means that all matrices have been processed
            matrices_done = true
            i             = i + 1
            continue
        end
        if matrices_done
            assert(typeof(A) <: AbstractVector)
            # For each vector, 1 permutation vector follows
            perm_vec_idx = array_tuple[i + 1]
            perm         = permutation_vectors[perm_vec_idx]
            i    = i + 2
            if inverse_reorder
                inverse_reorder_vector!(A, perm)
            else
                # reordering(v, ROW_PERM) is equivalent to
                # reverse_reordering(v, ROW_INV_PERM) but is slower. So change
                # to the second form if possible.
                if perm_vec_idx == ROW_PERM || perm_vec_idx == COL_PERM
                    perm = permutation_vectors[perm_vec_idx + 1]
                    inverse_reorder_vector!(A, perm)
                else
                    assert(perm_vec_idx == ROW_INV_PERM || perm_vec_idx == COL_INV_PERM)
                    perm = permutation_vectors[perm_vec_idx - 1]
                    inverse_reorder_vector!(A, perm)
                end
            end
        else
            assert(typeof(A) <: AbstractSparseMatrix)
            # For each matrix, 2 permutation vectors and 1 matrix knob follow
            perm1 = permutation_vectors[array_tuple[i + 1]]
            perm2 = permutation_vectors[array_tuple[i + 2]]
            mknob = array_tuple[i + 3]
            fill_matrix_info_once(mknob, A)
            i     = i + 4
            if inverse_reorder
                reorder_matrix!(A, perm2, perm1)
            else
                reorder_matrix!(A, perm1, perm2)
            end
        end
    end
end

@doc """
The first time an reordering decision maker (a function, represented by the
fknob) has decided that reordering should be done, this function should be
called to reorder the specified arrays.

Status is a vector of the following elements:
(1) reordering_done:       Bool, true if reordering has already been done
(2) row_perm:              row permutation vector
(3) row_inv_perm:          row inverse permutaiton vector
(4) col_perm:              column permutation vector
(5) col_inv_perm:          column inverse permutaiton vector
(6) reordering_time:       time spent on reordering so far (for statistics)
The vector is updated when the function returns.
"""
function reordering(
    fknob  :: Ptr{Void},
    status :: Vector,
    arrays...
)
    reordering_done = status[1]
    if !reordering_done
        row_perm, row_inv_perm, col_perm, col_inv_perm = get_reordering_vectors(fknob)
        if length(row_perm) > 0
            assert(length(row_inv_perm) > 0 && length(col_perm) > 0 && length(col_inv_perm) > 0)

            status[2] = row_perm
            status[3] = row_inv_perm
            status[4] = col_perm
            status[5] = col_inv_perm

            reordering_time  = status[6]
            reordering_time -= time()

            reorder_arrays(false, status[2 : 5], arrays)

            reordering_time += time()
            status[6]        = reordering_time
            status[1]        = true # reordering is done
        end
    end
end

@doc """
Reverse reorder the specified arrays.

Status is a vector of the following 4 elements. See the above reordering()
function for details.
"""
function reverse_reordering(
    status :: Vector,
    arrays...
)
    reordering_done = status[1]
    if reordering_done
        reordering_time  = status[6]
        reordering_time -= time()

        reorder_arrays(true, status[2 : 5], arrays)

        reordering_time += time()
        status[6]        = reordering_time
    end
end

@doc """
Read a matrix market file. Force it to be symmetric if force_symmetric is true.
Error when symmetric_only is true but the matrix is not symmetric (and not
forced to be symmetric). It calls SpMP library to do the actual work, which is
faster than a pure Julia version.
"""
function matrix_market_read(
    filename        :: AbstractString,
    symmetric_only  :: Bool = false,
    force_symmetric :: Bool = false
)
    # Read header to figure out if it's an array or sparse matrix
    mmfile = open(filename,"r")
    # Read first line
    firstline = chomp(readline(mmfile))
    tokens = split(firstline)
    if length(tokens) != 5
        throw(ParseError(string("Not enough words on first line: ", ll)))
    end
    if tokens[1] != "%%MatrixMarket"
        throw(ParseError(string("Not a valid MatrixMarket header:", ll)))
    end
    (head1, rep, field, symm) = map(lowercase, tokens[2:5])
    if head1 != "matrix"
        throw(ParseError("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)"))
    end
    close(mmfile)

    if rep == "array"
      m_ref = Ref{Cint}(0)
      n_ref = Ref{Cint}(0)
      v_ref = Ref{Ptr{Cdouble}}(C_NULL)

      ccall((:load_vector_matrix_market, LIB_PATH), Void,
           (Ptr{Cchar}, Ref{Ptr{Cdouble}}, Ref{Cint}, Ref{Cint}),
           filename, v_ref, m_ref, n_ref)
      assert(n_ref[] == 1)

      return pointer_to_array(v_ref[], m_ref[])
    else
      # Read into a COO array
      sizes::Vector{Cint} = [0, 0, 0, 0]
      ccall((:load_matrix_market_step1, LIB_PATH), Void,
          (Ptr{UInt8}, Ptr{Cint}, Bool, Bool), filename, pointer(sizes), force_symmetric, false)

      is_symmetric::Cint = sizes[1] # 0/1: true/false
      if symmetric_only && is_symmetric == 0
          # We cannot handle asymmetric matrix right now, since we need a matrix
          # to be symmetric, and thus we can use its CSC representation as CSR (The
          # SpMP library is based on CSR, while Julia is based on CSC.)
          throw(AsymmetricMatrixMarketFile(filename))
      end

      n::Cint = sizes[2]
      m::Cint = sizes[3]
      assert(is_symmetric == 0 || m == n)
      nnz::Cint = sizes[4]  # Note: if symmetric, nnz includes elements in both
                            # upper and lower triangle

      v = Array(Cdouble, nnz)
      i = Array(Cint, nnz)
      j = Array(Cint, n + 1)

      A = SparseMatrixCSC{Cdouble, Cint}(m, n, j, i, v)

      # Convert the COO array to CSC array.
      ccall((:load_matrix_market_step2, LIB_PATH), Void,
          (Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Bool),
          filename, pointer(j), pointer(i), pointer(v), pointer(sizes), true)

      return A
    end
end

function lbfgs_compute_direction!(
    dx,
    k     :: Int,
    it    :: Int,
    n     :: Int,
    S,
    Y,
    dfk)

    if use_SPMP
      ccall((:LBFGSComputeDirection, LIB_PATH), Void,
            (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            k, it, n, S, Y, dfk, dx)
    else
      a = zeros(k)

      mm = min(it - 1, k)
      begin_idx = it + k - mm - 2

      yk = 1
      if it <= k
      else
        skm1 = S[:,(begin_idx + mm)%k + 1]
        ykm1 = Y[:,(begin_idx + mm)%k + 1]
        yk = dot(skm1, ykm1)/dot(ykm1, ykm1)
      end

      p = zeros(mm, 1)
      for i = 1:mm
        p[i] = 1.0/dot(S[:,(begin_idx + i)%k + 1], Y[:,(begin_idx + i)%k + 1])
      end

      q = copy(dfk)
      for i = mm:-1:1
        si = S[:,(begin_idx + i)%k + 1]
        yi = Y[:,(begin_idx + i)%k + 1]
        a[i] = p[i]*dot(si, q)
        q = q - a[i]*yi
      end

      r = yk*q
      for i = 1:mm
        si = S[:,(begin_idx + i)%k + 1]
        yi = Y[:,(begin_idx + i)%k + 1]
        b = p[i]*dot(yi, r)
        r = r + si*(a[i] - b)
      end

      dx[:] = -r
    end
end

function lbfgs_loss_function1(yXw::Vector, w::Vector, lambda::Number)
  m = length(yXw)

  if use_SPMP
    ccall((:LBFGSLossFunction1, LIB_PATH), Cdouble,
          (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble),
          m, length(w), yXw, w, lambda)
  else
    s = sum(log1p(exp(-abs(yXw))) - min(yXw, 0))
    s/m + (lambda/2)*norm(w)^2
  end
end

function lbfgs_loss_function2!(temp::Vector, y::Vector, yXw::Vector)
  if use_SPMP
    ccall((:LBFGSLossFunction2, LIB_PATH), Void,
          (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}),
          temp, length(y), y, yXw)
  else
    temp[:] = y./(1 + exp(yXw))
  end
end

function get_spmp_spmv_time()
  ccall((:GetSpMPSpMVTime, LIB_PATH), Cdouble, ())
end

function reset_spmp_spmv_time()
  ccall((:ResetSpMPSpMVTime, LIB_PATH), Void, ())
end

function get_spmp_trisolve_time()
  ccall((:GetSpMPTriangularSolveTime, LIB_PATH), Cdouble, ())
end

function reset_spmp_trisolve_time()
  ccall((:ResetSpMPTriangularSolveTime, LIB_PATH), Void, ())
end

function get_knob_spmv_time()
  ccall((:GetKnobSpMVTime, LIB_PATH), Cdouble, ())
end

function reset_knob_spmv_time()
  ccall((:ResetKnobSpMVTime, LIB_PATH), Void, ())
end

function get_knob_trisolve_time()
  ccall((:GetKnobTriangularSolveTime, LIB_PATH), Cdouble, ())
end

function reset_knob_trisolve_time()
  ccall((:ResetKnobTriangularSolveTime, LIB_PATH), Void, ())
end

function set_knob_log_level(level)
  ccall((:SetLogLevel, LIB_PATH), Void, (Cint,), level)
end

function set_reuse_inspection(reuse)
  ccall((:SetReuseInspection, LIB_PATH), Void, (Bool,), reuse)
end

function set_collective_structure_prediction(collective)
  ccall((:SetCollectiveStructurePrediction, LIB_PATH), Void, (Bool,), collective)
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using .Glimpse
const Gl = Glimpse
using ColorSchemes
using .Glimpse: RGBf0

struct WfInterface <: Gl.System end

function Gl.update(::WfInterface, m::AbstractLedger)
    wfman = Gl.singleton(m, WfManager)
    curvis = wfman.entities[wfman.current]

    gui_func = () -> begin
        Gl.Gui.SetNextWindowPos((650, 200), Gl.Gui.ImGuiCond_FirstUseEver)
        Gl.Gui.SetNextWindowSize((550, 680), Gl.Gui.ImGuiCond_FirstUseEver)
        Gl.Gui.Begin("Wannier Function Manager")
        Gl.Gui.SetWindowFontScale(3.3f0)
        curid = Int32(wfman.current)
        Gl.Gui.@c Gl.Gui.InputInt("Id", &curid)
        if 0 < curid <= length(wfman.entities)
            wfman.current = Int(curid)
        elseif curid <= 0
            wfman.current = length(wfman.entities)
        else
            wfman.current = 1
        end
        Gl.Gui.@c Gl.Gui.InputDouble("dt", &wfman.dt,0.01,0.01,"%.3f")
        if wfman.dt <= 0.0
            wfman.dt += 0.01
        end
        if !wfman.iterating
            if Gl.Gui.Button("Iterate")
                wfman.iterating = true
            end
        else
            if Gl.Gui.Button("Stop Iterating")
                wfman.iterating = false
            end
        end
        Gl.Gui.End()
    end
    push!(Gl.singleton(m, Gl.GuiFuncs).funcs, gui_func)

    if wfman.iterating
        wfman.current_t += m[Gl.TimingData][1].dtime
        if wfman.current_t >= wfman.dt
            wfman.current += 1
            if wfman.current >= length(wfman.entities)
                wfman.current= 1
            end
            wfman.current_t = 0.0
        end
    end

    vis = m[Gl.Visible]
    for e in wfman.entities
        vis[e] = Gl.Visible(e == curvis)
    end
end

@component @with_kw mutable struct WfManager
    entities::Vector{Entity}
    wfuncs::Vector{<:WannierFunction}
    current::Int = 1
    dt::Float64 = 0.5
    current_t::Float64 = 0.0
    iterating::Bool = false
end

"""
    visualize_wfuncs(wfuncs::Vector{<:WannierFunction}, str::Structure;
                     iso_ratio = 1/4,
                     alpha     = 0.6,
                     material  = Gl.Material(),
                     phase_channel = Up())

Visualize the wannierfunctions in a Diorama.
The isosurface is determined from the charge, where the coloring signifies the phase of the wavefunction at that point in space.
The `iso_ratio` will be used to determine the isosurface values as the ratio w.r.t the maximum value.
The `phase_channel` denotes whether the phase of the spin up or down channel should be shown in the case of spinor Wannierfunctions.
"""
function visualize_wfuncs(wfuncs::Vector{<:WannierFunction}, str::Structure; kwargs...)
    dio = Diorama(background = RGBAf0(60/255, 60/255, 60/255, 1.0f0))
    visualize_wfuncs!(dio, wfuncs, str; kwargs...)
    return dio
end

function visualize_wfuncs!(dio::Diorama, wfuncs::Vector{<:WannierFunction}, str::Structure;
                          iso_ratio = 1/4,
                          alpha = 0.6,
                          material = Gl.Material(),
                          phase_channel = Up())

    DFControl.Display.add_structure!(dio, str)

    phase_id = length(wfuncs[1].values[1]) > 1 ? (phase_channel == Up() ? 1 : 2) : 1
    colfunc(x) = RGBf0(get(ColorSchemes.rainbow, (angle(x[phase_id]) .+ π)/2π))

    wfentities = fill(Entity(0), length(wfuncs))
    grid = Gl.Grid([Point3f0(w...) for w in wfuncs[1].points])

    for (i, w) in enumerate(wfuncs)
        d = Float32.(norm.(w.values))
        geom = Gl.DensityGeometry(d, iso_ratio * maximum(d))
        color = Gl.DensityColor(colfunc.(w.values))

        wfentities[i] =  Entity(dio, Gl.Spatial(), grid, geom, color, material, Gl.Alpha(alpha))

    end

    man = Entity(dio, WfManager(entities=wfentities, wfuncs=wfuncs))
    push!(dio, :simulation, WfInterface())
    update(WfInterface(), dio)
    return dio
end


text:
function build()
    p = pwd()
    cd(Pkg.dir("RNG/deps/Random123"))
    if is_windows()
        run(`mingw32-make`)
    else
        run(`make`)
    end
    cd(p)
end

function have_aesni()
    if VERSION < v"0.5-" || sizeof(Int) != 8
        return false
    end
    @eval begin # use `@eval` to avoid errors while compiling (with julia 0.4)
        ecx = Base.llvmcall(
        """%1 = call { i32, i32, i32, i32 } asm "xchgq  %rbx,\${1:q}\\0A  cpuid\\0A  xchgq  %rbx,\${1:q}",
        "={ax},=r,={cx},={dx},0,~{dirflag},~{fpsr},~{flags}"(i32 1)
        %2 = extractvalue { i32, i32, i32, i32 } %1, 2
        ret i32 %2""", UInt32, Tuple{})
        (ecx >> 25) & 1 == 1
    end
end

check_compiler() = success(`gcc --version`)

if have_aesni() && check_compiler()
    build()
else
    warn("AES-NI will not be compiled.")
end


text:
module GoToTargetDirectedModule

import ..GoToTargetUndirectedModule as GTTUM
import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase

#####
##### game logic
#####

const NUM_OBJECTS = GTTUM.NUM_OBJECTS
const AGENT = GTTUM.AGENT
const WALL = GTTUM.WALL
const TARGET1 = GTTUM.TARGET1
const TARGET2 = GTTUM.TARGET2
const NUM_ACTIONS = 4

mutable struct GoToTargetDirected{R, RNG} <: GW.AbstractGridWorld
    env::GTTUM.GoToTargetUndirected{R, RNG}
    agent_direction::Int
end

function GoToTargetDirected(; R = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    env = GTTUM.GoToTargetUndirected(R = R, height = height, width = width, rng = rng)
    agent_direction = rand(rng, 0:GW.NUM_DIRECTIONS-1)
    env = GoToTargetDirected(env, agent_direction)
    GW.reset!(env)
    return env
end

function GW.reset!(env::GoToTargetDirected)
    GW.reset!(env.env)
    env.agent_direction = rand(env.env.rng, 0:GW.NUM_DIRECTIONS-1)
    return nothing
end

function GW.act!(env::GoToTargetDirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    inner_env = env.env
    tile_map = inner_env.tile_map
    agent_position = inner_env.agent_position
    agent_direction = env.agent_direction

    if action in Base.OneTo(2)
        if action == 1
            new_agent_position = GW.move_forward(agent_position, agent_direction)
        else
            new_agent_position = GW.move_backward(agent_position, agent_direction)
        end

        if !tile_map[WALL, new_agent_position]
            tile_map[AGENT, agent_position] = false
            inner_env.agent_position = new_agent_position
            tile_map[AGENT, new_agent_position] = true
        end
    elseif action == 3
        env.agent_direction = GW.turn_left(agent_direction)
    else
        env.agent_direction = GW.turn_right(agent_direction)
    end

    agent_position = inner_env.agent_position
    if tile_map[TARGET1, agent_position] || tile_map[TARGET2, agent_position]
        inner_env.done = true
        if tile_map[2 + inner_env.target, agent_position]
            inner_env.reward = inner_env.terminal_reward
        else
            inner_env.reward = inner_env.terminal_penalty
        end
    else
        inner_env.done = false
        inner_env.reward = zero(inner_env.reward)
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::GoToTargetDirected) = GW.get_height(env.env)
GW.get_width(env::GoToTargetDirected) = GW.get_width(env.env)

GW.get_action_names(env::GoToTargetDirected) = (:MOVE_FORWARD, :MOVE_BACKWARD, :TURN_LEFT, :TURN_RIGHT)
GW.get_object_names(env::GoToTargetDirected) = GW.get_object_names(env.env)

function GW.get_pretty_tile_map(env::GoToTargetDirected, position::CartesianIndex{2})
    characters = ('☻', '█', '✖', '♦', '→', '↑', '←', '↓', '⋅')

    object = findfirst(@view env.env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    elseif object == AGENT
        return characters[NUM_OBJECTS + 1 + env.agent_direction]
    else
        return characters[object]
    end
end

function GW.get_pretty_sub_tile_map(env::GoToTargetDirected, window_size, position::CartesianIndex{2})
    tile_map = env.env.tile_map
    agent_position = env.env.agent_position
    agent_direction = env.agent_direction

    characters = ('☻', '█', '✖', '♦', '→', '↑', '←', '↓', '⋅')

    sub_tile_map = GW.get_sub_tile_map(tile_map, agent_position, window_size, agent_direction)

    object = findfirst(@view sub_tile_map[:, position])
    if isnothing(object)
        return characters[end]
    elseif object == AGENT
        return '↓'
    else
        return characters[object]
    end
end

function Base.show(io::IO, ::MIME"text/plain", env::GoToTargetDirected)
    characters = ('☻', '█', '✖', '♦', '→', '↑', '←', '↓', '⋅')

    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nsub_tile_map:\n"
    str = str * GW.get_pretty_sub_tile_map(env, GW.get_window_size(env))
    str = str * "\nreward: $(env.env.reward)"
    str = str * "\ndone: $(env.env.done)"
    str = str * "\ntarget = $(env.env.target) ($(characters[2 + env.env.target]))"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::GoToTargetDirected) = ('w', 's', 'a', 'd')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: GoToTargetDirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: GoToTargetDirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: GoToTargetDirected} = (env.env.env.tile_map, env.env.env.target, env.env.agent_direction)

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: GoToTargetDirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: GoToTargetDirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: GoToTargetDirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: GoToTargetDirected} = env.env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: GoToTargetDirected} = env.env.env.done

end # module


text:
module BandwidthBenchmark

# stdlibs
using Base.Threads: @threads, nthreads
using Printf
using Statistics

# packages
using DataFrames
using PrettyTables
using Requires

# LIKWID dummies
LIKWID_register(tag) = nothing
LIKWID_start(tag, nthreads) = nothing
LIKWID_stop(tag, nthreads) = nothing
LIKWID_init() = nothing
LIKWID_close() = nothing

# includes
include("allocate.jl")
include("affinity.jl")
include("benchmarks.jl")
include("bwbench.jl")
include("bwscaling.jl")
export bwbench, bwscaling, flopsscaling

function __init__()
    @require LIKWID="bf22376a-e803-4184-b2ed-56326e3bff83" begin
        LIKWID_register(tag) = LIKWID.Marker.registerregion(tag)
        LIKWID_start(tag, nthreads) = begin
            @threads :static for i in 1:nthreads
                LIKWID.Marker.startregion(tag)
            end
        end
        LIKWID_stop(tag, nthreads) = begin
            @threads :static for i in 1:nthreads
                LIKWID.Marker.stopregion(tag)
            end
        end
        LIKWID_init() = begin
            println("BandwidthBenchmark.jl: LIKWID Marker API is enabled!")
            LIKWID.Marker.init()
        end
        LIKWID_close() = LIKWID.Marker.close()
    end
end

end

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
export FOKet,FOBra

abstract type AbstractVanisher end
abstract type AbstractKet <: AbstractVanisher end
abstract type AbstractBra <: AbstractVanisher end

struct FOKet <: AbstractKet
    vanish :: Array{FO,1}
end

struct FOBra <: AbstractBra
    vanish :: Array{FO,1}
end

"""
    function checkpair(bra::AbstractBra,ket::AbstractKet)

    bra and ket should be a pair, which means that: for arbitary Fermionic Operator(FO)
    if vanish(bra,op) then vanish(ket,dag(op)).

    We should check this because if we call the (function simplify) and moving a FO to the right. It would met its adjoint.
    If its adjoint could be vanished by bra, this met will only produce a coefficient and it is free to exit the moving procedure.

    <0|...(no d and d dag)... d d+ ....|0> = <0|...(no d and d dag)... ....|0>
"""
function checkpair(bra::FOBra,ket::FOKet)
    for i in bra.vanish
        if !(dag(i) in ket.vanish)
            return false
        end
    end
    return true
end

export vanish
vanish(bra::FOBra,fo::FO) = fo in bra.vanish
vanish(ket::FOKet,fo::FO) = fo in ket.vanish
vanish(bra::FOBra,pfo::PFO) = pfo.fo in bra.vanish
vanish(ket::FOKet,pfo::PFO) = pfo.fo in ket.vanish