text:
# test parallel communication types

# import the API to make things easier
import PdePumiInterface: PeerData, ScatterData, allocateBuffer, pushKey!,
                         getPeerIdx, getRemotePeerIdx,
                         pushReceivePart!, exchangeKeys, pushValue!,
                         getDataColons, resetBuffers,
                         sendParallelData, receiveParallelData

"""
  Test the API in serial
"""
function test_ScatterData(mesh::PumiMesh{T}) where {T}

  @testset "ScatterData" begin
    data = ScatterData(T, (mesh.dim,), mesh.comm)

    # test the pushKey! API
    @test getPeerIdx(data, 1) == 0
    @test getRemotePeerIdx(data, 1) == 0

    idx = CartesianIndex(2, 3)
    peernum = 2
    e_local = C_NULL
    e_remote = C_NULL
    pushKey!(data, e_local, e_remote, peernum, idx)

    @test getPeerIdx(data, 1) == 0
    @test getRemotePeerIdx(data, 2) == 0
    @test getPeerIdx(data, 2) == 1
    @test getRemotePeerIdx(data, 2) == 0

    pushReceivePart!(data, 3)
    @test getPeerIdx(data, 1) == 0
    @test getRemotePeerIdx(data, 3) == 1
    @test getPeerIdx(data, 2) == 1
    @test getRemotePeerIdx(data, 2) == 0

    peernum = 4
    pushKey!(data, e_local, e_remote, peernum, idx)
    pushReceivePart!(data, 5)

    @test getPeerIdx(data, 1) == 0
    @test getRemotePeerIdx(data, 3) == 1
    @test getPeerIdx(data, 2) == 1
    @test getRemotePeerIdx(data, 2) == 0

    @test getPeerIdx(data, 4) == 2
    @test getRemotePeerIdx(data, 5) == 2

    pushKey!(data, e_local, e_remote, 4, [CartesianIndex(2, 4), CartesianIndex(3, 5)])
    @test length(data.send[2].local_indices) == 3

    for data_i in data.send
      allocateBuffer(data_i, data.dims)
    end

    # test pushValue!
    pushValue!(data, 2, [2.0, 3.0])
    pushValue!(data, 2, [4.0, 5.0])

    @test data.curridx[2] == 3
    @test maximum(abs.(data.send[2].vals[:, 1] - [2.0, 3.0])) < 1e-13
    @test maximum(abs.(data.send[2].vals[:, 2] - [4.0, 5.0])) < 1e-13

  end

  return nothing
end

"""
  Test `initSendToOwner` is correct by testing against the `mesh.vert_coords`
  array (parallel only)
"""
function test_initSendToOwner(mesh::PumiMesh{T}) where {T}

  data_recv = zeros(T, mesh.coord_numNodes*mesh.dim)
  # create lambda function for receiving data
  function calc_func(data::PdePumiInterface.PeerData)
    PdePumiInterface.receiveVecFunction(data, mesh, data_recv)
  end

  @testset "testing initSendToOwner" begin
    data = PdePumiInterface.initSendToOwner(mesh, mesh.coordshape_ptr, (mesh.dim,))

    shr = apf.getNormalSharing(mesh.m_ptr)
    # test that all entities owned by another process are present
    for (entity, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
      owner = apf.getOwner(shr, entity)
      if apf.isSharedShr(shr, entity) && owner != mesh.myrank && apf.countNodesOn(mesh.coordshape_ptr, dim) != 0
        owner_idx = getPeerIdx(data, owner)
        @test entity in data.send[owner_idx]._entities_local
      end
    end  # end iterator

    # test that the local indices were computed correctly
    coords = Array{Float64}(3)
    for data_i in data.send
      for j=1:length(data_i._entities_local)
        apf.getPoint(mesh.m_ptr, data_i._entities_local[j], 0, coords)

        for j=data_i.colptr[j]:(data_i.colptr[j+1]-1)
          idx = data_i.local_indices[j]
          @test maximum(abs.(coords - mesh.vert_coords[:, idx])) < 1e-13
        end
      end
    end

    # test that sending the data works correctly
    exchangeKeys(data)

    # check that all entities received are owned by this process
    idx = 1
    for data_i in data.recv
      for entity in data_i.entities
        @test apf.isSharedShr(shr, entity)
        @test apf.getOwner(shr, entity) == mesh.myrank
      end
    end

    # send the coordinates using a SumReduction, so the receiver will get the
    # n * coordinates, where n is the number of remotes + self
    reduce_op = PdePumiInterface.AssignReduction{T}()
    sendParallelData(data, mesh.vert_coords, reduce_op)

    # also test coords3DTo1D at same time
    coords_vec = zeros(T, mesh.dim*mesh.coord_numNodes)
    coords3DTo1D(mesh, mesh.vert_coords, coords_vec, reduce_op)

    receiveParallelData(data, calc_func)
    coords = Array{Float64}(3)
    for (entity, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
      owner = apf.getOwner(shr, entity)
      #typ = apf.getType(mesh.m_ptr, entity)

      if apf.isSharedShr(shr, entity)
        ncopies = apf.countCopies(shr, entity)
        nlocals = apf.countAdjacent(mesh.m_ptr, entity, mesh.dim)
        for j=1:mesh.coord_numNodesPerType[dim+1]
          apf.getPoint(mesh.m_ptr, entity, j-1, coords)
          for k=1:mesh.dim
            idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, j-1, k-1)
            if owner == mesh.myrank
              @test abs(data_recv[idx] - ncopies*coords[k]) < 1e-13
              @test abs(coords_vec[idx] - coords[k]) < 1e-13
            elseif owner != mesh.myrank
              idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, j-1, k-1)
              @test coords_vec[idx] == 0
            end  # end if
          end  # end k
        end  # end j
      end  # end if
    end  # end iterator

    # test scattering data
    coords = zeros(Float64, 3)
    coords_vec = zeros(T, mesh.dim*mesh.coord_numNodes)
    coords3DTo1D(mesh, mesh.vert_coords, coords_vec, reduce_op)
    for i=1:length(coords_vec)
      coords_vec[i] += mesh.myrank  # add rank-specific offset
    end
    coords_vec2 = copy(coords_vec)

    PdePumiInterface.sendParallelData_rev(mesh, data, coords_vec)
    calc_func2 = (data::PeerData) -> PdePumiInterface.receiveFromOwner(data, mesh, coords_vec2)
    PdePumiInterface.receiveParallelData_rev(data, calc_func2)
    # check the rank-specific offset
    for (entity, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
      for j=1:mesh.coord_numNodesPerType[dim+1]
        apf.getPoint(mesh.m_ptr, entity, j-1, coords)
        for k=1:mesh.dim
          idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, j-1, k-1)

          if apf.isOwned(mesh.normalshr_ptr, entity)
            @test abs(coords_vec2[idx] - (coords[k] + mesh.myrank)) < 1e-13
          else
            owner = apf.getOwner(mesh.normalshr_ptr, entity)
            @test abs(coords_vec2[idx] - (coords[k] + owner)) < 1e-13
          end
        end  # end k
      end  # end j
    end   # end iterator

  end  # end testset
#  error("stop here")

  return nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 987c5c4d-3257-41cf-b8a9-7b3e814ba61b
begin
	using Pkg
	Pkg.add("Plots")
	using Plots
end

# ╔═╡ 9215d698-38ed-11ec-1e52-45e96b1bc12e
md"""
# Chapter 1
## Section 1

Taylor series in the neighborhood of x = a

$$f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n$$

Near a = 0 and for x = $\pm$h we can write at the third order

$$f(h) = f_0 + hf' + \frac{h^2}{2}f''$$
$$f(-h) = f_0 - hf' + \frac{h^2}{2}f''$$

If we sub

$$f(h) - f(-h) = 2hf'$$
$$\Rightarrow f' = \frac{f(h) - f(-h)}{2h}$$

Centered difference
"""

# ╔═╡ c7732e76-19b9-4aa6-a61e-5b142dae1d3e
function funcdiff(f, x::T, dx::T)::T where T <: Number
	(f(x + dx) - f(x - dx)) / (2 * dx)
end

# ╔═╡ 5137159c-90a8-49a5-9aa3-5356f3137662
begin
	x = collect(0:0.01:7)
	y = funcdiff.(sin, x, 0.01)
	plot(x, y)
end

# ╔═╡ 72affc38-4df9-48af-98b7-0b1b92d9d980
function arraydiff(x::Vector{T}, dx::T)::Vector{T} where T <: Number
	(x[2:end] - x[1:end-1]) / (2 * dx)
end

# ╔═╡ 42d3c309-bb3a-479b-91fd-3d26adda0555
begin
	u = collect(0:0.01:7)
	v = arraydiff(sin.(x), 0.01)
	plot(u[2:end], v)
end

# ╔═╡ cb996ce0-a6f1-4176-b226-b1e4b322fabb
md"""
We can apply this method to other orders of derivations and order of the Taylor series.

Increasing the order of the series will increase the accuracy but it will also increase the calculation time.

The value of h also influences the precision

Exercice: For the same function and interval, evaluate the mean error with the analytic derivation.
"""

# ╔═╡ 7b6412e5-2f4f-409a-b959-6c9c4eb29a9e
begin
	allh = collect(0.01:0.01:1)
	alleps = Float64[]
	for h in allh
		x = 0:h:7
		cosine = cos.(x)
		sinus = funcdiff.(sin, x, h)
		eps = sum(sqrt.(cosine.^2 + sinus.^2)) / length(x)
		append!(alleps, eps)
	end
	plot(allh, alleps)
end

# ╔═╡ Cell order:
# ╠═987c5c4d-3257-41cf-b8a9-7b3e814ba61b
# ╟─9215d698-38ed-11ec-1e52-45e96b1bc12e
# ╠═c7732e76-19b9-4aa6-a61e-5b142dae1d3e
# ╠═5137159c-90a8-49a5-9aa3-5356f3137662
# ╠═72affc38-4df9-48af-98b7-0b1b92d9d980
# ╠═42d3c309-bb3a-479b-91fd-3d26adda0555
# ╟─cb996ce0-a6f1-4176-b226-b1e4b322fabb
# ╠═7b6412e5-2f4f-409a-b959-6c9c4eb29a9e

ext: jl

max_stars_repo_name: zazbone/ComputationalPhysics

max_issues_repo_name: zazbone/ComputationalPhysics

max_forks_repo_name: zazbone/ComputationalPhysics


text:
# arithmetic.jl

function Base.:+(qdt::QDate, y::Year)
    qdinfo = QREF.qref(qdt)
    ny = qdinfo.y + value(y)
    nqdinfo = QREF.rqref_strict(ny, qdinfo.m, false, 1)
    nqdinfo1 = QREF.nextmonth(nqdinfo)
    nqdinfo1.m == qdinfo.m && nqdinfo1.leap == qdinfo.leap && (nqdinfo = nqdinfo1)
    ld = QREF.daysinmonth(nqdinfo)
    return jdn2qdate(nqdinfo.j - 1 + (qdinfo.md <= ld ? qdinfo.md : ld))
end
function Base.:-(qdt::QDate, y::Year)
    qdinfo = QREF.qref(qdt)
    ny = qdinfo.y - value(y)
    nqdinfo = QREF.rqref_strict(ny, qdinfo.m, false, 1)
    nqdinfo1 = QREF.nextmonth(nqdinfo)
    nqdinfo1.m == qdinfo.m && nqdinfo1.leap == qdinfo.leap && (nqdinfo = nqdinfo1)
    ld = QREF.daysinmonth(nqdinfo)
    return jdn2qdate(nqdinfo.j - 1 + (qdinfo.md <= ld ? qdinfo.md : ld))
end

const _DAYS_IN_MONTH_F = 29.530589

function Base.:+(qdt::QDate, m::Month)
    oqdinfo = QREF.qref(qdt)
    nqdinfo = QREF._check_qdinfo(QREF.addmonth(oqdinfo, value(m)))
    ld = QREF.daysinmonth(nqdinfo)
    return jdn2qdate(nqdinfo.j - 1 + (oqdinfo.md <= ld ? oqdinfo.md : ld))
end
function Base.:-(qdt::QDate, m::Month)
    oqdinfo = QREF.qref(qdt)
    nqdinfo = QREF._check_qdinfo(QREF.addmonth(oqdinfo, -value(m)))
    ld = QREF.daysinmonth(nqdinfo)
    return jdn2qdate(nqdinfo.j - 1 + (oqdinfo.md <= ld ? oqdinfo.md : ld))
end

function Base.:+(qdt::QDate, d::Day)
    QDate(UTD(days(qdt) + value(d)))
end
function Base.:-(qdt::QDate, d::Day)
    QDate(UTD(days(qdt) - value(d)))
end


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

DF = load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_1.jld")
EV__NGG = DF["NGG"]

DF_NGG = DataFrame(Exp = first.(EV__NGG),Var =last.(EV__NGG), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_NGG.std = sqrt.(DF_NGG.Var)
DF_NGG.beta_log = log.(DF_NGG.beta)
DF_NGG

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_NGG
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line( alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_ngg <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,15)+ labs(y='Std', x='Expectation')+scale_colour_identity()+
  theme_classic()+ theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15), plot.margin = unit(c(0,0, 0, 0),'pt'))
"

R"p_ngg
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_6_1.pdf', width= 4, height = 4)
plot(p_ngg)
dev.off()"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
module Rcl

include(joinpath(@__DIR__, "..", "gen", "LibRcl.jl"))

# module internal
# include("internal.jl")
# end
#
using .LibRcl: RclError
export RclError
#
# module rcl
# include("base.jl")
# include("node_base.jl")
# end
#
# export rcl
# using .rcl: name, namespace
# export name, namespace

# abstract type Node end
# abstract type Logger end
# abstract type Publisher end
# abstract type Timer end
# abstract type Subscription end

# function create_node(name) end
# function create_publisher(node::Node, msg_type, topic_name) end
# function get_logger(node::Node) end
# function info(logger::Logger, message) end
# function publish(publisher::Publisher, msg) end
# function destroy_node(node::Node) end
# function create_timer(node::Node, callback) end
# function destroy_timer(timer::Timer) end
# function create_subscription(node::Node, msg_type, topic_name, callback) end
# function spin_once(node::Node) end

end # module


text:
# conj!: Overwriting operation
# -
# ~
# conj
# sign
# real
# imag
# !
# promote_array_type: Not exported
# ./
# .\
# .^
# +
# div
# mod
# &
# |
# $
# .+
# .-
# .*
# .%
# .<<
# .>>
# rem
# slicedim
# flipdim
# rotl90
# rotr90
# rot180
# transpose!: Overwriting
# transposeblock!: Not exported
# ctranspose!: Overwriting
# ctransposeblock!: Not exported
# ccopy!: Not exported
# transpose
@primitive transpose(x),dy  transpose(dy)
addtest(transpose,rand(2,2))
# ctranspose
@primitive ctranspose(x),dy ctranspose(dy)
addtest(ctranspose,rand(2,2))
# _cumsum_type: Not exported
# cumsum
# cumsum!
# cumsum_pairwise!
# cummin
# cummax


text:
# merge!(bdblk1::BranchDataBlock{T}, bdblk2::BranchDataBlock{T}) where T =
# 	BranchDataBlock{T}()

# function _init_datablocks!(tdblk::TreeDataBlock{T}, t::Tree) where T
# 	push!(t.datablocks, tdblk)
# 	tdblk = length(t.datablocks)
# 	for (p, q) in PreorderIterator(t)
# 		push!(q.datablocks, NodeDataBlock{T}())
# 		p == q && continue
# 		push!(getbranch(p, q), BranchDataBlock{T}())
# 	end
# end

function _delete_datablock!(p::AbstractNode, q::AbstractNode, ind::Int)
	for link in q.links
		link.to == p && continue
		_delete_datablock!(q, link.to, ind)
		deleteat!(link.branch.datablocks, ind)
		for dblk in link.branch.datablocks[ind:end]
			dblk.ind -= 1
		end
	end

	deleteat!(q.datablocks, ind)
	# for dblk in q.datablocks[ind:end]
	# 	dblk.ind -= 1
	# end

	return nothing
end

"""
	delete_datablock!(tree, datablock)
	delete_datablock!(tree, index)

Delete a `datablock` from `tree`, including all its node data blocks and branch data blocks.

The methods of this function accept either giving the `TreeDataBlock` object or its index in the `datablocks` array of `tree`.
"""
function delete_datablock!(tree::Tree, ind::Int)
	@assert 0 < ind ≤ length(tree.datablocks)

	_delete_datablock!(tree.origin, tree.origin, ind)

	deleteat!(tree.datablocks, ind)

	@inbounds for dblk in tree.datablocks[ind:end]
		dblk.ind -= 1
	end

	return nothing
end

delete_datablock!(tree::Tree, datablock::AbstractTreeDataBlock) =
	delete_datablock!(t, datablock.ind)


text:
using QuasiGeostrophy, LinearAlgebra, Test, FFTW, BenchmarkTools

# boiler plate definitions
Ωxy = S¹(0,2π) × S¹(0,4π)
Nx = 2^8; Ny = 2^8;
fourier_grid = FourierGrid((Nx, Ny), Ωxy)
x, y = fourier_grid.grid
kx, ky = fourier_grid.wavenumbers
f = @. 2 * sin(x) + sin(2*y) + 0im
∂ˣf = @. 2 * cos(x) + 0im * y
∂ʸf = @. 2 * cos(2*y) + 0im * x
Δf = @. - 2 * sin(x) - 4 * sin(2*y)
f̂ = copy(f)
FFTW.set_num_threads(Threads.nthreads())
P = plan_fft(f) # typeof(iP) <: AbstractFFTs.ScaledPlan
iP = plan_ifft(f)
mul!(f̂, P, f)
∂x = FourierOperator(im .* kx)
∂y = FourierOperator(im .* ky)
Δ = ∂x^2 + ∂y^2
∫𝒢dxdy = inv(Δ)

@testset "FourierTests" begin
    tolerance = eps(1.0) * 1e4
    bool = norm(∂ˣf - iP * ∂x(f̂)) ./ norm(∂ˣf ) < tolerance
    @test bool
    bool = norm(∂ʸf - iP * ∂y(f̂)) ./ norm(∂ʸf ) < tolerance
    @test bool
    bool = norm(Δf  - iP * Δ(f̂) ) ./ norm(Δf) < tolerance
    @test bool
    bool = norm(f - iP * ∫𝒢dxdy(Δ(f̂))) ./ norm(f) < tolerance
    @test bool
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
#Simple example of a Dirichlet problem:
#
#           y
#             |
#             |     u = sin(omega x)
#          Ly ----------------------------------
#             |                                |
#             |                                |
#             |                                |
#       u = 0 |        u_xx + u_yy = 0         | u = 0
#             |                                |
#             |                                |
#             |                                |
#           -----------------------------------|-----  x
#                          u = 0              Lx
#
#Here, omega = 2 n pi / Lx, and the exact solution is
#
#      u = sin(omega x) * sinh(omega y) / sinh(omega Ly)
#

using FinElt
using FinElt.PlanarPoisson
using Printf

include("params.jl")

const omega = 2 * pi / Lx

function exact_u(x)
    return sin(omega*x[1]) * sinh(omega*x[2]) / sinh(omega*Ly)
end

essential_bc = [ "Top", "Bottom", "Left", "Right" ]

maxerr = zeros(refinements+1)
@printf("%10s  %12s  %8s  %8s\n\n",
        "N", "max error", "rate", "seconds")
for k = 0:refinements
    start = time()
    mesh = read_msh_file("rect$k.msh")
    vp = VariationalProblem(mesh, essential_bc)
    assign_bdry_vals!(vp, "Top", exact_u)
    add_bilin_form!(vp, "Omega", grad_dot_grad!)
    A, b = assembled_linear_system(vp)
    ufree = A \ b
    uh = complete_soln(ufree, vp)
    u = get_nodal_vals(exact_u, mesh)
    finish = time()
    maxerr[k+1] = maximum(abs.(uh-u))
    N = length(ufree)
    if k == 0
        @printf("%10d  %12.4e\n", N, maxerr[k+1])
    else
        rate = log2(maxerr[k]/maxerr[k+1])
        elapsed = finish - start
        @printf("%10d  %12.4e  %8.4f  %8.4f\n",
                N, maxerr[k+1], rate, elapsed)
    end
end


text:
using Test, BenchmarkTools

function num_allocs(model)
    dt = 0.1
    x, u = rand(model)
    n,m = size(model)
    z = KnotPoint(x, u, dt)
    ∇c  = zeros(n,n+m)
    ∇cd = zeros(n,n+m+1)
    dynamics(model, x, u)
    jacobian!(∇c, model, z)
    discrete_dynamics(RK3, model, x, u, z.t, dt)
    discrete_jacobian!(RK3, ∇cd, model, z)
    allocs  = @allocated dynamics(model, x, u)
    allocs += @allocated jacobian!(∇c, model, z)
    allocs += @allocated discrete_dynamics(RK3, model, x, u, z.t, dt)
    allocs += @allocated discrete_jacobian!(RK3, ∇cd, model, z)
end

# Double Integrator
dim = 3
di = Dynamics.DoubleIntegrator(dim)
n,m = size(di)
@test (n,m) == (6,3)
x,u = rand(di)
@test num_allocs(di) == 0

# Pendulum
pend = Dynamics.Pendulum()
@test size(pend) == (2,1)
@test num_allocs(pend) == 0

# Car
car = Dynamics.DubinsCar()
@test size(car) == (3,2)
@test num_allocs(car) == 0

# Cartpole
cartpole = Dynamics.Cartpole()
@test size(cartpole) == (4,1)
@test num_allocs(cartpole) == 0

# Quadrotor
quad = Dynamics.Quadrotor2()
@test size(quad) == (13,4)
@test num_allocs(cartpole) == 0

# Infeasible
model = Dynamics.DubinsCar()
inf = InfeasibleModel(model)
@test inf._u == 1:2
@test inf._ui == 3:5
@test size(inf) == (3,5)
x = rand(inf)[1]
u0 = @SVector rand(2)
ui = @SVector rand(3)
u = [u0; ui]
dt = 0.1
@test length(rand(inf)[2]) == 5
@test_throws ErrorException dynamics(inf, x, u)
@test discrete_dynamics(RK3, inf, KnotPoint(x, u, dt)) ==
    (discrete_dynamics(RK3, model, x, u0, 0.0, dt) + ui)

function inf_allocs(inf)
    x,u = rand(inf)
    n,m = size(inf)
    dt = 0.1
    z = KnotPoint(x,u,0.1)
    # allocs = @allocated discrete_dynamics(RK3, inf, x, u, dt)
    ∇c = zeros(n,n+m+1)
    allocs = @allocated discrete_dynamics(RK3, inf, z)
    allocs += @allocated discrete_jacobian!(RK3, ∇c, inf, z)
end
@test inf_allocs(inf) == 0

# Test other functions
car = Dynamics.DubinsCar()
n,m = size(car)
@test zeros(car) == (zeros(n), zeros(m))
@test zeros(Int,car)[1] isa SVector{n,Int}
@test fill(car,0.1) == (fill(0.1,n), fill(0.1,m))
@test ones(Float32,car)[2] isa SVector{m,Float32}

# Test default integrator
x,u = rand(car)
z = KnotPoint(x,u,dt)
@test discrete_dynamics(car, z) == discrete_dynamics(RK3, car, z)


text:
export tensorcontract, TensorNetwork, contract_label!, PlotMeta, contract, ContractionTree
export contract_tree

function _align_eltypes(xs::AbstractArray...)
    T = promote_type(eltype.(xs)...)
    return map(x->eltype(x)==T ? x : T.(x), xs)
end

function _align_eltypes(xs::AbstractArray{T}...) where T
    xs
end

# batched routines
@inline _indexpos(iAs, i)::Int = findfirst(==(i), iAs)
_isunique(x) = length(unique(x)) == length(x)

# can be used in either static or dynamic invoke
@noinline function analyse_dim(iAs, iBs, iOuts)
    # check indices
    @assert _isunique(iAs) "indices in A matrix is not unique $iAs"
    @assert _isunique(iBs) "indices in B matrix is not unique $iBs"
    @assert _isunique(iOuts) "indices in C matrix is not unique $iOuts"
    allinds = [iAs..., iBs..., iOuts...]
    @assert all(i -> count(==(i), allinds) == 2, allinds) "indices does not appear in pairs! got $iAs, $iBs and $iOuts"

    mdims = setdiff(iAs, iBs)
    ndims = setdiff(iBs, iAs)
    kdims = iAs ∩ iBs
    iABs = mdims ∪ ndims
    lA = mdims ∪ kdims
    lB = kdims ∪ ndims
    lOut = mdims ∪ ndims
    pA = indexin(lA, iAs)
    pB = indexin(lB, iBs)
    pOut = indexin(iOuts, lOut)
    pA, pB, pOut, length(kdims)
end

@noinline function analyse_size(pA, sA, pB, sB, nc)
    nA = length(sA)
    nB = length(sB)
    sA1 = Int[sA[pA[i]] for i=1:nA-nc]
    sB2 = Int[sB[pB[i]] for i=nc+1:nB]
    K = mapreduce(i->sB[pB[i]], *, 1:nc, init=1)
    return prod(sA1), K, prod(sB2), [sA1..., sB2...]
end

function tensorcontract(iAs, A::AbstractArray, iBs, B::AbstractArray, iOuts)
    A, B = _align_eltypes(A, B)
    pA, pB, pOut, nc = analyse_dim([iAs...], [iBs...], [iOuts...])
    M, K, N, sOut = analyse_size(pA, size(A), pB, size(B), nc)

    Apr = reshape(_conditioned_permutedims(A, pA), M, K)
    Bpr = reshape(_conditioned_permutedims(B, pB), K, N)
    AB = Apr * Bpr
    AB = _conditioned_permutedims(reshape(AB, sOut...), (pOut...,))
end

function _conditioned_permutedims(A::AbstractArray{T,N}, perm) where {T,N}
    if any(i-> (@inbounds perm[i]!=i), 1:N)
        return permutedims(A, (perm...,))
    else
        return A
    end
end

"""
    TensorNetwork{T,TT<:AbstractTensor{T}}
    TensorNetworks(tensors)

A tensor network, its only field is `tensors` - a vector of `AbstractTensor` (e.g. `LabeledTensor`).
"""
struct TensorNetwork{T,TT<:AbstractTensor{T}}
    tensors::Vector{TT}
end

function TensorNetwork(tensors)
    T = promote_type([eltype(x.array) for x in tensors]...)
    TensorNetwork(AbstractTensor{T}[tensors...])
end

Base.copy(tn::TensorNetwork) = TensorNetwork(copy(tn.tensors))
Base.length(tn::TensorNetwork) = length(tn.tensors)

struct ContractionTree
    left
    right
    function ContractionTree(left::Union{Integer, ContractionTree}, right::Union{Integer, ContractionTree})
        new(left, right)
    end
end
function Base.getindex(ct::ContractionTree, i::Int)
     Base.@boundscheck i==1 || i==2
    i==1 ? ct.left : ct.right
end

function contract_tree(tn::TensorNetwork{T}, ctree; normalizer) where T
    if ctree isa Integer
        t = tn.tensors[ctree]
        if normalizer !== nothing
            normt = normalizer(t)
            return normt, rmul!(copy(t), T(1/normt))
        else
            return one(T), t
        end
    else
        normt1, t1 = contract_tree(tn, ctree[1]; normalizer=normalizer)
        normt2, t2 = contract_tree(tn, ctree[2]; normalizer=normalizer)
        normt = normt1 * normt2
        t = t1 * t2
        if normalizer !== nothing
            normt_ = normalizer(t)
            normt = normt * normt_
            rmul!(t, T(1/normt_))
        end
        return normt, t
    end
end

function contract(tn::TensorNetwork{T}, ctree; normalizer=nothing) where T
    normt, t = contract_tree(tn, ctree; normalizer=normalizer)
    if normt != one(T)
        rmul!(t.array, normt)
    end
    return t
end

function contract_label!(tn::TensorNetwork{T}, label) where {T}
    ts = findall(x->label ∈ x.labels, tn.tensors)
    @assert length(ts) == 2 "number of tensors with the same label $label is not 2, find $ts"
    t1, t2 = tn.tensors[ts[1]], tn.tensors[ts[2]]
    tout = t1 * t2
    deleteat!(tn.tensors, ts)
    push!(tn.tensors, tout)
    return tout, t1.labels ∩ t2.labels
end

using Base.Cartesian
c2l(size::NTuple{0, Int}, c::NTuple{0,Int}) = 1
@generated function c2l(size::NTuple{N, Int}, c::NTuple{N,Int}) where N
    quote
        res = c[1]
        stride = size[1]
        @nexprs $(N-1) i->begin
            res += (c[i+1]-1) * stride
            stride *= size[i+1]
        end
        return res
    end
end

l2c(size::NTuple{0, Int}, l::Int) = ()
@generated function l2c(size::NTuple{N, Int}, l::Int) where N
    quote
        l -= 1
        @nexprs $(N-1) i->begin
            @inbounds l, s_i = divrem(l, size[i])
        end
        $(Expr(:tuple, [:($(Symbol(:s_, i))+1) for i=1:N-1]..., :(l+1)))
    end
end

function Base.show(io::IO, tn::TensorNetwork)
    print(io, "$(typeof(tn).name):\n  $(join(["$(t.meta === nothing ? i : dispmeta(t.meta)) => $t" for (i, t) in enumerate(tn.tensors)], "\n  "))")
end

function Base.show(io::IO, ::MIME"plain/text", tn::TensorNetwork)
    Base.show(io, tn)
end

function adjacency_matrix(tn::TensorNetwork)
    indx = Int[]
    indy = Int[]
    data = Int[]
    for i=1:length(tn)
        for j=1:length(tn)
            if i!=j && any(l->l ∈ tn.tensors[i].labels, tn.tensors[j].labels)
                push!(indx, i)
                push!(indy, j)
                push!(data, 1)
            end
        end
    end
    return sparse(indx, indy, data)
end


text:
module Theano

export tensor, grad, TheanoVariable, TheanoFunction, theanofunction

import Base: (.*), (./), (.^), (+), (-), (*), (/), (^), log, exp, transpose, sum, log1p, mean, size, getindex, ndims

using PyCall
@pyimport theano
@pyimport theano.tensor as T

# TODO we should tighten this up with parameteric types for Theano variables (dimensions, eltype, maybe broadcastability)
typealias TheanoVariable{N} PyObject
typealias TheanoFunction    PyObject

(+)(x::TheanoVariable, y::TheanoVariable) = x[:__add__](y)
(-)(x::TheanoVariable, y::TheanoVariable) = x[:__sub__](y)
(.*)(x::TheanoVariable, y::TheanoVariable) = x[:__mul__](y)
(./)(x::TheanoVariable, y::TheanoVariable) = x[:__div__](y)
(.^){N<:Integer}(x::TheanoVariable, n::N)  = x[:__pow__](n)

(*){N<:Number}(x::TheanoVariable, y::N) = x[:__mul__](y)
(*){N<:Number}(x::N, y::TheanoVariable) = y[:__rmul__](x)

for fn in [:exp, :log, :log1p, :mean, :sum]
  @eval $fn(x::TheanoVariable) = T.$fn(x)
end

transpose(x::TheanoVariable) = ndims(x) = 1 ? T.shape_padleft(x) : T.transpose(x)

dot(x::TheanoVariable, y::TheanoVariable) = T.dot(x, y)
#tensordot(x::TheanoVariable, y::TheanoVariable, axes::Array{(Int,Int)}) = T.tensordot(x, y, axes)

sum(x::TheanoVariable, region::Int) = T.sum(x, region-1)
getindex(x::TheanoVariable, ix::Int) = x[:__getitem__](ix)
#getindex(x::TheanoVariable, ix::Tuple{}) = x[:__getitem__](())

grad(cost::TheanoVariable{1}, wrt::TheanoVariable...) = T.grad(cost, wrt)

theanofunction(inputs::Vector{TheanoVariable}, output::TheanoVariable) = theano.pymember(:function)(inputs, output, allow_input_downcast=true)

ndims(x :: TheanoVariable) = x[:ndim]

size(x::TheanoVariable) = T.shape(x)
# Follow the convention from Mocha.blobs that size(x, -1) gives size(x)[end], etc.
function size(x :: TheanoVariable, dim :: Int)
  N = ndims(x)
  if dim < 0
    dim = N+1 + dim
  end
  return T.shape(x)[dim-1]
end

DTYPES = [Float32 => "float32",
          Float64 => "float64"]

tensor(eltype::DataType, ndims::Int) = T.TensorType(DTYPES[eltype], repmat([false],ndims))()

end

max_stars_repo_name: benmoran/MochaTheano

max_issues_repo_name: benmoran/MochaTheano

max_forks_repo_name: benmoran/MochaTheano


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

import Test

Test.@testset "Package MGVI" begin

include("test_mgvi_impl.jl")
include("test_jacobians.jl")
include("information/test_information.jl")
include("test_samplers.jl")
include("test_normal_mvnormal.jl")

end # testset


text:
"""
    generate_default_project(path_to_project_dir::String)

    Generates a default project structure which contains an empty model file and Defaults.toml file.

    Inputs:
    path_to_project_dir = path to where you want model code to be generated
"""
function generate_default_project(path_to_project_dir::String)

    # before we get too far along, we need to check if the user already has code in the location that we want to generate at -
    # if they do, then move it -
    if (isdir(path_to_project_dir) == true)

        # ok, looks like we may have a conflict - mv the offending code
        if (move_existing_project_at_path(path_to_project_dir) == false)

            # Something happend ... the world is ending ...
            throw(ArgumentError("automatic directory conflict resolution failed. Unable to move existing directory $(path_to_project_dir)"))
        end
    else

        # ok, no code at our target directory -
        mkpath(path_to_project_dir)
    end

    # ok, if we get here, then we have a clean place to generate the default project structure -
    # We need to two things for a project, the defaults file, and a blank network file with all the sections -

    # Transfer distrubtion files to the output -> these files are shared between model types
    transfer_distribution_files("$(path_to_package)/distribution/julia/default_project_files", "$(path_to_project_dir)",".vff")
    transfer_distribution_files("$(path_to_package)/distribution/julia/default_project_files", "$(path_to_project_dir)",".toml")
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
using Test
using POMDPs
using Random

let
    struct M <: POMDP{Int, Int, Char} end
    @test_throws MethodError generate_s(M(), 1, 1, MersenneTwister(4))
    POMDPs.transition(::M, ::Int, ::Int) = [1]
    @test generate_s(M(), 1, 1, MersenneTwister(4)) == 1
    @test_throws MethodError generate_sor(M(), 1, 1, MersenneTwister(4))
    @test_throws MethodError generate_sr(M(), 1, 1, MersenneTwister(4))
    POMDPs.reward(::M, ::Int, ::Int, ::Int) = 0.0
    @test generate_sr(M(), 1, 1, MersenneTwister(4)) == (1, 0.0)
    @test_throws MethodError generate_sor(M(), 1, 1, MersenneTwister(4))
    POMDPs.generate_o(::M, ::Int, ::Int, ::Int, ::AbstractRNG) = `a`
    @test generate_sor(M(), 1, 1, MersenneTwister(4)) == (1, `a`, 0.0)
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>1.0, :alg=>TMJets(abs_tol=1e-15, orderT=10, orderQ=2, max_steps=50_000)),
    "X0"   => Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0]),
    "info" => Dict(
        "name" => "lorenz",
        "dim" => 3,
        "linear" => false,
        "hybrid" => false
    )
])


text:
@testset "400.nth-digit.jl" begin
    @test find_nth_digit(3) == 3
    @test find_nth_digit(11) == 0
    for (i, j) in
        zip([500, 501, 502, 503, 504, 505, 506, 507, 508, 513, 514, 515, 516, 517, 518, 519,
             520], [0, 3, 2, 0, 4, 2, 0, 5, 2, 7, 2, 0, 8, 2, 0, 9, 2])
        @test find_nth_digit(i) == j
    end
end


text:
#
# Run this file to generate Figure 7 in
#
#   D. Neumann, T. Wiese, and W. Utschick, Learning the MMSE Channel Estimator,
#   IEEE Transactions on Signal Processing, 2018.
#
# (It may be necessary to split the loop over the different SNR values into
# several smaller loops as the TensorFlow graph is not reset after each
# loop and memory usage may build up.)
#
push!(LOAD_PATH,".")
using DataFrames
using CSV
import SCM3GPP; const scm = SCM3GPP
import CondNormalMMSE; const mmse = CondNormalMMSE
import CondNormalTF; const cntf = CondNormalTF

include("sim_helpers.jl")
include("more_estimators.jl") # OMP and Genie MMSE

verbose = true
#-------------------------------------
# Simulation parameters
#
write_file = true
filename   = "results/figure7.csv"
nBatches   = 100
nBatchSize = 100

#-------------------------------------
# Channel Model
#
snrs       = 15 #-15:5:15 # [dB]
nAntennas  = 64
AS         = 2.0 # standard deviation of Laplacian (angular spread)
nCoherence = 1
Channel    = scm.SCMMulti(pathAS=AS, nPaths=3)
# method that generates "nBatches" channel realizations
get_channel(nAntennas, nCoherence, nBatches) = scm.generate_channel(Channel, nAntennas, nCoherence=nCoherence, nBatches = nBatches)
# method that samples C_delta from delta prior
get_cov(nAntennas) = scm.toeplitzHe( scm.generate_channel(Channel, nAntennas, nCoherence=1)[2][:] )
# get circulant vector that generates all covariance matrices for arbitrary delta (here: use delta=0)
get_circ_cov_generator(nAntennas) = real(scm.best_circulant_approximation(scm.scm_channel([0.0],[1.0],nAntennas,AS=AS)[2]))

#-------------------------------------
# Learning Algorithm parameters
#
nLayers = 2
nLearningBatchSize = 20
nLearningBatches   = 6000

init_params = Dict{Symbol,Any}()

init_params[:nAntennas]   = [   8;  16;  32;  64]
init_params[:nBatches]    = [1000;1000;1000]
init_params[:nBatchSize]  = [  20;  20;  20]
init_params[:get_channel] = (nAntennas,nBatches) -> get_channel(nAntennas, nCoherence, nBatches)

init_params[:learning_rates] = Dict{Symbol,Vector{Float64}}()
init_params[:learning_rates][:ToepReLU] = 1e-4*64./init_params[:nAntennas]
init_params[:learning_rates][:CircReLU] = 1e-4*64./init_params[:nAntennas]
init_params[:learning_rates][:ToepSoftmax] = 1e-3*ones(init_params[:nAntennas])
init_params[:learning_rates][:CircSoftmax] = 1e-3*ones(init_params[:nAntennas])

results = DataFrame()
# read results from previous run
if isfile(filename)
    results = CSV.read(filename)
end

srand(size(results,1)) # use number of entries in results as seed for random number generator
for iSNR in 1:length(snrs)
    snr = snrs[iSNR]
    init_params[:snrs] = snr*ones(length(init_params[:nBatches]))

    verbose && println("Simulating with ", snr, " [dB] SNR")

    # Conditionally normal estimators
    cn_est = Dict{Symbol,Any}()
    cn_est[:FastMMSE] = mmse.FastMMSE(snr, get_circ_cov_generator(nAntennas))
    cn_est[:CircMMSE] = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = circ_trans)
    cn_est[:ToepMMSE] = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = toep_trans)
    cn_est[:CircML]   = mmse.MLEst(snr, transform = circ_trans)

    # Neural network estimators
    nn_est = Dict{Symbol,Any}()
    nn_est[:ToepReLU]    = cntf.ConvNN(nLayers, init_params[:nAntennas][1],
                                       transform = toep_trans,
                                       learning_rate = init_params[:learning_rates][:ToepReLU][1])
    nn_est[:CircReLU]    = cntf.ConvNN(nLayers, init_params[:nAntennas][1],
                                       transform = circ_trans,
                                       learning_rate = init_params[:learning_rates][:CircReLU][1])
    nn_est[:ToepSoftmax] = cntf.ConvNN(nLayers, init_params[:nAntennas][1],
                                       transform = toep_trans,
                                       learning_rate = init_params[:learning_rates][:ToepSoftmax][1],
                                       activation = cntf.nn.softmax)
    nn_est[:CircSoftmax] = cntf.ConvNN(nLayers, init_params[:nAntennas][1],
                                       transform = circ_trans,
                                       learning_rate = init_params[:learning_rates][:CircSoftmax][1],
                                       activation = cntf.nn.softmax)

    init_hier!(nn_est, init_params, verbose = verbose)
    train!(nn_est, snr = snr, nBatches = nLearningBatches, get_channel = () -> get_channel(nAntennas, nCoherence, nLearningBatchSize), verbose = verbose)

    algs = Dict{Symbol,Any}()
    algs[:GenieMMSE] = (y,h,h_cov) -> mmse_genie(y, h_cov, snr)
    algs[:GenieOMP]  = (y,h,h_cov) -> omp_genie(y, h)
    for (alg,cn) in cn_est
        algs[alg] = (y,h,h_cov) -> mmse.estimate(cn, y)
    end
    for (alg,nn) in nn_est
        algs[alg] = (y,h,h_cov) -> cntf.estimate(nn, y)
    end

    (errs,rates) = evaluate(algs, snr = snr, nBatches = nBatches, get_channel = () -> get_channel(nAntennas, nCoherence, nBatchSize), verbose = verbose)

    for alg in keys(algs)
        new_row = DataFrame(MSE        = errs[alg],
                            rate       = rates[alg],
                            Algorithm  = String(alg),
                            SNR        = snr,
                            nAntennas  = nAntennas,
                            nCoherence = nCoherence)

        if isempty(results)
            results = new_row
        else
            results = vcat(results,new_row)
        end
    end

    if write_file
        CSV.write(filename, results)
    end
    @show results
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
abstract type Meta end

adjlist(m::T) where {T<:Meta} = m.adjlist
update_edge(m::T; kwargs...) where {T<:Meta} = ifelse(haskey(kwargs, :E), kwargs[:E], nothing)
update_vertex(m::T; kwargs...) where {T<:Meta} = ifelse(haskey(kwargs, :X), kwargs[:X], nothing)
update_global(m::T; kwargs...) where {T<:Meta} = ifelse(haskey(kwargs, :u), kwargs[:u], nothing)

aggregate_neighbors(m::T, aggr::Symbol; kwargs...) where {T<:Meta} = error("not implement")
aggregate_edges(m::T, aggr::Symbol; kwargs...) where {T<:Meta} = error("not implement")
aggregate_vertices(m::T, aggr::Symbol; kwargs...) where {T<:Meta} = error("not implement")

all_vertices_data(x) = ifelse(haskey(x, :X), (X=x[:X],), NamedTuple())
all_edges_data(x) = ifelse(haskey(x, :E), (E=x[:E],), NamedTuple())

function adjacent_vertices_data(x, i::Integer, ne::AbstractArray{<:Integer})
    if haskey(x, :X)
        return (x_i=view(x[:X],:,i), x_j=view(x[:X],:,ne))
    else
        return NamedTuple()
    end
end

function incident_edges_data(x, i::Integer, ne::AbstractArray{<:Integer})
    if haskey(x, :E)
        return (e_ij=view(x[:E],:,i,ne), )
    else
        return NamedTuple()
    end
end

function propagate(meta::T; adjl=adjlist(meta), kwargs...) where {T<:Meta}
    gi = GraphInfo(adjl)

    newE = update_edge(meta; gi=gi, kwargs...)

    if haskey(kwargs, :naggr)
        Ē = aggregate_neighbors(meta, kwargs[:naggr]; E=newE, cluster=generate_cluster(newE, gi))
        kwargs = (kwargs..., Ē=Ē)
    end

    newV = update_vertex(meta; kwargs...)

    if haskey(kwargs, :eaggr)
        Ē = aggregate_edges(meta, kwargs[:eaggr]; kwargs...)
        kwargs = (kwargs..., Ē=Ē)
    end
    if haskey(kwargs, :vaggr)
        V̄ = aggregate_vertices(meta, kwargs[:vaggr]; kwargs...)
        kwargs = (kwargs..., V̄=V̄)
    end

    new_u = update_global(meta; kwargs...)

    (newE, newV, new_u)
end

function generate_cluster(M::AbstractMatrix, gi::GraphInfo)
    cluster = similar(M, Int, gi.E)
    @inbounds for i = 1:gi.V
        j = gi.edge_idx[i]
        k = gi.edge_idx[i+1]
        cluster[j+1:k] .= i
    end
    cluster
end


text:
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MIN_SENSE
        model.inner.objective_sense = :Min
    elseif sense == MOI.MAX_SENSE
        model.inner.objective_sense = :Max
    elseif sense == MOI.FEASIBILITY_SENSE
        model.inner.objective_sense = :Feasibility
    else
        error("Unsupported objective sense: $sense")
    end
    return
end

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SAF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SQF}) = true

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction{F}, obj::F) where {F<:Union{SV, SAF, SQF}}
    model.inner.objective_expr = to_expr(obj)
    return
end


text:
# ==========================================================================================================
# Backend funcion

formula(trm::TableRegressionModel) = trm.mf.f

fe_intercept(f::FormulaTerm) = fe_intercept(f.rhs)
fe_intercept(term::StatsModels.TupleTerm) = map(fe_intercept, term)
fe_intercept(term::FunctionTerm) = first(term.exorig.args) == :fe
fe_intercept(term) = false

width(term::FunctionTerm) = first(term.exorig.args) == :fe ? 0 : 1

# Variable dispersion
dof(model::FixedEffectModel) = model.nobs - model.dof_residual + 1


text:
import Random

function logsumexp(arr::AbstractArray{T}) where {T <: Real}
    max_arr = maximum(arr)
    max_arr == -Inf ? -Inf : max_arr + log(sum(exp.(arr .- max_arr)))
end

function logsumexp(x1::Real, x2::Real)
    m = max(x1, x2)
    m == -Inf ? m : m + log(exp(x1 - m) + exp(x2 - m))
end

export logsumexp

# mcmc
include("kernel_dsl.jl")
include("mh.jl")
include("hmc.jl")
include("mala.jl")
include("elliptical_slice.jl")
include("involution_dsl.jl")

include("importance.jl")
include("particle_filter.jl")
include("map_optimize.jl")
include("train.jl")
include("variational.jl")


text:
module CallStack

export @callstack

using Cassette: Cassette

mutable struct CallStackState
    output::IO
    serial::Int
    depth::Int
    CallStackState(output::IO) = new(output, 0, 0)
end

const INDENT = "  "

Cassette.@context CallStackContext

function Cassette.prehook(ctx::CallStackContext, f, args...)
    state = ctx.metadata
    if !(f isa Core.Builtin)
        state.serial += 1
        method = first(methods(f, typeof.(args)))
        println(
            # output stream
            state.output,
            # serial
            lpad(state.serial, 4), "  ",
            # indent
            INDENT ^ state.depth,
            # method info
            method,
        )
    end
    state.depth += 1
end

function Cassette.posthook(ctx::CallStackContext, f, args...)
    state = ctx.metadata
    state.depth -= 1
end

function gencallstack(fcall; output = stderr)
    quote
        state = CallStackState($(esc(output)))
        Cassette.overdub(CallStackContext(metadata = state), () -> $(esc(fcall)))
    end
end

macro callstack(fcall)
    gencallstack(fcall)
end

macro callstack(output, fcall)
    gencallstack(fcall, output = output)
end

end # module


text:
push!(LOAD_PATH,"./src");
using Documenter, EntanglementDist;

makedocs(format = :html, sitename = "Entanglement Distillation",
    pages = [
        "index.md",
        "intro.md",
        "optimize.md",
        "quantum.md",
        "schemes.md",
        "states.md",
        "examplePPT.md",
        "exampleKEXT.md",
        "exampleSeesaw.md"
        ]
)


text:
# on vérifie simplement que le constructeur par défaut fait son travail
n = 2
test_meta = MPCCModelMeta(n)

@test test_meta.nvar == n       # number of variables
@test test_meta.x0   == zeros(n)    # initial guess
@test size(test_meta.lvar) == size(zeros(n))   # vector of lower bounds
@test size(test_meta.uvar) == size(zeros(n))   # vector of upper bounds

@test test_meta.ncon == 0            # number of general constraints
@test test_meta.y0   == zeros(0)    # initial Lagrange multipliers
@test test_meta.lcon == zeros(0)    # vector of constraint lower bounds
@test test_meta.ucon == zeros(0)    # vector of constraint upper bounds

@test test_meta.ncc   == 0             # number of complementarity constraints
@test test_meta.lccG  == zeros(0)   # vector of constraint lower bounds of the complementarity constraint
@test test_meta.lccH  == zeros(0)   # vector of constraint upper bounds of the complementarity constraint

@test test_meta.minimize == true        # true if optimize == minimize
@test test_meta.nlo  == 1             # number of nonlinear objectives
@test test_meta.islp == false             # true if the problem is a linear program

try
    MPCCModelMeta(0)
    @test false
catch
    @test true
end

try
    show(test_meta)
    print(test_meta)
    @test true
catch
    @test false
end

@test MPCC.bound_constrained(test_meta)           == false
@test MPCC.unconstrained(test_meta)               == true
@test MPCC.linearly_constrained(test_meta)        == false
@test MPCC.inequality_constrained(test_meta)        == false
@test MPCC.equality_constrained(test_meta)        == false
@test complementarity_constrained(test_meta) == false


text:
using JustJoshing
using Test
using ForwardDiff: derivative

const EPS = 0.0001

@testset "Call option consistency" begin
    # I: see code

    # II: Price of call option monotonically decreases with strike price
    @test all(derivative.(x->C(1,0.5;K=x),rand(100_000)*10) .<= EPS)

    # III: Call option should be between S and S - Ke(-r(T-t)) for all inputs
    #
    # (if call option is more expensive than stock, it is clearly better
    # to just buy the stock; if it is cheaper than the difference between
    # the current value of the stock and the discounted strike price, we
    # can arbitrage by selling the stock now, buying the option, putting the
    # strike price in a risk free bond, and then exercising the option
    # (so we end up with the stock back and some profit))
    #
    # NB: Joshi refers to time-to-maturity, i.e. t_joshi = T-t :)
    @test [
        let (S,t,K,r,σ) = rand(5), T = t + rand()
            S + EPS >= C(S,t;T=T,K=K,r=r,σ=σ) >= (S - K*exp(-r*(T-t))) - EPS
        end
    for _ in 1:100] |> all

    # IV: Price of call option should monotonically increase with volatility
    @test all(derivative.(x->C(1,0.5;σ=x),rand(1000)) .>= -EPS)

    # V: if d=0 (always true in our model) then price should increase as we get closer to expiry
    @test all(derivative.(T->C(1,0;T=T),rand(1000)) .>= -EPS)

    # VI: Call option should be a convex function of strike - i.e. gradient is an increasing function of strike
    #   TODO: work out if there's a more ergonomic way of calculating higher order derivatives
    @test all(derivative.(x->derivative.(K->C(1,0;K=K),x),rand(1000)) .>= -EPS)

    # VII: The price of a call-spread should approximate the price of a digital-call option
    # TODO - NOT IMPLEMENTED - 1) find out what call-spread is and price it; 2) find out to calculate price of digital call (it was really simple IIRC); 3) check they are roughly the same

    # VIII: The price of a digital-call option plus a digital-put option is equal to the price of a zero-coupon bond
    # TODO - Not implemented: see VII
end

@testset "Monte-Carlo validation" begin
    # Validation with Monte-Carlo
    # S-K: forward contract
    let S_t=100, t=0.01, K=100, r=0.02, σ=0.05, T=1, trials=100_000_000
        exact = S_t - K*exp(-r*(T-t))

        mc = mc_pricer(S->S-K;S_0=S_t,t=t,T=T,K=K,r=r,σ=σ,trials=trials) |> first

        @test mc ≈ exact rtol=1e-3
    end

    # max(S-K,0): call option
    let S_t=100, t=0.01, K=100, r=0.02, σ=0.05, T=1, trials=100_000_000
        bs = C(S_t,t;T=T,K=K,r=r,σ=σ)

        mc = mc_pricer(S->max(S-K,0);S_0=S_t,t=t,T=T,K=S_t,r=r,σ=σ,trials=trials) |> first

        @test mc ≈ bs rtol=1e-3
    end

    # Int(S>E): binary option
    let S_t=100, t=0.01, E=100, r=0.02, σ=0.05, T=1, trials=100_000_000
        exact = binary(S_t,t;E=E,r=r,σ=σ,T=T)

        mc = mc_pricer(S->Int(S>E);S_0=S_t,t=t,T=T,K=E,r=r,σ=σ,trials=trials) |> first

        @test mc ≈ exact rtol=1e-3
    end
end


text:
module SnpParser

using LinearAlgebra;
using Plots;

include("parser.jl")
include("utils.jl")
include("plotter.jl")

export TouchstoneSnP, plot_snp

end # module


text:
using StipplePlotly
using Test

@testset "StipplePlotly.jl" begin
    # Write your tests here.
end


text:
"""
    struct Table{T,D,P} <: AbstractVector{Vector{T}}
      data::D
      ptrs::P
    end

Type representing a list of lists (i.e., a table) in
compressed format.
"""
struct Table{T,D,P} <: AbstractVector{Vector{T}}
  data::D
  ptrs::P

  @doc """
  """
  function Table(data::AbstractVector,ptrs::AbstractVector)
    T = eltype(data)
    D = typeof(data)
    P = typeof(ptrs)
    new{T,D,P}(data,ptrs)
  end
end

size(a::Table) = (length(a.ptrs)-1,)

IndexStyle(::Type{<:Table}) = IndexLinear()

function array_cache(a::Table)
  if length(a.ptrs) > 1
    pini = a.ptrs[1]
    pend = a.ptrs[2]
  else
    pend = a.ptrs[1]
    pini = pend
  end
  T = eltype(a.data)
  l = pend - pini
  r = zeros(T,l)
  CachedArray(r)
end

function getindex!(c,a::Table,i::Integer)
  pini = a.ptrs[i]
  l = a.ptrs[i+1] - pini
  setsize!(c,(l,))
  pini -= 1
  r = c.array
  for j in 1:l
     @inbounds r[j] = a.data[pini+j]
  end
  r
end

function getindex(a::Table,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end


text:
abstract type PrimitiveMeasure <: AbstractMeasure end

struct LebesgueMeasure <: PrimitiveMeasure end

export CountingMeasure

struct CountingMeasure <: PrimitiveMeasure end

struct TrivialMeasure <: PrimitiveMeasure end

logdensity_def(::PrimitiveMeasure, x) = static(0.0)

basemeasure(μ::PrimitiveMeasure) = μ

basemeasure_type(μ::M) where {M<:PrimitiveMeasure} = M

basemeasure_type(::Type{M}) where {M<:PrimitiveMeasure} = M

basemeasure_depth(::Type{M}) where {M<:PrimitiveMeasure} = static(0)


text:
include("loglinspace.jl")
include("../src/ttv.jl")
maxabs(x) = maximum(abs.(x))
if VERSION >= v"0.7"
  using Test
  using LinearAlgebra
  using Statistics
  using DelimitedFiles
else
  using Base.Test
end

include("test_kepler_init.jl")
include("test_init_nbody.jl")
include("test_elliptic_derivative.jl")
include("test_keplerij.jl")
include("test_kepler_drift.jl")
include("test_kepler_driftij.jl")
include("test_kepler_drift_gamma.jl")
include("test_kepler_driftij_gamma.jl")
include("test_kickfast.jl")
include("test_phisalpha.jl")
include("test_dh17.jl")
include("test_ah18.jl")
include("test_findtransit3.jl")
include("test_ttv_cartesian.jl")
include("test_ttvbv_cartesian.jl")
include("test_ttv_elements.jl")
include("test_ttvbv_elements.jl")


text:
# Julia Interface to Pips-NLP
# Feng Qiang, Cosmin G. Petra
# Argonne National Laboratory, 2016

# Wrapper for the parallel/structured PIPS-NLP interface
module PipsNlpSolver
import MPI
using Libdl

try
  sharedLib=ENV["PIPS_NLP_PAR_SHARED_LIB"]

  # Explicitly check if the file exists. dlopen sometimes does not throw an
  # error for invalid filenames, resulting in a seg fault)
  if(!isfile(sharedLib))
    error(string("The specified shared library ([", sharedLib, "]) does not exist"))
  end
  global const libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_PAR_SHARED_LIB",""))
catch
  @warn("Could not load PIPS-NLP shared library. Make sure the ENV variable 'PIPS_NLP_PAR_SHARED_LIB' points to its location, usually in the PIPS repo at PIPS/build_pips/PIPS-NLP/libparpipsnlp.so")
  rethrow()
end
abstract type ModelInterface end
#######################
mutable struct PipsModel <: ModelInterface
    sense::Symbol
    status::Int
    nscen::Int
    rowmap::Dict{Int,Int}
    colmap::Dict{Int,Int}
    eq_rowmap::Dict{Int,Int}
    inq_rowmap::Dict{Int,Int}

    get_num_scen::Function
    get_sense::Function

    get_status::Function
    get_num_rows::Function
    get_num_cols::Function
    get_num_eq_cons::Function
    get_num_ineq_cons::Function

    set_status::Function
    set_num_rows::Function
    set_num_cols::Function
    set_num_eq_cons::Function
    set_num_ineq_cons::Function

    str_init_x0::Function
    str_prob_info::Function
    str_eval_f::Function
    str_eval_g::Function
    str_eval_grad_f::Function
    str_eval_jac_g::Function
    str_eval_h::Function
    str_write_solution::Function

    function PipsModel(sense::Symbol,status::Int,nscen::Int, str_init_x0, str_prob_info, str_eval_f, str_eval_g,str_eval_grad_f,str_eval_jac_g, str_eval_h, str_write_solution)
        instance = new(sense,status,nscen,Dict{Int,Int}(),Dict{Int,Int}(),Dict{Int,Int}(),Dict{Int,Int}())
        instance.str_init_x0 = str_init_x0
        instance.str_prob_info = str_prob_info
        instance.str_eval_f = str_eval_f
        instance.str_eval_g = str_eval_g
        instance.str_eval_grad_f = str_eval_grad_f
        instance.str_eval_jac_g = str_eval_jac_g
        instance.str_eval_h = str_eval_h
	    instance.str_write_solution = str_write_solution

        instance.get_num_scen = function()
            return instance.nscen
        end
        instance.get_sense = function()
            return instance.sense
        end
        instance.get_status = function()
            return instance.status
        end
        instance.get_num_rows = function(id::Integer)
            return instance.rowmap[id]
        end
        instance.get_num_cols = function(id::Integer)
            return instance.colmap[id]
        end
        instance.get_num_eq_cons = function(id::Integer)
            return instance.eq_rowmap[id]
        end
        instance.get_num_ineq_cons = function(id::Integer)
            return instance.inq_rowmap[id]
        end
        instance.set_status = function(s::Integer)
            instance.status = s
        end
        instance.set_num_rows = function(id::Integer, v::Integer)
            return instance.rowmap[id] = v
        end
        instance.set_num_cols = function(id::Integer, v::Integer)
            return instance.colmap[id] = v
        end
        instance.set_num_eq_cons = function(id::Integer, v::Integer)
            return instance.eq_rowmap[id] = v
        end
        instance.set_num_ineq_cons = function(id::Integer, v::Integer)
            return instance.inq_rowmap[id] = v
        end
        return instance
    end
end

#######################
mutable struct PipsNlpProblemStruct
    ref::Ptr{Nothing}
    model::ModelInterface
    comm::MPI.Comm
    prof::Bool

    n_iter::Int
    t_jl_init_x0::Float64
    t_jl_str_prob_info::Float64
    t_jl_eval_f::Float64
    t_jl_eval_g::Float64
    t_jl_eval_grad_f::Float64

    t_jl_eval_jac_g::Float64
    t_jl_str_eval_jac_g::Float64
    t_jl_eval_h::Float64
    t_jl_str_eval_h::Float64
    t_jl_write_solution::Float64

    t_jl_str_total::Float64
    t_jl_eval_total::Float64

    function PipsNlpProblemStruct(comm, model, prof)
        prob = new(C_NULL, model, comm, prof,-3
            ,0.0,0.0,0.0,0.0,0.0
            ,0.0,0.0,0.0,0.0,0.0
            ,0.0,0.0
            )
        #finalizer(prob, freeProblemStruct)
        finalizer(freeProblemStruct, prob)
        return prob
    end
end

struct CallBackData
	prob::Ptr{Nothing}
	row_node_id::Cint
    col_node_id::Cint
    flag::Cint  #this wrapper ignore this flag as it is only for problems without linking constraint
end

export  ModelInterface, PipsModel,
        createProblemStruct, solveProblemStruct

###########################################################################
# Callback wrappers
###########################################################################

function str_init_x0_wrapper(x0_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    data = unsafe_load(cbd)
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    @assert(rowid == colid)
    n0 = prob.model.get_num_cols(colid)
    x0 = unsafe_wrap(Array,x0_ptr,n0)

    # if prob.prof
    #     t1 = time()
    # end
    t = @elapsed prob.model.str_init_x0(colid,x0)
    if prob.prof
        prob.t_jl_init_x0 += t
    end
    return Int32(1)
end

# prob info (prob_info)
function str_prob_info_wrapper(n_ptr::Ptr{Cint}, col_lb_ptr::Ptr{Float64}, col_ub_ptr::Ptr{Float64}, m_ptr::Ptr{Cint}, row_lb_ptr::Ptr{Float64}, row_ub_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})

    data = unsafe_load(cbd)
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    flag = data.flag
    @assert(rowid == colid)

	mode = (col_lb_ptr == C_NULL&&col_ub_ptr==C_NULL&&row_lb_ptr==C_NULL&&row_ub_ptr==C_NULL) ? (:Structure) : (:Values)
    #@show flag
    if flag != 1
        #@show mode
    	if(mode==:Structure)
            col_lb = unsafe_wrap(Array,col_lb_ptr,0)
            col_ub = unsafe_wrap(Array,col_ub_ptr,0)
            row_lb = unsafe_wrap(Array,row_lb_ptr,0)
            row_ub = unsafe_wrap(Array,row_ub_ptr,0)
            # if prob.prof
            #     tic()
            # end
    		t = @elapsed (n,m) = prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            if prob.prof
                prob.t_jl_str_prob_info += t
            end

    		unsafe_store!(n_ptr,convert(Cint,n)::Cint)
    		unsafe_store!(m_ptr,convert(Cint,m)::Cint)
            # @show typeof(colid), typeof(m)
    		prob.model.set_num_rows(colid, m)
    		prob.model.set_num_cols(colid, n)
    	else
    		n = unsafe_load(n_ptr)
    		m = unsafe_load(m_ptr)
    		col_lb = unsafe_wrap(Array,col_lb_ptr,n)
    		col_ub = unsafe_wrap(Array,col_ub_ptr,n)
    		row_lb = unsafe_wrap(Array,row_lb_ptr,m)
    		row_ub = unsafe_wrap(Array,row_ub_ptr,m)

            # if prob.prof
            #     tic()
            # end
    		t = @elapsed prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            if prob.prof
                prob.t_jl_str_prob_info += t
            end

    		neq = 0
    		nineq = 0
    		for i = 1:length(row_lb)
    			if row_lb[i] == row_ub[i]
    				neq += 1
    			else
    				nineq += 1
    			end
    		end
    		@assert(neq+nineq == length(row_lb) == m)
    		prob.model.set_num_eq_cons(colid,neq)
    		prob.model.set_num_ineq_cons(colid,nineq)
    	end
    else
	#println("prob_info: ",mode,flag)
        @assert flag == 1
        if mode == :Structure
            col_lb = unsafe_wrap(Array,col_lb_ptr,0)
            col_ub = unsafe_wrap(Array,col_ub_ptr,0)
            row_lb = unsafe_wrap(Array,row_lb_ptr,0)
            row_ub = unsafe_wrap(Array,row_ub_ptr,0)
            (n,m) = prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
            # @show n,m
            unsafe_store!(m_ptr,convert(Cint,m)::Cint)
        else
	# edited by yankai
            m = unsafe_load(m_ptr)
	        col_lb = unsafe_wrap(Array,col_lb_ptr,0)
            col_ub = unsafe_wrap(Array,col_ub_ptr,0)
            row_lb = unsafe_wrap(Array,row_lb_ptr,m)
            row_ub = unsafe_wrap(Array,row_ub_ptr,m)
            prob.model.str_prob_info(colid,flag,mode,col_lb,col_ub,row_lb,row_ub)
        end
    end
    # @show "exit  julia - str_prob_info_wrapper "
	return Int32(1)
end
# Objective (eval_f)
function str_eval_f_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, obj_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia - eval_f_wrapper "
    data = unsafe_load(cbd)
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    @assert(rowid == colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(colid)
    # Calculate the new objective
    x0 = unsafe_wrap(Array,x0_ptr, n0)
    x1 = unsafe_wrap(Array,x1_ptr, n1)

    # if prob.prof
    #     tic()
    # end
    t = @elapsed new_obj = convert(Float64, prob.model.str_eval_f(colid,x0,x1))::Float64
    if prob.prof
        prob.t_jl_eval_f += t
    end
    # Fill out the pointer
    unsafe_store!(obj_ptr, new_obj)
    # Done
    # @show "exit julia - eval_f_wrapper "
    return Int32(1)
end

# Constraints (eval_g)
function str_eval_g_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, eq_g_ptr::Ptr{Float64}, inq_g_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia - eval_g_wrapper "
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    @assert(rowid == colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(colid)
    x0 = unsafe_wrap(Array,x0_ptr, n0)
    x1 = unsafe_wrap(Array,x1_ptr, n1)
    # Calculate the new constraint values

    #println("rowid",rowid)
    #println(prob.model)
    # if MPI.Comm_rank(MPI.COMM_WORLD) == 1
    #     println("rank: ",MPI.Comm_rank(MPI.COMM_WORLD),"rowid: ",rowid)
    # end
    neq = prob.model.get_num_eq_cons(rowid)
    nineq = prob.model.get_num_ineq_cons(rowid)
    new_eq_g = unsafe_wrap(Array,eq_g_ptr,neq)
    new_inq_g = unsafe_wrap(Array,inq_g_ptr, nineq)

    # if prob.prof
    #     tic()
    # end
    t = @elapsed prob.model.str_eval_g(colid,x0,x1,new_eq_g,new_inq_g)
    if prob.prof
        prob.t_jl_eval_g += t
    end
    # Done
    # @show " exit julia - eval_g_wrapper "
    return Int32(1)
end

# Objective gradient (eval_grad_f)
function str_eval_grad_f_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, grad_f_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    # @show " julia -  eval_grad_f_wrapper "
    # Extract Julia the problem from the pointer
    # @show cbd
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(rowid)
    # @show n0,n1
    x0 = unsafe_wrap(Array,x0_ptr, n0)
    x1 = unsafe_wrap(Array,x1_ptr, n1)
    # Calculate the gradient
    grad_len = prob.model.get_num_cols(colid)
    new_grad_f = unsafe_wrap(Array,grad_f_ptr, grad_len)

    # if prob.prof
    #     tic()
    # end
    t = @elapsed prob.model.str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
    if prob.prof
        prob.t_jl_eval_grad_f += t
    end
    if prob.model.get_sense() == :Max
        new_grad_f *= -1.0
    end
    # Done
    # @show " julia -  eval_grad_f_wrapper "
    return Int32(1)
end

# Jacobian (eval_jac_g)
function str_eval_jac_g_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64},
	e_nz_ptr::Ptr{Cint}, e_values_ptr::Ptr{Float64}, e_row_ptr::Ptr{Cint}, e_col_ptr::Ptr{Cint},
	i_nz_ptr::Ptr{Cint}, i_values_ptr::Ptr{Float64}, i_row_ptr::Ptr{Cint}, i_col_ptr::Ptr{Cint},
	cbd::Ptr{CallBackData}
	)
    # @show " julia -  eval_jac_g_wrapper "
    # Extract Julia the problem from the pointer
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    flag = data.flag
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(rowid) #we can do this because of 2-level and no linking constraint
    # @show n0, n1
    x0 = unsafe_wrap(Array,x0_ptr, n0)
    x1 = unsafe_wrap(Array,x1_ptr, n1)
    # @show x0
    # @show x1
    nrow = prob.model.get_num_rows(rowid)
    ncol = prob.model.get_num_cols(colid)
    #@show prob
    # Determine mode
    mode = (e_row_ptr == C_NULL &&e_col_ptr == C_NULL&&e_values_ptr == C_NULL && i_values_ptr == C_NULL&&i_row_ptr == C_NULL &&i_col_ptr == C_NULL) ? (:Structure) : (:Values)
    #@show mode
    if flag != 1
        if(mode == :Structure)
        	e_values = unsafe_wrap(Array,e_values_ptr,0)
    		e_colptr = unsafe_wrap(Array,e_col_ptr,0)
    		e_rowidx = unsafe_wrap(Array,e_row_ptr,0)
    		i_values = unsafe_wrap(Array,i_values_ptr,0)
    		i_colptr = unsafe_wrap(Array,i_col_ptr,0)
    		i_rowidx = unsafe_wrap(Array,i_row_ptr,0)
            # if prob.prof
            #     tic()
            # end
            t = @elapsed (e_nz,i_nz) = prob.model.str_eval_jac_g(rowid,colid,flag,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            if prob.prof
                prob.t_jl_str_eval_jac_g += t
            end
    		unsafe_store!(e_nz_ptr,convert(Cint,e_nz)::Cint)
    		unsafe_store!(i_nz_ptr,convert(Cint,i_nz)::Cint)
    		# @show "structure - ",(e_nz,i_nz)
        else
        	e_nz = unsafe_load(e_nz_ptr)
        	e_values = unsafe_wrap(Array,e_values_ptr,e_nz)
        	e_rowidx = unsafe_wrap(Array,e_row_ptr, e_nz)
        	e_colptr = unsafe_wrap(Array,e_col_ptr, ncol+1)
        	i_nz = unsafe_load(i_nz_ptr)
        	# @show "values - ",(e_nz,i_nz), ncol
        	i_values = unsafe_wrap(Array,i_values_ptr,i_nz)
        	i_rowidx = unsafe_wrap(Array,i_row_ptr, i_nz)
        	i_colptr = unsafe_wrap(Array,i_col_ptr, ncol+1)

            # if prob.prof
            #     tic()
            # end
        	t = @elapsed prob.model.str_eval_jac_g(rowid,colid,flag,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            if prob.prof
                prob.t_jl_eval_jac_g += t
            end
        end
    else
        @assert flag == 1
        if mode == :Structure
                e_values = unsafe_wrap(Array,e_values_ptr,0)
                e_colptr = unsafe_wrap(Array,e_col_ptr,0)
                e_rowidx = unsafe_wrap(Array,e_row_ptr,0)
                i_values = unsafe_wrap(Array,i_values_ptr,0)
                i_colptr = unsafe_wrap(Array,i_col_ptr,0)
                i_rowidx = unsafe_wrap(Array,i_row_ptr,0)
            # if prob.prof
            #     tic()
            # end
            t = @elapsed (e_nz,i_nz) = prob.model.str_eval_jac_g(rowid,colid,flag, x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            if prob.prof
                prob.t_jl_str_eval_jac_g += t
            end
            unsafe_store!(e_nz_ptr,convert(Cint,e_nz)::Cint)
            unsafe_store!(i_nz_ptr,convert(Cint,i_nz)::Cint)
        else
            e_nz = unsafe_load(e_nz_ptr)
            e_values = unsafe_wrap(Array,e_values_ptr,e_nz)
            e_rowidx = unsafe_wrap(Array,e_row_ptr, e_nz)
            e_colptr = unsafe_wrap(Array,e_col_ptr, ncol+1)
            i_nz = unsafe_load(i_nz_ptr)
            i_values = unsafe_wrap(Array,i_values_ptr,i_nz)
            i_rowidx = unsafe_wrap(Array,i_row_ptr, i_nz)
            i_colptr = unsafe_wrap(Array,i_col_ptr, ncol+1)
            # if prob.prof
            #     tic()
            # end
            t = @elapsed prob.model.str_eval_jac_g(rowid,colid,flag, x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
            if prob.prof
                prob.t_jl_eval_jac_g += t
            end
        end
    end
    # Done
    #@show "exit julia -  eval_jac_g_wrapper "
    return Int32(1)
end

# Hessian
function str_eval_h_wrapper(x0_ptr::Ptr{Float64}, x1_ptr::Ptr{Float64}, lambda_ptr::Ptr{Float64}, nz_ptr::Ptr{Cint}, values_ptr::Ptr{Float64}, row_ptr::Ptr{Cint}, col_ptr::Ptr{Cint}, cbd::Ptr{CallBackData})
    # @show " julia - eval_h_wrapper "
    # Extract Julia the problem from the pointer
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    # @show prob.prof
    rowid = data.row_node_id
    colid = data.col_node_id

    high = max(rowid,colid)
    low  = min(rowid,colid)
    n0 = prob.model.get_num_cols(0)
    n1 = prob.model.get_num_cols(high)
    x0 = unsafe_wrap(Array,x0_ptr, n0)
    x1 = unsafe_wrap(Array,x1_ptr, n1)

    ncol = prob.model.get_num_cols(low)
    g0 = prob.model.get_num_rows(high)

    lambda = unsafe_wrap(Array,lambda_ptr, g0)
    obj_factor = 1.0
    if prob.model.get_sense() == :Max
        obj_factor *= -1.0
    end
    # Did the user specify a Hessian
    mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
    if(mode == :Structure)
    	values = unsafe_wrap(Array,values_ptr,0)
		colptr = unsafe_wrap(Array,col_ptr,0)
		rowidx = unsafe_wrap(Array,row_ptr,0)
        # if prob.prof
        #     tic()
        # end
		t = @elapsed nz = prob.model.str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        if prob.prof
            prob.t_jl_str_eval_h += t
        end
		unsafe_store!(nz_ptr,convert(Cint,nz)::Cint)
		# @show "structure - ", nz
    else
    	nz = unsafe_load(nz_ptr)
    	values = unsafe_wrap(Array,values_ptr, nz)
    	rowidx = unsafe_wrap(Array,row_ptr, nz)
    	colptr = unsafe_wrap(Array,col_ptr, ncol+1)
    	# @show "value - ", nz
        # if prob.prof
        #     tic()
        # end
    	t = @elapsed prob.model.str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        if prob.prof
            prob.t_jl_eval_h += t
        end
    end
    # Done
    if prob.prof
        if rowid == colid == 0
            prob.n_iter += 1
            if prob.n_iter == 0
                prob.t_jl_str_total = t_reset(prob)
            end
        end
    end
    # @show "exit  julia - eval_h_wrapper "
    return Int32(1)
end

#write solution
function str_write_solution_wrapper(x_ptr::Ptr{Float64}, y_eq_ptr::Ptr{Float64}, y_ieq_ptr::Ptr{Float64}, cbd::Ptr{CallBackData})
    data = unsafe_load(cbd)
    # @show data
    userdata = data.prob
    prob = unsafe_pointer_to_objref(userdata)::PipsNlpProblemStruct
    rowid = data.row_node_id
    colid = data.col_node_id
    @assert rowid == colid

    nx = prob.model.get_num_cols(rowid)
    neq = prob.model.get_num_eq_cons(rowid)
    nieq = prob.model.get_num_ineq_cons(rowid)
    x = unsafe_wrap(Array,x_ptr, nx)
    y_eq = unsafe_wrap(Array,y_eq_ptr,neq)
    y_ieq = unsafe_wrap(Array,y_ieq_ptr,nieq)
    # if prob.prof
    #     tic()
    # end
    t = @elapsed prob.model.str_write_solution(rowid,x,y_eq,y_ieq)
    if prob.prof
        prob.t_jl_write_solution += t
    end

    return Int32(1)
end
###########################################################################
# C function wrappers
###########################################################################

function createProblemStruct(comm::MPI.Comm, model::ModelInterface, prof::Bool)
	# println(" createProblemStruct  -- julia")
	str_init_x0_cb = @cfunction(str_init_x0_wrapper, Cint, (Ptr{Float64}, Ptr{CallBackData}))
    str_prob_info_cb = @cfunction(str_prob_info_wrapper, Cint, (Ptr{Cint}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_f_cb = @cfunction(str_eval_f_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_g_cb = @cfunction(str_eval_g_wrapper,Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_grad_f_cb = @cfunction(str_eval_grad_f_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}) )
    str_eval_jac_g_cb = @cfunction(str_eval_jac_g_wrapper, Cint, (Ptr{Float64}, Ptr{Float64},
    	Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint},
    	Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint},
    	Ptr{CallBackData}))
    str_eval_h_cb = @cfunction(str_eval_h_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}, Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{CallBackData}))
    str_write_solution_cb = @cfunction(str_write_solution_wrapper, Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{CallBackData}))

    # println(" callback created ")
    prob = PipsNlpProblemStruct(comm, model, prof)
    # @show prob
    libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_PAR_SHARED_LIB",""))
    ret = ccall(Libdl.dlsym(libparpipsnlp,:CreatePipsNlpProblemStruct),Ptr{Nothing},
            (MPI.Comm,
            Cint, Ptr{Nothing}, Ptr{Nothing},Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing},Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing},Any
            # ,Ptr{Void}, Ptr{Void}  #comply with link interface from yankai
            ),
            comm,
            model.get_num_scen(),
            str_init_x0_cb,
            str_prob_info_cb,
            str_eval_f_cb,
            str_eval_g_cb,
            str_eval_grad_f_cb,
            str_eval_jac_g_cb,
            str_eval_h_cb,
            str_write_solution_cb,
            prob
            # ,Ptr{Void}(0), Ptr{Void}(0)
            )
    if ret == C_NULL
        error("PIPS-NLP: Failed to construct problem.")
    else
        prob.ref = ret
    end
    return prob
end

function solveProblemStruct(prob::PipsNlpProblemStruct)
    # println("solveProblemStruct - julia")
    # @show prob

    libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_PAR_SHARED_LIB",""))
    ret = ccall(Libdl.dlsym(libparpipsnlp,:PipsNlpSolveStruct), Cint,
            (Ptr{Nothing},),
            prob.ref)
    # @show ret
    prob.model.set_status(Int(ret))

    prob.t_jl_eval_total = report_total_now(prob)
    # @show prob
    return prob.model.get_status()
end

function freeProblemStruct(prob::PipsNlpProblemStruct)
    # @show "freeProblemStruct"
    libparpipsnlp=Libdl.dlopen(get(ENV,"PIPS_NLP_PAR_SHARED_LIB",""))
    ret = ccall(Libdl.dlsym(libparpipsnlp,:FreePipsNlpProblemStruct),
            Nothing, (Ptr{Nothing},),
            prob.ref)
    # @show ret
    return ret
end

function report_total_now(prob::PipsNlpProblemStruct)
    total = 0.0
    total += prob.t_jl_init_x0
    total += prob.t_jl_str_prob_info
    total += prob.t_jl_eval_f
    total += prob.t_jl_eval_g
    total += prob.t_jl_eval_grad_f
    total += prob.t_jl_eval_jac_g
    total += prob.t_jl_str_eval_jac_g
    total += prob.t_jl_eval_h
    total += prob.t_jl_str_eval_h
    total += prob.t_jl_write_solution
    return total
end

function t_reset(prob::PipsNlpProblemStruct)
    total = report_total_now(prob)
    prob.t_jl_init_x0  = 0.0
    prob.t_jl_str_prob_info  = 0.0
    prob.t_jl_eval_f = 0.0
    prob.t_jl_eval_g = 0.0
    prob.t_jl_eval_grad_f = 0.0
    prob.t_jl_eval_jac_g = 0.0
    prob.t_jl_str_eval_jac_g = 0.0
    prob.t_jl_eval_h = 0.0
    prob.t_jl_str_eval_h = 0.0
    prob.t_jl_write_solution = 0.0
    return total
end

end


text:
module PrePostCall

import MacroTools: splitdef,combinedef,postwalk

export @pre,@post

function get_qargs_for_call(x::Array{Any,1})
    args = map(x) do a
        isa(a,Expr) && return Expr(:quote,:($(a.args[1])))
        return Expr(:quote,:($a))
    end
end

function get_args_for_call(x::Array{Any,1})
    args = map(x) do a
        isa(a,Expr) && return (a.args[1])
        return a
    end
end

"""
    @pre function name ...

Create a macro `@name [[variable] variable ...] function other ...` for inserting a call to `name` before each call to `other`.

`variable` defines the variable names passed to `name`. If `variable` is omitted the names of the attributes of `name` are used.
If multiple variables are given `name` is called on each of them, for example:

- `@name x y z` is just a short notation for `@name x @name y @name z` and calls `name(x)`, `name(y)`, `name(z)`
- `@name x,y,z` calls `name(x,y,z)`

# Examples
```jldoctest
julia> @pre function nonzero(x::Int)
           @assert x!=0
       end
@nonzero (macro with 2 methods)

julia> @nonzero @nonzero y function foo(x::Int,y::Int)
           x+y
       end
foo (generic function with 1 method)
```
The outer `@nonzero` uses `x` as the attribute due to the function definition of `nonzero`.
```jldoctest
julia> foo(1,2)
3

julia> foo(1,0)
ERROR: AssertionError: x != 0

julia> foo(0,1)
ERROR: AssertionError: x != 0
```

"""
macro pre(cfun)
    def = splitdef(cfun)
    mname = QuoteNode(Symbol("@$(def[:name])"))
    return esc(quote
               # Define the actual function
               $cfun
               # Define the macro without explicit variable definition
               macro $(def[:name])(fun)
                   # If the head is a macro call evaluate this macro first.
                   # In order to function properly the pre and post macros
                   # must be evaluated from insided to outside
                   if fun.head == :macrocall
                       fun = macroexpand(@__MODULE__,fun)
                   end
                   _def = PrePostCall.splitdef(fun)
                   # Insert the function call as first element in the body
                   insert!(_def[:body].args,1,Expr(:call,$(def[:name]),$(PrePostCall.get_qargs_for_call(def[:args])...)))
                   esc(PrePostCall.combinedef(_def))
               end
               # Define the macro with explicit variable definition
               macro $(def[:name])(arg,fun)
                   # If multiple explicit variables are defined expand them on the function call
                   arg = isa(arg,Symbol) ? arg : Expr(:...,arg)
                   if fun.head == :macrocall
                       fun = macroexpand(@__MODULE__,fun)
                   end
                   _def = PrePostCall.splitdef(fun)
                   # The only difference to the macro without the explicit definition is,
                   # that the passed variables are used for the function call
                   insert!(_def[:body].args,1,Expr(:call,$(def[:name]),:($arg)))
                   esc(PrePostCall.combinedef(_def))
               end
               macro $(def[:name])(args...)
                   # The second attribute to a macrocall expression is always a
                   # LineNumberNode. Here I replaced it by nothing.
                   clast = Expr(:macrocall,$mname, nothing,args[end-1], args[end])
                   for i in length(args)-2:-1:1
                       clast = Expr(:macrocall,$mname, nothing,args[i], clast)
                   end
                   return esc(clast)
               end
               end)
end

"""
    @post function name ...

Create a macro `@name [[variable] variable ...] function other ...` for inserting a call to `name` after each call to `other`.

`variable` defines the variable names passed to `name`. If `variable` is omitted, `name` is called on the return argument of `other`.
If `variable` is used, the call to `other` is inserted before each `return`, or if non present, as last expression in `other`.
If multiple variables are given `name` is called on each of them, for example:

- `@name x y z` is just a short notation for `@name x @name y @name z` and calls `name(x)`, `name(y)`, `name(z)`
- `@name x,y,z` calls `name(x,y,z)`

# Examples
```jldoctest
julia> @post function nonzero(x::Int)
           @assert x!=0
       end
@nonzero (macro with 2 methods)

julia> @nonzero function foo(x::Int,y::Int)
           x*y
       end
foo (generic function with 1 method)

julia> foo(1,2)
2

julia> foo(1,0)
ERROR: AssertionError: x != 0

julia> @nonzero @nonzero a function foo(x::Int,y::Int)
           a = x-1
           return a*y
       end
foo (generic function with 1 method)

julia> foo(1,2)
ERROR: AssertionError: x != 0
```
Failes because `a` must be nonzero
```jldoctest
julia> foo(2,2)
2

julia> foo(2,0)
ERROR: AssertionError: x != 0
```
Failes because the return value must be nonzero

"""
macro post(cfun)
    def = splitdef(cfun)
    # In case the check function has multiple arguments the return value
    # must be expanded
    # This is only used for the macro without explicit variable definition.
    retExpr = length(def[:args])==1 ? :(Expr(:call,$(def[:name]),:($ret))) : :(Expr(:call,$(def[:name]),Expr(:...,:($ret))))
    mname = QuoteNode(Symbol("@$(def[:name])"))
    return esc(quote
               $cfun
               # Define the macro without explicit variable definition
               # This macro nests the function call insided a new equally named funciton.
               macro $(def[:name])(fun)
                   # If the head is a macro call evaluate this macro first.
                   # In order to function properly the pre and post macros
                   # must be evaluated from insided to outside
                   if fun.head == :macrocall
                       fun = macroexpand(@__MODULE__,fun)
                   end
                   # Dict for the new function
                   _def = PrePostCall.splitdef(fun)
                   # Dict for the original function
                   _fun = copy(_def)
                   # Get a unique name for the return value variable
                   ret = gensym("ret")
                   # Get a unique name for the now nested original function
                   name = gensym(_def[:name])
                   _fun[:name] = name
                   # Define the function body of the new function
                   _def[:body] = Expr(:block,
                                      # Insert original function now as local and renamed
                                      PrePostCall.combinedef(_fun),
                                      # Store the return value in ret
                                      Expr(Symbol("="),:($ret),Expr(:call,:($name),:($(PrePostCall.get_args_for_call(_def[:args])...)))),
                                      # Evaluate the ckeck function
                                      $retExpr,
                                      # Return the return value
                                      Expr(:return,:($ret))
                                      )
                   esc(PrePostCall.combinedef(_def))
               end
               # Define the macro with explicit variable definition
               # This macro add the check function call before each return statement,
               # or if none where found as the last expression of the original function
               macro $(def[:name])(arg,fun)
                   # If multiple explicit variables are defined expand them on the function call
                   arg = isa(arg,Symbol) ? arg : Expr(:...,arg)
                   if fun.head == :macrocall
                       fun = macroexpand(@__MODULE__,fun)
                   end
                   _def = PrePostCall.splitdef(fun)
                   # Counter to check if the function has return statements
                   found = 0
                   # Walk through the expression an search for return statements
                   _def[:body] = PrePostCall.postwalk(_def[:body]) do l
                       # Add the check function call before the return statement
                       if isa(l,Expr) && l.head == :return
                           found += 1
                           return Expr(:block,Expr(:call,$(def[:name]),:($arg)),l)
                       end
                       return l
                   end
                   # No return statement was found add the call after the last expression
                   if found == 0
                       push!(_def[:body].args,Expr(:block,Expr(:call,$(def[:name]),:($arg)),:nothing))
                   end
                   esc(PrePostCall.combinedef(_def))
               end
               macro $(def[:name])(args...)
                   # The second attribute to a macrocall expression is always a
                   # LineNumberNode. Here I replaced it by nothing.
                   clast = Expr(:macrocall,$mname, nothing,args[end-1], args[end])
                   for i in length(args)-2:-1:1
                       clast = Expr(:macrocall,$mname, nothing,args[i], clast)
                   end
                   return esc(clast)
               end
               end)
end

end # module


text:
#=
#define FMPR_RND_DOWN  0    RoundingMode{:Down}()
#define FMPR_RND_UP    1
#define FMPR_RND_FLOOR 2
#define FMPR_RND_CEIL  3
#define FMPR_RND_NEAR  4
=#

import Base.Rounding: rounding_raw, setrounding_raw, rounding, setrounding;

const ARB_ROUNDING_MODE = Ref{Cint}(0)

to_arb(::RoundingMode{:Nearest}) = Cint(0)
to_arb(::RoundingMode{:ToZero}) = Cint(1)
to_arb(::RoundingMode{:Up}) = Cint(2)
to_arb(::RoundingMode{:Down}) = Cint(3)
to_arb(::RoundingMode{:FromZero}) = Cint(4)

function from_arb(c::Integer)
    if c == 0
        return RoundNearest
    elseif c == 1
        return RoundToZero
    elseif c == 2
        return RoundUp
    elseif c == 3
        return RoundDown
    elseif c == 4
        return RoundFromZero
    else
        throw(ArgumentError("invalid MPFR rounding mode code: $c"))
    end
    RoundingMode(c)
end

rounding_raw(::Type{ArbFloat}) = ARB_ROUNDING_MODE[]
setrounding_raw(::Type{ArbFloat},i::Integer) = ARB_ROUNDING_MODE[] = (Cint)(i)

rounding(::Type{ArbFloat}) = from_arb(rounding_raw(ArbFloat))
setrounding(::Type{ArbFloat},r::RoundingMode) = setrounding_raw(ArbFloat,to_arb(r))

rounding_raw{P}(::Type{ArbFloat{P}}) = ARB_ROUNDING_MODE[]
setrounding_raw{P}(::Type{ArbFloat{P}},i::Integer) = ARB_ROUNDING_MODE[] = (Cint)(i)

rounding{P}(::Type{ArbFloat{P}}) = from_arb(rounding_raw(ArbFloat))
setrounding{P}(::Type{ArbFloat{P}},r::RoundingMode) = setrounding_raw(ArbFloat,to_arb(r))


text:
module EmpiricalPseudopotential

include("Database.jl")
using .Database
const db = Database

# Constant value
const HBAR = 1.054571800e-34
const M = 9.10938356e-31
const Q = 1.60217662e-19

using DataStructures: DefaultDict
using LinearAlgebra
using Plots
plotlyjs()

const specialpaths = Dict(
    "sc" => "GXMGRX",
    "fcc" => "GXWKGLUWLK",
    "bcc" => "GHNGPH",
)

const specialpoints = Dict(
    "sc" => Dict(
        "G" => [0 0 0],
        "M" => [1/2 1/2 0],
        "R" => [1/2 1/2 1/2],
        "X" => [0 1/2 0],
    ),
    "fcc" => Dict(
        "G" => [0 0 0],
        "K" => [3/8 3/8 3/4],
        "L" => [1/2 1/2 1/2],
        "U" => [5/8 1/4 5/8],
        "W" => [1/2 1/4 3/4],
        "X" => [1/2 0 1/2],
    ),
    "bcc" => Dict(
        "G" => [0 0 0],
        "H" => [1/2 -1/2 1/2],
        "P" => [1/4 1/4 1/4],
        "N" => [0 0 1/2],
    ),
)

struct Magnitude
    m1::Float64
    m2::Float64
    m3::Float64
    function Magnitude(i, j, t)
        new(i, j, t)
    end
end

const magnitude = [ Magnitude(i, j, t) for i in -3:3 for j in -3:3 for t in -3:3 ]

# Crystal Lattice Structure
struct BasicLatticeVector
    b1::Matrix{Float64}
    b2::Matrix{Float64}
    b3::Matrix{Float64}
    function BasicLatticeVector(latticetype="fcc")
        if latticetype == "fcc"
            return new(
                [-1 1 1],
                [1 -1 1],
                [1 1 -1]
            )
        elseif latticetype == "bcc"
            return new(
                [0 1 1],
                [1 0 1],
                [1 1 0]
            )
        elseif latticetype == "sc"
            return new(
                [1 0 0],
                [0 1 0],
                [0 0 1]
            )
        else
            println("Invalid crystalstructure")
            exit()
        end
    end
end

const latticetypes = ["fcc", "bcc", "sc"]

# Reciprocal lattice vector
struct Gm
    magnitude::Float64
    vector::Matrix{Float64}
    function Gm(b::BasicLatticeVector, m::Magnitude)
        gm = b.b1 * m.m1 + b.b2 * m.m2 + b.b3 * m.m3
        n = gm[1]^2 + gm[2]^2 + gm[3]^2
        new(n, gm)
    end
end

mutable struct Strain
    vertical::Float64
    parallel::Float64
end

mutable struct Material
    name::String
    Fs::Dict{Float64, Float64}
    Fa::Dict{Float64, Float64}
    elasticconstant::Dict{String, Float64}
    latticetype::String
    latticeconstant::Float64
    strainedlatticeconstant::Strain
    vector::BasicLatticeVector
    reciprocalvector::DefaultDict{Float64, Vector{Matrix{Float64}}}
    function Material(name, latticetype)
        if !(latticetype in latticetypes)
            println("Lattice Type is invalid")
        elseif !(name in keys(db.Material))
            println("Invalid material name")
        end
        obj = new(
            name,
            db.Material[name]["pseudopotential"]["Fs"],
            db.Material[name]["pseudopotential"]["Fa"],
            db.Material[name]["elasticconstant"],
            latticetype,
            db.Material[name]["latticeconstant"],
            Strain(
                db.Material[name]["latticeconstant"],
                db.Material[name]["latticeconstant"]
            ),
            BasicLatticeVector(latticetype),
            DefaultDict{Float64, Vector{Matrix{Float64}}}(Vector{Matrix{Float64}})
        )
        reciprocalvector(obj)
        return obj
    end
end

function mix(name::String, material1::Material, alloy1::Float64,  material2::Material, alloy2::Float64)::Material
    if material1.latticetype !== material2.latticetype
        println("Crystal Lattice type is unmutched between given tow materials")
        exit()
    end
    if alloy1 + alloy2 !== 1.0
        println("Total alloy composition should be 1.0")
        exit()
    end
    ternary = Material(material1.name, material1.latticetype)
    ternary.latticeconstant = material1.latticeconstant * alloy1 + material2.latticeconstant * alloy2
    ternary.strainedlatticeconstant.parallel = ternary.latticeconstant
    ternary.strainedlatticeconstant.vertical = ternary.latticeconstant
    for k in keys(ternary.Fs)
        ternary.Fs[k] = material1.Fs[k] * alloy1 + material2.Fs[k] * alloy2
        ternary.Fa[k] = material1.Fa[k] * alloy1 + material2.Fs[k] * alloy2
    end
    for k in keys(ternary.elasticconstant)
        ternary.elasticconstant[k] = material1.elasticconstant[k] * alloy1 + material2.elasticconstant[k] * alloy2
    end
    return ternary
end

function strain(material::Material, by::Float64=material.latticeconstant; percent=1000)
    if percent !== 1000
        if 100.0 < abs(percent)
            println("Strain rate should be in range of -100 ~ 100")
            exit()
        end
        material.strainedlatticeconstant.parallel = material.latticeconstant * (100 + percent)/100
    else
        material.strainedlatticeconstant.parallel = by
    end

    material.strainedlatticeconstant.vertical =
        material.latticeconstant * (
            1 - 2 * ( material.elasticconstant["C12"]/material.elasticconstant["C11"] *
            (material.strainedlatticeconstant.parallel - material.latticeconstant) / material.latticeconstant
            )
        )
end

# Reciprocal lattice vector
function reciprocalvector(material::Material)
    for m in magnitude
        gm = Gm(material.vector, m)
        if gm.magnitude <= 11
            push!(material.reciprocalvector[gm.magnitude], gm.vector)
        end
    end
end

# tau
function tau(material::Material)::Matrix{Float64}
    return [material.strainedlatticeconstant.parallel/8 material.strainedlatticeconstant.parallel/8 material.strainedlatticeconstant.vertical/8]
end

# Ss(Km-Kn)
function Ss(Km::Matrix{Float64}, Kn::Matrix{Float64}, material::Material)
    tmp = (Km - Kn)
    tmp[1, 1] *= 2pi/material.strainedlatticeconstant.parallel
    tmp[1, 2] *= 2pi/material.strainedlatticeconstant.parallel
    tmp[1, 3] *= 2pi/material.strainedlatticeconstant.vertical
    tmp *= tau(material)'
    return cos(tmp)
end

# Sa(Km-Kn)
function Sa(Km::Matrix{Float64}, Kn::Matrix{Float64}, material::Material)
    tmp = (Km - Kn)
    tmp[1, 1] *= 2pi/material.strainedlatticeconstant.parallel
    tmp[1, 2] *= 2pi/material.strainedlatticeconstant.parallel
    tmp[1, 3] *= 2pi/material.strainedlatticeconstant.vertical
    tmp *= tau(material)'
    return sin(tmp)
end

# V(Km-Kn)
function V(Km::Matrix{Float64}, Kn::Matrix{Float64}, material::Material)
    K = Km - Kn
    m = K[1]^2 + K[2]^2 + K[3]^2
    strained = material.latticeconstant^3 / (
        material.strainedlatticeconstant.parallel^2 * material.strainedlatticeconstant.vertical
    )
    if m in [0,3,4,8,11]
        return (Ss(Km, Kn, material) * material.Fs[m] * strained - im * Sa(Km, Kn, material) * material.Fa[m])[1, 1] * strained
    else
        return Complex(0.0)
    end
end

# Hamiltonian matrix
function Hamiltonian(material::Material)::Matrix{ComplexF64}
    N = 51
    hamiltonian = zeros(Complex{Float64}, N, N)
    column = 1
    for i in [0,3,4,8,11]
        for k in material.reciprocalvector[i]
            row  = 1
            for j in [0,3,4,8,11]
                for kk in material.reciprocalvector[j]
                    if column !== row
                        hamiltonian[row, column] =  V(kk, k, material)
                    end
                    row += 1
                end
            end
            column += 1
        end
    end
    return hamiltonian
end

struct BandStructure
    energy::Matrix{Float64}
    kpoints::Vector{Float64}
end

function EigenEnergy(material::Material, kpoints=40)::BandStructure
    sorted = sort(collect(material.reciprocalvector), by=x->x[1])
    diagonal = Vector{Matrix{Float64}}()

    kvector = Vector{Float64}()

    A = [1/material.strainedlatticeconstant.parallel 0 0;
        0 1/material.strainedlatticeconstant.parallel 0;
        0 0 1/material.strainedlatticeconstant.vertical;]

    # Hamiltonian's diagnonal factor
    for K in sorted
        for KK in K[2]
            push!(diagonal, KK * A' * 2pi)
        end
    end

    strained = material.latticeconstant^3 / (
        material.strainedlatticeconstant.parallel^2 * material.strainedlatticeconstant.vertical
    )

    # setup Hamiltonian
    H = Hamiltonian(material)

    E = zeros(Float64, 51, kpoints)

    G = trunc(Int, kpoints/2)
    for k in 1:kpoints
        # L -> G
        if k < G+1
            dk = strained * (pi/material.latticeconstant) / G
            ddk = dk * (k * strained - 1)
            push!(kvector, strained*pi/material.latticeconstant - ddk)
            for i in 1:51
                H[i, i] = (HBAR^2/2M * ( [(pi/material.strainedlatticeconstant.parallel)-ddk (pi/material.strainedlatticeconstant.parallel)-ddk (pi/material.strainedlatticeconstant.vertical)-ddk] + diagonal[i] ) * ( [(pi/material.strainedlatticeconstant.parallel)-ddk (pi/material.strainedlatticeconstant.parallel)-ddk (pi/material.strainedlatticeconstant.vertical)-ddk] + diagonal[i] )')[1,1]
            end
            ε , ψ = LinearAlgebra.eigen(H)
            E[:, k] = sort((ε/Q))
        # G -> X
        else
            dk = strained * (2pi/material.latticeconstant) / (kpoints - G)
            ddk = dk * (k*strained - G - 1)
            push!(kvector, ddk)
            for i in 1:51
                H[i, i] = (HBAR^2/2M * ( [ddk 0 0] + diagonal[i] ) * ( [ddk 0 0] + diagonal[i] )')[1,1]
            end
            ε , ψ = LinearAlgebra.eigen(H)
            E[:, k] = sort((ε/Q))
        end
    end

    e0 = E[4, G]
    f(x) = real(x - e0)
    E = f.(E)

    return BandStructure(E, kvector)
end

function Plot(x::Vector{Float64}, y::Vector{Float64}, title="Strained Electron Effective Mass";show=false)
    try
        if show == false
            plot(x, y, title=title, legend=false)
            tmp = split(title)
            filename = join(tmp, "_")
            savefig("$(filename).png")
        else
            plot(x, y, title=title, legend=false)
        end
    catch e
        println(e)
    end
end

function BandPlot(band::BandStructure; title="Band Structure", show=false)
    try
        plot([iknum for iknum in 1:length(band.energy[1, :])], band.energy[1, :], title=title, legend=false)
        for ii in 2:9
            plot!([iknum for iknum in 1:length(band.energy[ii, :])], band.energy[ii, :], title=title, legend=false)
        end
        if show == false
            tmp = split(title)
            filename = join(tmp, "_")
            savefig("$filename.png")
        end
    catch e
        println(e)
    end
end

function BandGap(band::BandStructure)
    valenceband = band.energy[4, :]
    conductionband = band.energy[5, :]
    G = trunc(Int, length(conductionband)/2)
    gap = conductionband[G] - valenceband[G]
    println("Band Gap : ", gap, " [eV]")
    return gap
end

# y = ax^2 + bx = c
function quadraticInterpolation(x::Vector{Float64}, y::Vector{Float64})
    quartic(x) = x^4
    cubic(x) = x^3
    quadratic(x) = x^2
    X =
    [
        sum(quartic.(x)) sum(cubic.(x)) sum(quadratic.(x));
        sum(cubic.(x)) sum(quadratic.(x)) sum(x);
        sum(quadratic.(x)) sum(x) length(x);
    ]
    Y =
    [
        sum( (quadratic.(x))' * y );
        sum(x'*y);
        sum(y);
    ]
    a,b,c = X' * Y
    return a,b,c
end

# return electron effective mass
function EffectiveMass(band::BandStructure, from, to)
    conductionband = band.energy[5, :]
    a,b,c = quadraticInterpolation(band.kpoints[from:to], conductionband[from:to])
    emass = M^2 * a
    println("Effective Mass : ", emass)
    return emass
end

end