text:
using ElectronicStructure
using ElectronicStructure.PsPot

function test_main()
  #psp1 = PsPot_HGH(1, "H" , "LDA_HGH/1h.1.hgh")
  #psp2 = PsPot_HGH(2, "Ti", "LDA_HGH/22ti.12.hgh")
  #psp = PsPot_HGH(3, "Sm", "LDA_HGH/62sm.16.hgh")
  #psp4 = PsPot_HGH(4, "Pt", "LDA_HGH/78pt.18.hgh")
  psp = PsPot_HGH(1, "Ge", "../LDA_HGH/32ge.4.hgh")

  info_PsPot_HGH(psp)

end

test_main()


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

using Markdown
using InteractiveUtils

# ╔═╡ 61dc86f0-238f-11eb-054f-41a45b78e94d
begin
	using  MLJ
	import RDatasets: dataset
    import DataFrames: describe, select, Not, rename!,head
	import Plots: plot,plot!,scatter,scatter!,histogram,histogram!
end

# ╔═╡ ff5bd480-29c0-11eb-3112-cb35c17acc60
md" ### Libs"

# ╔═╡ 8e226310-29c0-11eb-1093-f95e53c029b6
md"### Carga de Dados"

# ╔═╡ 114c91a0-2a90-11eb-365b-4d360bb51ebe
md" É interessante sempre mostrarmos os primeiros exemplos dos dados para ter uma ideia dos valores. Entretanto, como visto no curso, há várias análises possíveis para se fazer como identificar nulos, outliers, ... Tudo isto, é bom ser tratado."

# ╔═╡ b6faf1b0-294a-11eb-38ae-0924c4943d81
begin
	boston = dataset("MASS", "Boston")
	head(boston)
end

# ╔═╡ a2efe1c0-29c8-11eb-287d-49b90e597876
md"### Sobre o Problema

url: https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html
"

# ╔═╡ a04c5f12-29d3-11eb-12a6-3f1b8acc88dc
md"#### Atributos"

# ╔═╡ a417e700-29c8-11eb-3ee5-6d2c302d95ac
md"There are 14 attributes in each case of the dataset. They are:

CRIM - per capita crime rate by town

ZN - proportion of residential land zoned for lots over 25,000 sq.ft.

INDUS - proportion of non-retail business acres per town.

CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise)

NOX - nitric oxides concentration (parts per 10 million)

RM - average number of rooms per dwelling

AGE - proportion of owner-occupied units built prior to 1940

DIS - weighted distances to five Boston employment centres

RAD - index of accessibility to radial highways

TAX - full-value property-tax rate per $10,000

PTRATIO - pupil-teacher ratio by town

B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town

LSTAT - % lower status of the population

MEDV - Median value of owner-occupied homes in $1000's"

# ╔═╡ 9a810b20-29c0-11eb-1ca3-d155108c6e52
md"### Estatísticas dos Dados"

# ╔═╡ 91187080-2a8f-11eb-3001-1973d885e3f0
md"É sempre bom tiramos umas estatísticas dos dados para entender melhor o comportamento das variáveis do nosso problema. Qual a variância de cada variável? Elas estão na mesma ordem de grandeza? Máixmo e minímo ?"

# ╔═╡ 2033e842-294f-11eb-2b38-1969d7440b88
describe(boston, :mean, :std, :min, :max, :eltype)

# ╔═╡ 2f06c4a0-294f-11eb-2cb1-2da452f98b0f
data = coerce(boston, autotype(boston, :discrete_to_continuous));

# ╔═╡ d2a9dfe0-29c0-11eb-36a9-e7e3c1d17717
md"### Selecionando Dados"

# ╔═╡ 49b99b50-2a90-11eb-2b54-2704b914334d
md"É nesta fase que você vai escolher seu target e vai selecionar os atributos para trabalhar."

# ╔═╡ 47a80140-294f-11eb-2880-ef7cfb0a8997
begin
	y = data.MedV;
    X = select(data, Not(:MedV));
end

# ╔═╡ e3d83ea0-29a3-11eb-10f5-b5f94102e1f1
md"### Verificando Algoritmos Disponíveis"

# ╔═╡ 60083560-2a90-11eb-373f-3b8a1ba538a0
md"Essa é uma característica bem legal do MLJ. Como nosso dados possuem atributos numéricos e o target numérico, ele me mostra todos os algoritmos que conseguem rodar para este tipo de dado."

# ╔═╡ 868b73b0-29a9-11eb-0ef8-971446a3c919
models(matching(X, y))

# ╔═╡ 87cedb30-2a90-11eb-26dd-6162d49edd88
md"Abaixo realizamos uma carga no nosso regressor linear."

# ╔═╡ 6b75f5a0-294f-11eb-0835-cfe4a412cffe
begin
	@load LinearRegressor pkg=MLJLinearModels
	mdl = LinearRegressor();
end

# ╔═╡ 9abff260-2a90-11eb-026f-59a85ec6b64c
md"O MLJ nos permite montar um pipeline. Pipeline é uma forma de montar uma sequência de operações poupando-lhe código. No nosso caso, queremos padronizar X e y antes de realizarmos o treino do regressor. Podemos fazer isso em uma única linha."

# ╔═╡ fcde2920-29a6-11eb-37f7-4befad3ee3d2
predictor = @pipeline Standardizer mdl target=Standardizer

# ╔═╡ 323d2620-29c0-11eb-258a-4703b4f5a0d2
md"### Partição dos Dados"

# ╔═╡ e57688a2-2a90-11eb-2272-3d9898bf4375
md"Realizamos aqui o que conhecemos como hold-out. Particionamos em treino e teste, que é a forma mais simples de se fazer uma validação. Sabemos que existem outras no curso. No caso, deixei 70% para treino e 30% para teste. Como podem notar, a função gera índices qu correspondem às linhas dos exemplos. Ele não tira cópia dos dados."

# ╔═╡ fca39122-29a6-11eb-3b25-c7f0df7fc68f
train, test  = partition(eachindex(y), 0.7, shuffle=true)

# ╔═╡ 287ae7e0-2a91-11eb-0d5b-b138443bdc29
md"Um conceito interessante do MLJ é o machine. Esta função faz uma mapeamento do preditor ao dado. Quando aplicamos machine, estamos gerando um nó num grafo de computação e um estado associado. Existem algmas vantagens neste paradigma como: associar um mesmo preditor a vários dados distintos que podem ser processados em diferentes recursos; evitar re execução de um mesmo preditor e um mesmo dado."

# ╔═╡ fc474050-29a6-11eb-1fbd-87a88d38525a
pred_machine = machine(predictor, X, y);

# ╔═╡ 4449f8d0-29b5-11eb-3f83-c352d034ecb2
md"### Treino da regressão linear "

# ╔═╡ 91918630-2a91-11eb-1113-43b18086f208
md" Aqui temos enfim, o treino do nosso regressor para o data X, apenas no conjunto de treino. Isto porque estamos passando os índices do exemplos."

# ╔═╡ 0de5e6e0-29a7-11eb-0558-2bb10a5fe1a5
fit!(pred_machine, rows=train);

# ╔═╡ 50f36580-29b5-11eb-06f2-cf4b9e142979
md"### Análise da Predição"

# ╔═╡ ef445730-2a91-11eb-1552-c3e7d6e5ff34
md"Aqui analisamos cada predição e cada target. Este tipo de análise é mais apropriada quando nosso conjunto de dados possui uma relação de ordem temporal entre os exemplos (que não é o caso). Entretanto, é uma análise boa a se fazer. Podemos notar que há um atendência de subestimar os valores mais altos do que valores mais baixos. Você pode analisar o que pode estar gerando dificldade nisto analisando alguns exemplos individualmente."

# ╔═╡ 152f7100-29a7-11eb-2300-735ca72b15c3
ŷ = predict(pred_machine, rows=test)

# ╔═╡ c29c1380-29b5-11eb-254c-2d026e15f903
begin
	plot(ŷ,label="Target")
	scatter!(y[test],label="Label")

end

# ╔═╡ 89d06380-29bf-11eb-0d5b-4533c94745d5
md" Distribuição da Predição"

# ╔═╡ 8b3ac8f0-29bf-11eb-274e-f7b07eb35e7c
begin
   histogram(ŷ,label="Target")
end

# ╔═╡ 1e35c150-29de-11eb-0069-cb75699b6d61
histogram(y[test],label="Label")

# ╔═╡ 76b4c570-29b5-11eb-32d5-75da0001a894
md"### Parâmetros encontrados"

# ╔═╡ 46a9fc40-2a8e-11eb-301f-077f9a6d50fe
md" Muitas vezes é interessante analisarmos a relevância dos coeficientes por váriosmotivos. Isto porque possibilita descobrir quais são as variáveis que ajudam mais num problema, ajudando inclusive na tomada de decisão de um ser humano. Adicionalmente, pode ser usado para compactar o modelo (ambientes com poucos recursos), eliminando atributos pouco relevantes. Finalmente, pode-se utilizar para verificar se seu modelo usa alguma variável que seja discrimnatória em termos de sociedade (ex. idade para tomada de crédito). "

# ╔═╡ 8d50cc90-294f-11eb-012b-3b973cbd1f1e
begin
	fp = fitted_params(pred_machine)
	coefs = fp.fitted_params_given_machine[fp.machines[2]]
	#fp.machines
end

# ╔═╡ 2d7233c0-29c4-11eb-0e0f-2f9b5773a062
begin
	names = [string(e[1]) for e in coefs[1] ]
	vals  = [e[2] for e in coefs[1] ]
	plot(names,vals,linetype = :bar,label="Coefs")
end

# ╔═╡ dd923000-2a8e-11eb-1100-11a1ff5510bc
md"Coeficientes positivos significam que contribuem positivamente para o incremento de y. Os negativos, de forma oposta. Já os com valores próximos de zero, são pouco relevantes ao probelma."

# ╔═╡ eaad8980-29ab-11eb-059a-4fe1c6b0bb2d
md"### Análise de Erro"

# ╔═╡ 7cc2ab20-2a8d-11eb-2a73-e384f6d3b477
md" Separei aqui duas métricas que já vimos no curso: mae e rmse. Com o mae, conseguimos ter uma ideia do quanto erramos para cima ou para baixo com relação ao $y$. Erramos $3.2$ em termo de dinheiro para cima e para baixo do gabarito. Já o rmse indica uma variabilidade de $4.4$. Note que podemos ter um mae baixo, porém, se tivermos muito outliers, o rmse aumenta."

# ╔═╡ 2dadd5c0-2a8b-11eb-077d-d174e5c3e4dc
md"$\frac{1}{n}\sum_{i=1}^{n}{|ŷ - y|}$"

# ╔═╡ 1b9e5290-29a7-11eb-1c8f-4f12216b4391
md"mae = $(round(mae(ŷ, y[test]), sigdigits=4))"

# ╔═╡ a0b1a6f0-2a8b-11eb-1c89-d534b1671eec
md"$\sqrt{\frac{1}{n}\sum_{i=1}^{n}{(ŷ - y)^2}}$"

# ╔═╡ c3fb2bc0-295c-11eb-0061-c3875f5cb01b
md"rmse = $(round(rms(ŷ, y[test]), sigdigits=4))"

# ╔═╡ 0ebfd9e0-29ac-11eb-055f-615e53896367
md" Abaixo uma análise muito utilizada quando há uma relação de ordem temporal entre os exemplos. Apesar de neste problema não termos uma ordem temporal, deixo aqui este tipo de análise.Neste gráfico podemos notar que existem certos exemplos muito difíceis que talvez seja interessante analisar no futuro olhando seus atributos."

# ╔═╡ e102cfc0-295c-11eb-2c5b-834000ca2077
begin
	res = ŷ .- y[test]
	plot(res)
	scatter!(res)
end

# ╔═╡ 2b478eb0-2a8c-11eb-0fe4-17dad212a408
md"Também é interessante analisarmos a distribuição do erro a partir de um histograma. Desta forma, temos uma noção maior da magnitude do erro e sinal. Da mesma forma que na análise anterior, vemos que existem exemplos que são bem mais difíceis que o normal. Cabe uma análise individual dos exemplos."

# ╔═╡ 845c9050-29aa-11eb-3587-1ff2b98c8905
histogram(res,label="Erro")

# ╔═╡ Cell order:
# ╟─ff5bd480-29c0-11eb-3112-cb35c17acc60
# ╠═61dc86f0-238f-11eb-054f-41a45b78e94d
# ╟─8e226310-29c0-11eb-1093-f95e53c029b6
# ╟─114c91a0-2a90-11eb-365b-4d360bb51ebe
# ╠═b6faf1b0-294a-11eb-38ae-0924c4943d81
# ╟─a2efe1c0-29c8-11eb-287d-49b90e597876
# ╟─a04c5f12-29d3-11eb-12a6-3f1b8acc88dc
# ╟─a417e700-29c8-11eb-3ee5-6d2c302d95ac
# ╟─9a810b20-29c0-11eb-1ca3-d155108c6e52
# ╟─91187080-2a8f-11eb-3001-1973d885e3f0
# ╠═2033e842-294f-11eb-2b38-1969d7440b88
# ╠═2f06c4a0-294f-11eb-2cb1-2da452f98b0f
# ╟─d2a9dfe0-29c0-11eb-36a9-e7e3c1d17717
# ╟─49b99b50-2a90-11eb-2b54-2704b914334d
# ╠═47a80140-294f-11eb-2880-ef7cfb0a8997
# ╟─e3d83ea0-29a3-11eb-10f5-b5f94102e1f1
# ╟─60083560-2a90-11eb-373f-3b8a1ba538a0
# ╠═868b73b0-29a9-11eb-0ef8-971446a3c919
# ╟─87cedb30-2a90-11eb-26dd-6162d49edd88
# ╠═6b75f5a0-294f-11eb-0835-cfe4a412cffe
# ╟─9abff260-2a90-11eb-026f-59a85ec6b64c
# ╠═fcde2920-29a6-11eb-37f7-4befad3ee3d2
# ╟─323d2620-29c0-11eb-258a-4703b4f5a0d2
# ╟─e57688a2-2a90-11eb-2272-3d9898bf4375
# ╠═fca39122-29a6-11eb-3b25-c7f0df7fc68f
# ╟─287ae7e0-2a91-11eb-0d5b-b138443bdc29
# ╠═fc474050-29a6-11eb-1fbd-87a88d38525a
# ╟─4449f8d0-29b5-11eb-3f83-c352d034ecb2
# ╟─91918630-2a91-11eb-1113-43b18086f208
# ╠═0de5e6e0-29a7-11eb-0558-2bb10a5fe1a5
# ╟─50f36580-29b5-11eb-06f2-cf4b9e142979
# ╟─ef445730-2a91-11eb-1552-c3e7d6e5ff34
# ╠═152f7100-29a7-11eb-2300-735ca72b15c3
# ╠═c29c1380-29b5-11eb-254c-2d026e15f903
# ╠═89d06380-29bf-11eb-0d5b-4533c94745d5
# ╠═8b3ac8f0-29bf-11eb-274e-f7b07eb35e7c
# ╠═1e35c150-29de-11eb-0069-cb75699b6d61
# ╟─76b4c570-29b5-11eb-32d5-75da0001a894
# ╟─46a9fc40-2a8e-11eb-301f-077f9a6d50fe
# ╟─8d50cc90-294f-11eb-012b-3b973cbd1f1e
# ╠═2d7233c0-29c4-11eb-0e0f-2f9b5773a062
# ╟─dd923000-2a8e-11eb-1100-11a1ff5510bc
# ╟─eaad8980-29ab-11eb-059a-4fe1c6b0bb2d
# ╟─7cc2ab20-2a8d-11eb-2a73-e384f6d3b477
# ╟─2dadd5c0-2a8b-11eb-077d-d174e5c3e4dc
# ╟─1b9e5290-29a7-11eb-1c8f-4f12216b4391
# ╟─a0b1a6f0-2a8b-11eb-1c89-d534b1671eec
# ╟─c3fb2bc0-295c-11eb-0061-c3875f5cb01b
# ╟─0ebfd9e0-29ac-11eb-055f-615e53896367
# ╟─e102cfc0-295c-11eb-2c5b-834000ca2077
# ╟─2b478eb0-2a8c-11eb-0fe4-17dad212a408
# ╟─845c9050-29aa-11eb-3587-1ff2b98c8905

ext: jl

max_stars_repo_name: lalvim/machinelearning

max_issues_repo_name: lalvim/machinelearning

max_forks_repo_name: lalvim/machinelearning


text:
# setup
using ZipFile

function get_coastline_shapefile()
    dir = mktempdir()
    zip = Base.download("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_coastline.zip")

    # read the Zip file
    r = ZipFile.Reader(zip);

    # populate the folder
    cd(dir) do
        for f in r.files
            open(f.name, "w") do file
                write(file, read(f, String))
            end
        end
    end

    return dir
end

using Proj4, Shapefile, GeoInterface, Makie, MakieLayout, PlotUtils

# path to shapefile of coastline
path = "ne_110m_coastline/ne_110m_coastline.shp"

# plain projection
src_proj = Projection("+proj=longlat +datum=WGS84 +no_defs")

# Robinson projection
dest_proj = Projection("+proj=robin")

# load the shapefile
table = Shapefile.Table(path)

# remove any missing elements (of which there are none)
shapes = collect(skipmissing(Shapefile.shapes(table)))

# index into all arrays once - there is one extra layer of arrays
data = getindex.(GeoInterface.coordinates.(shapes), 1)

# compute the length of the eventual padded array
length_of_data = sum(length.(data)) + length(data)

# preallocate an array
linevec_transformed = Vector{Point2f0}(undef, length_of_data)
linevec_original = Vector{Point2f0}(undef, length_of_data)

function to_nansep_vec(f::Function, data::Vector{Vector{T}}) where T

    length_of_data = sum(length.(data)) + length(data)

    lvec = Vector{T}(undef, length_of_data)

    pos = 1

    for (i, datum) in enumerate(data)
        lvec[pos:(pos+length(datum)-1)] .= f(datum)
        pos += length(datum)
        lvec[pos] = Point2f0(NaN, NaN)
        pos += 1
    end

    return lvec

end

to_nansep_vec(data::Vector{Vector{T}}) where T = to_nansep_vec(identity, data)

linevec_transformed = to_nansep_vec(data) do datum
    transform.(src_proj, dest_proj, datum)
end

plv = Point2f0.(linevec_transformed)

linevec_original = to_nansep_vec(identity, data) .|> Point2f0

# loop through all linestrings, transform, and merge
# (you can also skip the merging, and make individual lineplots)
# for (i, datum) in enumerate(data)
#     global pos
#     linevec_transformed[pos:(pos+length(datum)-1)] .= transform.(src_proj, dest_proj, datum)
#     linevec_original[pos:(pos+length(datum)-1)]    .= datum
#
#     pos += length(datum)
#     linevec_transformed[pos] = Point2f0(NaN, NaN)
#     linevec_original[pos]    = Point2f0(NaN, NaN)
#     pos += 1
# end

# rudimentary comparison plot
scene, layout = layoutscene()

sc1 = layout[1, 1] = LAxis(scene, title = "Lat/Long", xticks = AutoLinearTicks(100f0), yticks = AutoLinearTicks(100f0), aspect = AxisAspect(1))
sc2 = layout[1, 2] = LAxis(scene, title = "Robinson", xticks = AutoLinearTicks(500f0), yticks = AutoLinearTicks(500f0), aspect = AxisAspect(1))

lines!(sc1, linevec_original;    color = Makie.Colors.colorant"#add8e6", linewidth = 2)
lines!(sc2, linevec_transformed; color = Makie.Colors.colorant"#add8e6", linewidth = 2)

display(scene)

"""
    axis_points_for_geo(xmin, xmax, ymin, ymax)

Returns a vector of NaN separated Point2 series which represent the sampled axis
for transformation.
"""
function axis_grid_for_geo(xmin, xmax, ymin, ymax; density = 100, x_k_min=5, x_k_max=10, x_k_ideal=7, y_k_min=5, y_k_max=10, y_k_ideal=7)

    xticks = PlotUtils.optimize_ticks(xmin, xmax; k_min = x_k_min, k_max = x_k_max, k_ideal = x_k_ideal, extend_ticks = true)[1]
    yticks = PlotUtils.optimize_ticks(ymin, ymax; k_min = y_k_min, k_max = y_k_max, k_ideal = y_k_ideal, extend_ticks = true)[1]

    numticks = length(xticks) + length(yticks)

    linevec = Vector{Point2f0}()

    for xtick in [xticks]
        append!(linevec, transform.(src_proj, dest_proj, Point2f0.(xtick, LinRange(ymin, ymax, density))))
        push!(linevec, Point2f0(NaN))
    end
    for ytick in yticks
        append!(linevec, transform.(src_proj, dest_proj, Point2f0.(LinRange(xmin, xmax, density), ytick)))
        push!(linevec, Point2f0(NaN))
    end

    return linevec

end

"""
    axis_points_for_geo(xmin, xmax, ymin, ymax)

Returns a vector of NaN separated Point2 series which represent the sampled axis
frame which has been transformed..
"""
function axis_frame_for_geo(xmin, xmax, ymin, ymax; density = 100)
    linevec = Vector{Point2f0}(undef, (density+1)*4)

    append!(linevec, transform.(src_proj, dest_proj, Point2f0.(xmin, LinRange(ymin, ymax, density))))
    push!(linevec, Point2f0(NaN))

    append!(linevec, transform.(src_proj, dest_proj, Point2f0.(xmax, LinRange(ymin, ymax, density))))
    push!(linevec, Point2f0(NaN))

    append!(linevec, transform.(src_proj, dest_proj, Point2f0.(LinRange(xmin, xmax, density), ymin)))
    push!(linevec, Point2f0(NaN))

    append!(linevec, transform.(src_proj, dest_proj, Point2f0.(LinRange(xmin, xmax, density), ymax)))
    push!(linevec, Point2f0(NaN))

    return linevec

end

xmin, xmax = (-180, 180)
ymin, ymax = (-90, 90)

sc = lines(transform.(src_proj, dest_proj, linevec_original); color = Makie.Colors.colorant"#add8e6", linewidth = 2, show_axis = false)
lines!(sc, axis_points_for_geo(xmin, xmax, ymin, ymax); color = RGBA(Makie.Colors.colorant"light grey", 0.7))
lines!(sc, axis_frame_for_geo(xmin, xmax, ymin, ymax))

display(sc) # projected thingamajig

function draw_geo_axis!(
        scene,
        data,
        attributes::Attributes;
        source = crs(data) === nothing ? Projection("+proj=longlat +datum=LonLat +no_defs") : crs(data),
        dest = Projection("+proj=longlat +datum=LonLat +no_defs")
    )

    limits = boundingbox(data)

end


text:
#include("constants.jl")

# schemes
cpl=1;
moist=1;

#boundary conditions
TAup=270; #top model temperature
UAup=(0.2/karman).*(log((ZAF[end]+ΔZA[end]/2)/1e-4)); #top model temperature
UOdown=0; #lower boundary model current

# numerics
Aimp=0.5; # Atmospheric implicitness
Oimp=0.5; # Oceanic implicitness

@with_kw struct snparameters{R}
  #soil numerical parameters and domain
  Zs::R = 2.; #soil's depth [m]
  Ns::Int = 201; #number of cells in soil colomn
  dz_s::R = Zs/(Ns-1); #soil cell size [m]
  z_s = Zs.-((1:Ns).-1).*dz_s; #soil depth at each cell [m]
  eps_h::R = 1e-5; # matric head tolerance
  eps_T::R = 1e-5; #soil temperature tolerance
end

soil_numerical_parameters=snparameters{Float64}()


text:
@testset "Maximal Independent Set" begin

    g0 = SimpleGraph(0)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g0)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            @test isempty(z)
        end
    end

    g1 = SimpleGraph(1)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g1)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            @test (z == [1,])
        end
    end

    add_edge!(g1, 1, 1)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g1)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            isempty(z)
        end
    end

    g3 = star_graph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g3)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            @test (length(z)== 1 || length(z)== 4)
        end
    end

    g4 = complete_graph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g4)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            @test length(z)== 1 #Exactly one vertex
        end
    end

    g5 = path_graph(5)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g5)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            sort!(z)
            @test (z == [2, 4] || z == [2, 5] || z == [1, 3, 5] || z == [1, 4])
        end
    end

    add_edge!(g5, 2, 2)
    add_edge!(g5, 3, 3)
    for parallel in [:threads, :distributed]
        for g in testgraphs(g5)
            z = @inferred(Parallel.independent_set(g, 4, MaximalIndependentSet(); parallel=parallel))
            sort!(z)
            @test (z == [4,] || z == [5,] || z == [1, 5] || z == [1, 4])
        end
    end
end


text:
using .Plots
using .StatsPlots
function plot_im(mean_img)
    p = plot()
    heatmap!(p,
              imagepixels(mean_img)...,
              mean_img,
              aspect_ratio=:equal,
              size=(500,400),
              xflip=true
            )
    xlims!(p,-80.0, 80.0)
    ylims!(p,-80.0, 80.0)
    title!(p,"Mean Image")
    return p
end

function plot_samples(sims)
    anim = @animate for s in sims[1:end]
        heatmap(imagepixels(s)..., s, aspect_ratio=:equal, size=(500,400), xflip=true)
        xlims!(-80.0,80.0)
        ylims!(-80.0, 80.0)
    end
    return anim
end

function plot_vis_comp(visobs, gamps, gphase, m)
    u = Comrade.getdata(visobs, :u)
    v = Comrade.getdata(visobs, :v)
    visr = Comrade.getdata(visobs, :visr)
    visi = Comrade.getdata(visobs, :visi)
    error = Comrade.getdata(visobs, :error)
    bl = Comrade.getdata(visobs, :baseline)

    mvis = Comrade.visibility.(Ref(m), u, v)
    for i in eachindex(u)
        s1,s2 = bl[i]
        ga1 = gamps[s1]
        ga2 = gamps[s2]
        gθ1 = gphase[s1]
        gθ2 = gphase[s2]
        s,c = sincos(gθ1 - gθ2)
        vr = ga1*ga2*(real(mvis[i])*c + imag(mvis[i])*s)
        vi = ga1*ga2*(-real(mvis[i])*s + imag(mvis[i])*c)
        mvis[i] = vr + vi*1im
    end

    p = plot()

    scatter!(p,hypot.(u,v)*rad2μas/1e9, visr,
                  yerr=error,
                  color=:cornflowerblue, label="Real Data", markershape=:square, alpha=0.5)
    scatter!(p,hypot.(u,v)*rad2μas/1e9, visi,
                  yerr=error,
                  color=:orange, label="Imag Data", markershape=:square, alpha=0.5)
    scatter!(p,hypot.(u,v)*rad2μas/1e9, real.(mvis),
                  color=:blue, label="Real Model")
    scatter!(p,hypot.(u,v)*rad2μas/1e9, imag.(mvis),
                  color=:red, label="Imag Model")
    ylabel!("Visibility")
    xlabel!("uvdist Gλ")
    return p

end

function plot_amp_comp(ampobs, gains, m)
    u = Comrade.getdata(ampobs, :u)
    v = Comrade.getdata(ampobs, :v)
    amp = Comrade.getdata(ampobs, :amp)
    error = Comrade.getdata(ampobs, :error)
    bl = Comrade.getdata(ampobs, :baseline)

    vamp = Comrade.amplitude.(Ref(m), u, v)
    for i in eachindex(u)
        s1,s2 = bl[i]
        g1 = gains[s1]
        g2 = gains[s2]
        vamp[i] =  g1*g2*vamp[i]
    end

    p = plot()

    scatter!(p,hypot.(u,v)*rad2μas/1e9, amp,
                  yerr=error,
                  color=:cornflowerblue, label="Data", markershape=:square, alpha=0.5)
    scatter!(p,hypot.(u,v)*rad2μas/1e9, vamp,
                  color=:blue, label="Model")
    ylabel!("V Amplitude")
    xlabel!("uvdist Gλ")
    return p

end

function plot_cp_comp(obscp, m)
    u1 = Comrade.getdata(obscp, :u1)
    v1 = Comrade.getdata(obscp, :v1)
    u2 = Comrade.getdata(obscp, :u2)
    v2 = Comrade.getdata(obscp, :v2)
    u3 = Comrade.getdata(obscp, :u3)
    v3 = Comrade.getdata(obscp, :v3)
    cp = Comrade.getdata(obscp, :phase)
    error = Comrade.getdata(obscp, :error)

    mcp = Comrade.closure_phase.(Ref(m),
                               u1, v1,
                               u2, v2,
                               u3, v3
                            )

    p = plot()

    tridist = @. sqrt(u1^2 + v1^2 + u2^2 + v2^2 + u3^2 + v3^2)*rad2μas/1e9
    scatter!(p, tridist, cp,
                  yerr=error,
                  color=:cornflowerblue, label="Data", markershape=:square, alpha=0.5)
    scatter!(p,tridist, mcp,
                  color=:blue, label="Model")
    ylabel!("Closure Phase")
    xlabel!("tridist Gλ")
    return p

end

function plot_res_density(anres, cpnres)
    p = plot()
    density!(p,anres, label="Amplitude")
    density!(p,cpnres, label="CP")
    plot!(p,x->pdf(Normal(), x), label="Std. Normal")
    xlabel!(p,"Normalized Residuals")
    return p
end


text:
# # Fractional Maxwell
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_maxwell.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_maxwell.ipynb)

using RHEOS
## include a helper function for plotting
include("assets/plothelper.jl");
#-

Fract_Maxwell

# #### Constitutive Equation
# ```math
# \sigma(t) + \frac{c_\alpha}{c_\beta} \frac{d^{\alpha-\beta} \sigma(t)}{dt^{\alpha-\beta}}= c_{\alpha} \frac{d^\alpha \epsilon(t)}{dt^\alpha}
# ```

# ```math
# \text{for}\; \ 0 \leq \beta \leq \alpha \leq 1
# ```

# #### Relaxation Modulus
# ```math
# G(t) = c_\beta t^{-\beta} E_{\alpha-\beta,1-\beta}\left(-\frac{c_\beta}{c_\alpha} t^{\alpha-\beta}\right)
# ```

# #### Creep Modulus
# ```math
# J(t) = \frac{1}{c_\alpha \Gamma(1+\alpha)}t^\alpha+\frac{1}{c_\beta \Gamma(1+\beta)}t^\beta
# ```

# #### Storage Modulus
# ```math
# G^{\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \cos(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \cos(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})}
# ```

# #### Loss Modulus
# ```math
# G^{\prime\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \sin(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \sin(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})}
# ```

# # Fractional (Spring) Maxwell

FractS_Maxwell
#-

models = Vector{RheoModel}()

## plot moduli for varying α
for alpha in [0.1, 0.3, 0.5, 0.7, 0.9]

    push!(models, RheoModel(FractS_Maxwell, (cₐ = 1.0, a = alpha, k = 1.0)))

end

plotmodel(models)

# # Fraction (Dashpot) Maxwell

FractD_Maxwell
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.3, 0.5, 0.7, 0.9]

    push!(models, RheoModel(FractD_Maxwell, (η = 10, cᵦ= 1.0, β = beta)))

end

plotmodel(models, ymaxG = 2.0)

# # Maxwell

Maxwell
#-

models = Vector{RheoModel}()

## plot moduli for varying k
for k in [5.0, 10.0, 20.0]

    push!(models, RheoModel(Maxwell, (η = 10, k = k)))

end

plotmodel(models)


text:
export
    left_env,
    right_env,
    dot!

# --------------------------- Conventions -----------------------
#
#      MPS          MPS*         MPO       left env     right env
#       2            2            2           - 1          2 -
#   1 - A - 3    1 - B - 3    1 - W - 3      L               R
#                                 4           - 2          1 -
# ---------------------------------------------------------------
#

function LinearAlgebra.dot(ϕ::AbstractMPS, ψ::AbstractMPS)
    T = promote_type(eltype(ψ), eltype(ϕ))
    C = ones(T, 1, 1)
    for (A, B) ∈ zip(ψ, ϕ)
        @tensor C[x, y] := conj(B)[β, σ, x] * C[β, α] * A[α, σ, y] order = (α, β, σ)
    end
    tr(C)
end

"""
Creates left environment (ϕ - bra, ψ - ket)
"""
function left_env(ϕ::AbstractMPS, ψ::AbstractMPS)
    T = promote_type(eltype(ψ), eltype(ϕ))
    S = typeof(similar(ψ[1], T, (1, 1)))
    L = Vector{S}(undef, length(ψ)+1)
    L[1] = similar(ψ[1], T, (1, 1))
    L[1][1, 1] = one(T)
    for (i, (A, B)) ∈ enumerate(zip(ψ, ϕ))
        C = L[i]
        @tensor C[x, y] := conj(B)[β, σ, x] * C[β, α] * A[α, σ, y] order = (α, β, σ)
        L[i+1] = C
    end
    L
end

# TODO: remove it (after SpinGlassEngine is updated)
@memoize Dict function left_env(ϕ::AbstractMPS, σ::Vector{Int})
    l = length(σ)
    if l == 0 return ones(eltype(ϕ), 1) end
    m = σ[l]
    L̃ = left_env(ϕ, σ[1:l-1])
    M = ϕ[l]
    @matmul L[x] := sum(α) L̃[α] * M[α, $m, x]
    L
end

"""
Creates right environment (ϕ - bra, ψ - ket)
"""
function right_env(ϕ::AbstractMPS, ψ::AbstractMPS)
    L = length(ψ)
    T = promote_type(eltype(ψ), eltype(ϕ))
    S = typeof(similar(ψ[1], T, (1, 1)))
    R = Vector{S}(undef, L+1)
    R[end] = similar(ψ[1], T, (1, 1))
    R[end][1, 1] = one(T)
    for i ∈ L:-1:1
        M = ψ[i]
        M̃ = ϕ[i]
        D = R[i+1]
        @tensor D[x, y] := M[x, σ, α] * D[α, β] * conj(M̃)[y, σ, β] order = (β, α, σ)
        R[i] = D
    end
    R
end

# TODO: remove it (after SpinGlassEngine is updated)
@memoize Dict function right_env(ϕ::AbstractMPS{T}, W::AbstractMPO{T}, σ::Union{Vector, NTuple}) where {T}
    l = length(σ)
    if l == 0
        R = similar(ϕ[1], T, (1, 1))
        R[1, 1] = one(T)
        return R
    end
    k = length(W)
    R̃ = right_env(ϕ, W, σ[2:l])
    M = ϕ[k-l+1]
    M̃ = W[k-l+1]
    K = @view M̃[:, σ[1], :, :]
    @tensor R[x, y] := K[y, β, γ] * M[x, γ, α] * R̃[α, β] order = (β, γ, α)
    R
end

"""
$(TYPEDSIGNATURES)

Calculates the norm of an MPS \$\\ket{\\phi}\$
"""
LinearAlgebra.norm(ψ::AbstractMPS) = sqrt(abs(dot(ψ, ψ)))

"""
$(TYPEDSIGNATURES)

Calculates \$\\bra{\\phi} O \\ket{\\psi}\$

# Details

Calculates the matrix element of \$O\$
```math
\\bra{\\phi} O \\ket{\\psi}
```
in one pass, utlizing `TensorOperations`.
"""
function LinearAlgebra.dot(ϕ::AbstractMPS, O::Union{Vector, NTuple}, ψ::AbstractMPS) #where T <: AbstractMatrix
    S = promote_type(eltype(ψ), eltype(ϕ), eltype(O[1]))
    C = similar(ψ[1], S, (1, 1))
    C[1, 1] = one(S)
    for (A, W, B) ∈ zip(ϕ, O, ψ)
        @tensor C[x, y] := conj(A)[β, σ, x] * W[σ, η] * C[β, α] * B[α, η, y] order = (α, η, β, σ)
    end
    tr(C)
end

function LinearAlgebra.dot(O::AbstractMPO, ψ::AbstractMPS)
    S = promote_type(eltype(ψ), eltype(O))
    T = typeof(ψ)
    ϕ = T.name.wrapper(S, length(ψ))
    for (i, (A, B)) ∈ enumerate(zip(O, ψ))
        @matmul N[(x, a), σ, (y, b)] := sum(η) A[x, σ, y, η] * B[a, η, b]
        ϕ[i] = N
    end
    ϕ
end

function LinearAlgebra.dot(O1::AbstractMPO, O2::AbstractMPO)
    S = promote_type(eltype(O1), eltype(O2))
    T = typeof(O1)
    O = T.name.wrapper(S, length(O1))
    for (i, (A, B)) ∈ enumerate(zip(O1, O2))
        @matmul V[(x, a), σ, (y, b), η] := sum(γ) A[x, σ, y, γ] * B[a, γ, b, η]
        O[i] = V
    end
    O
end

Base.:(*)(A::AbstractTensorNetwork, B::AbstractTensorNetwork) = dot(A, B)


text:
# export    UpdatableQR,
#           add_column!,
#           add_column_householder!,
#           remove_column!,
#           update_views,
#           qraddcol!
#
# # Credit to https://github.com/mpf/QRupdate.jl/blob/master/src/QRupdate.jl
#
# # function qraddcol(A::AbstractMatrix{T}, Rin::AbstractMatrix{T}, a::Vector{T}, β::T = 0.0) where {T <:Real}
#
# function qraddcol!(A::Array{Float64,2}, R::UpperTriangular{Float64,Array{Float64,2}}, a::Array{Float64,1}, β::Float64 = 0.0)
#
#   m, n = size(A)
#   anorm  = norm(a)
#   anorm2 = anorm^2
#   β2  = β^2
#   if β != 0
#       anorm2 = anorm2 + β2
#       anorm  = sqrt(anorm2)
#   end
#
#   if n == 0
#       return reshape([anorm], 1, 1)
#   end
#
#   # R = UpperTriangular(Rin)
#
#   c      = A'*a           # = [A' β*I 0]*[a; 0; β]
#   u      = R'\c
#   unorm2 = norm(u)^2
#   d2     = anorm2 - unorm2
#
#   if d2 > anorm2 #DISABLE 0.01*anorm2     # Cheap case: γ is not too small
#       γ = sqrt(d2)
#   else
#       z = R\u          # First approximate solution to min ||Az - a||
#       r = a - A*z
#       c = A'*r
#       if β != 0
#           c = c - β2*z
#       end
#       du = R'\c
#       dz = R\du
#       z  += dz          # Refine z
#     # u  = R*z          # Original:     Bjork's version.
#       u  += du          # Modification: Refine u
#       r  = a - A*z
#       γ = norm(r)       # Safe computation (we know gamma >= 0).
#       if β != 0
#           γ = sqrt(γ^2 + β2*norm(z)^2 + β2)
#       end
#   end
#
#   # This seems to be faster than concatenation, ie:
#   # [ Rin         u
#   #   zeros(1,n)  γ ]
#   # Rout = Array{T}(n+1, n+1)
#   # Rout[1:n,1:n] = R
#   # Rout[1:n,n+1] = u
#   # Rout[n+1,n+1] = γ
#   # Rout[n+1,1:n] = 0.0
#   R = UpperTriangular(hcat(vcat(R, zeros(1,n)),zeros(n+1,1)))
#   @show size(R)
#   # Rout = UpperTriangular(zeros(n+1, n+1))
#
#   # view(Rout,1:n,1:n) .= R.data
#   view(R,1:n,n+1) .= u
#   R[n+1,n+1] = γ
#   # fill!(view(Rout, n+1,1:n), 0.0)
#   # Rout[n+1,1:n] = 0.0
#
#   return R
#   # nothing
# end
#
#
# # Credit to https://github.com/oxfordcontrol/GeneralQP.jl/blob/master/src/linear_algebra.jl
# mutable struct UpdatableQR{T} <: Factorization{T}
#     """
#     Gives the qr factorization an (n, m) matrix as Q1*R1
#     Q2 is such that Q := [Q1 Q2] is orthogonal and R is an (n, n) matrix where R1 "views into".
#     """
#     Q::Matrix{T}
#     R::Matrix{T}
#     n::Int
#     m::Int
#
#     Q1::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
#     Q2::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
#     R1::UpperTriangular{T, SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int},UnitRange{Int}}, false}}
#
#     function UpdatableQR(A::AbstractMatrix{T}) where {T}
#         n, m = size(A)
#         @assert(m <= n, "Too many columns in the matrix.")
#
#         F = qr(A)
#         Q = F.Q*Matrix(I, n, n)
#         R = zeros(T, n, n)
#         R[1:m, 1:m] .= F.R
#
#         new{T}(Q, R, n, m,
#             view(Q, :, 1:m), view(Q, :, m+1:n),
#             UpperTriangular(view(R, 1:m, 1:m)))
#     end
# end
#
# function add_column!(F::UpdatableQR{T}, a::AbstractVector{T}) where {T}
#     a1 = F.Q1'*a;
#     a2 = F.Q2'*a;
#
#     x = copy(a2)
#     for i = length(x):-1:2
#         G, r = givens(x[i-1], x[i], i-1, i)
#         lmul!(G, x)
#         lmul!(G, F.Q2')
#     end
#
#     F.R[1:F.m, F.m+1] .= a1
#     F.R[F.m+1, F.m+1] = x[1]
#
#     F.m += 1; update_views!(F)
#
#     return a2
# end
#
# function add_column_householder!(F::UpdatableQR{T}, a::AbstractVector{T}) where {T}
#     a1 = F.Q1'*a;
#     a2 = F.Q2'*a;
#
#     Z = qr(a2)
#     LAPACK.gemqrt!('R','N', Z.factors, Z.T, F.Q2) # Q2 .= Q2*F.Q
#     F.R[1:F.m, F.m+1] .= a1
#     F.R[F.m+1, F.m+1] = Z.factors[1, 1]
#     F.m += 1; update_views!(F)
#
#     return Z
# end
#
# function remove_column!(F::UpdatableQR{T}, idx::Int) where {T}
#     Q12 = view(F.Q, :, idx:F.m)
#     R12 = view(F.R, idx:F.m, idx+1:F.m)
#
#     for i in 1:size(R12, 1)-1
#         G, r = givens(R12[i, i], R12[i + 1, i], i, i+1)
#         lmul!(G, R12)
#         rmul!(Q12, G')
#     end
#
#     for i in 1:F.m, j in idx:F.m-1
#         F.R[i, j] = F.R[i, j+1]
#     end
#     F.R[:, F.m] .= zero(T)
#
#     F.m -= 1; update_views!(F)
#
#     return nothing
# end
#
# function update_views!(F::UpdatableQR{T}) where {T}
#     F.R1 = UpperTriangular(view(F.R, 1:F.m, 1:F.m))
#     F.Q1 = view(F.Q, :, 1:F.m)
#     F.Q2 = view(F.Q, :, F.m+1:F.n)
# end


text:
struct op_M1
    bi::Int64
    i::Int64
    f::Int64
    phase::Bool
    fac::Float64
end
struct op_E2
    bi::Int64
    i::Int64
    f::Int64
    phase::Bool
    fac::Float64
end

struct op_M1_nd
    bl::Int64
    Nl::Int64
    br::Int64
    Nr::Int64
    phase::Bool
    fac::Float64
end
struct op_E2_nd
    bl::Int64
    Nl::Int64
    br::Int64
    Nr::Int64
    phase::Bool
    fac::Float64
end

function tri_relation(ja,jb,jc)
    TF= true
    if ja+jb < jc ; TF=false;end
    if jb+jc < ja ; TF=false;end
    if jc+ja < jb ; TF=false;end
    if abs(ja-jb) > jc ; TF=false;end
    if abs(jb-jc) > ja ; TF=false;end
    if abs(jc-ja) > jb ; TF=false;end
    return TF
end

function init_ho_by_mass(A,ihwprm=1)
    nmass = [938.27231, 939.56563]
    amass = (nmass[1]+nmass[2])/2
    hc = 197.32696
    if ihwprm == 1 ## 41A^(-1/3) MeV
        hw = 41.0 * (A^(-1.0/3.0))
    elseif ihwprm == 2 ## 45A^(-1/3)-25A^(-2/3) MeV
        #J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).
        hw = 45.0 * (A^(-1.0/3.0)) -25.0 * (A^(-2.0/3.0))
    end
    bpar = hc/sqrt(hw*amass)
    return hw,bpar
end

function radius_power(k,n1,l1,n2,l2,bpar)
    s=0.0
    ll  = l1+l2+k
    ll1 = l2-l1+k
    ll2 = l1-l2+k
    ll1 = div(ll1,2)
    ll2 = div(ll2,2)

    imin = maximum([0, n1-ll1, n2-ll2])
    imax = minimum([n1, n2])
    for i = imin:imax
        s +=  binomial(ll1, n1-i) * binomial(ll2, n2-i) * (
            doublefactorial(ll+2*i+1)/doublefactorial(2*i))
    end
    s *= sqrt(doublefactorial(2*n1)/doublefactorial(2*n1+2*l1+1))
    s *= sqrt(doublefactorial(2*n2)/doublefactorial(2*n2+2*l2+1))
    s = Float64(s) * (1.0/sqrt(2^k)) * (-1)^(n1-n2) * bpar^k
    return s
end

function prep_obtr_diag(p_sps,n_sps,
                        mstates_p,mstates_n,tdims,
                        jocc_p,jocc_n,
                        pbits,nbits,bpar,
                        gfactors,eff_charge)
    ## Note that mu is in units of muN (= 0.105 efm)
    lp=length(p_sps); ln=length(n_sps); lblock = length(nbits)
    ec_p,ec_n = eff_charge
    c_mom_mu = sqrt(3.0/(4.0*pi))
    c_mom_Q = sqrt(5.0/(16.0*pi))

    ## M1 (mu; magnetic dipole moment)
    pjumps_mu = op_M1[  ]
    njumps_mu = op_M1[  ]
    ## E2 (Q; electric quadropole moment)
    pjumps_Q = op_E2[ ]
    njumps_Q = op_E2[ ]

    vec_p_ani = [ false for i=1:length(mstates_p)]
    vec_p_cre = [ false for i=1:length(mstates_p)]
    ex_r2 = 0.0 ; coeff= 0.0
    for i=1:lp
        ni,li,ji,tzi = p_sps[i]
        for f = i:lp
            nf,lf,jf,tzf = p_sps[f]
            ex_r2 = radius_power(2,ni,li,nf,lf,bpar)
            l6j = wigner6j(Float64,lf,li,1, ji//2,jf//2,1//2)
            s6j = wigner6j(Float64,1//2,1//2,1,ji//2,jf//2,lf)
            Q3j = wigner3j(Float64,ji//2,2,jf//2,1//2,0,-1//2)
            c,mc_s,mc_idxs = possible_mz(p_sps[i],mstates_p)
            a,ma_s,ma_idxs = possible_mz(p_sps[f],mstates_p)

            for mode in ["M1", "E2"]
                if mode == "M1" # selection rule
                    if (ni !=nf || li != lf); continue;end
                end
                @inbounds for (ic,mc) in enumerate(mc_s)
                    vec_p_ani .= false
                    vec_p_ani[mc_idxs[ic]] = true
                    bit_c = bitarr_to_int(vec_p_ani)
                    @inbounds for (ia,ma) in enumerate(ma_s)
                        if ma != mc; continue;end #M conservation
                        vec_p_cre .= false
                        vec_p_cre[ma_idxs[ia]] = true
                        bit_a = bitarr_to_int(vec_p_cre)
                        pfac = sqrt((ji+1)*(jf+1))
                        if mode =="M1"
                            lam = 1
                            sign1 = (-1)^((jf-ma)/2+lf+(ji+3)/2)
                            sign2 = (-1)^((jf-ma)/2+lf+(jf+3)/2)
                            fac1 = sign1 * l6j .* sqrt(lf*(lf+1)*(2*lf+1))
                            fac2 = sign2 * s6j .* sqrt(1.5)
                            m3j = pfac * wigner3j(jf//2,lam,ji//2,-ma//2,0,mc//2)
                            coeff = m3j * (fac1*gfactors[1]+fac2*gfactors[3]) *c_mom_mu
                        elseif mode =="E2"
                            lam = 2
                            signQ = (-1)^(jf + (1-ma)//2) * (1+(-1)^(li+lf))
                            if signQ == 0; continue;end
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,0,mc//2)
                            coeff = signQ * m3j *Q3j * pfac * ec_p * c_mom_Q * ex_r2
                        end
                        @inbounds for bi = 1:lblock
                            bf = bi # M-conservation
                            TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1]
                            idim = tdims[bi]
                            l_Nn = length(nbits[bi])
                            @inbounds for (Npi,pPhi) in enumerate(pbits[bi])
                                TF_connectable_1(pPhi,bit_a,bit_c,TF)
                                if TF[1]==false; continue;end
                                calc_phase_1!(pPhi,bit_a,bit_c,ret)
                                bisearch!(pbits[bf],ret[2],ridx)
                                if ridx[1] == 0;continue;end
                                Npf = ridx[1]
                                if mode =="M1"
                                    push!(pjumps_mu,
                                          op_M1(bi,Npi,Npf,ret[1]==-1,
                                                coeff .*ifelse(Npi==Npf,0.5,1.0)))
                                elseif mode =="E2"
                                    push!(pjumps_Q,
                                          op_E2(bi,Npi,Npf,ret[1]==-1,
                                                coeff .*ifelse(Npi==Npf,0.5,1.0)))
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    vec_n_ani = [ false for i=1:length(mstates_n)]
    vec_n_cre = [ false for i=1:length(mstates_n)]
    for i=1:ln
        ni,li,ji,tzi = n_sps[i]
        for f = i:ln
            nf,lf,jf,tzf = n_sps[f]
            ex_r2 = radius_power(2,ni,li,nf,lf,bpar)
            l6j = wigner6j(Float64,lf,li,1,ji//2,jf//2,1//2)
            s6j = wigner6j(Float64,1//2,1//2,1,ji//2,jf//2,lf)
            Q3j = wigner3j(Float64,ji//2,2,jf//2,1//2,0,-1//2)
            jc,mc_s,mc_idxs = possible_mz(n_sps[i],mstates_n)
            ja,ma_s,ma_idxs = possible_mz(n_sps[f],mstates_n)
            for mode in ["M1", "E2"]
                if mode == "M1"
                    if (ni !=nf || li != lf)
                        continue
                    end
                end
                @inbounds for (ic,mc) in enumerate(mc_s)
                    vec_n_ani .= false
                    vec_n_ani[mc_idxs[ic]] = true
                    bit_c = bitarr_to_int(vec_n_ani)
                    if bit_c == 0; continue; end
                    @inbounds for (ia,ma) in enumerate(ma_s)
                        if ma != mc; continue;end #M conservation
                        #if a == c && ma > mc; continue;end
                        vec_n_cre .= false
                        vec_n_cre[ma_idxs[ia]] = true
                        bit_a = bitarr_to_int(vec_n_cre)
                        if bit_a == 0; continue; end

                        pfac = sqrt((ji+1)*(jf+1))
                        sign1 = (-1)^((jf-ma)/2+lf+(ji+3)/2)
                        sign2 = (-1)^((jf-ma)/2+lf+(jf+3)/2)
                        if mode =="M1"
                            fac1 = sign1 * pfac .* sqrt(lf*(lf+1)*(2*lf+1)) * l6j
                            fac2 = sign2 * pfac * sqrt(1.5) * s6j
                            lam = 1
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,0,mc//2)
                            coeff = m3j * (fac1*gfactors[2]+fac2*gfactors[4]) * c_mom_mu
                        elseif mode=="E2"
                            lam = 2
                            signQ = (-1)^(jf + (1-ma)//2) * (1+(-1)^(li+lf))
                            if signQ == 0; continue;end
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,0,mc//2)
                            coeff = signQ * m3j *Q3j * pfac * ec_n * c_mom_Q * ex_r2
                            #println("fE2 $coeff sign $signQ 3j $m3j jf $jf ji $ji ma $ma mc $mc")
                            if coeff == 0.0; continue;end
                        end
                        @inbounds for bi = 1:lblock
                            bf = bi # M-conservation
                            TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1]
                            nbit = nbits[bi]
                            idim = tdims[bi]; l_Nn = length(nbit)
                            @inbounds for (Nni,nPhi) in enumerate(nbit)
                                TF_connectable_1(nPhi,bit_a,bit_c,TF)
                                if TF[1]==false; continue;end
                                calc_phase_1!(nPhi,bit_a,bit_c,ret)
                                bisearch!(nbits[bf],ret[2],ridx)
                                if ridx[1] == 0;continue;end
                                Nnf = ridx[1]
                                if mode =="M1"
                                    push!(njumps_mu,
                                          op_M1(bi,Nni,Nnf,ret[1]==-1,
                                                coeff .*ifelse(Nni==Nnf,0.5,1.0)))
                                elseif mode =="E2"
                                    push!(njumps_Q,
                                          op_E2(bi,Nni,Nnf,ret[1]==-1,
                                                coeff .*ifelse(Nni==Nnf,0.5,1.0)))
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return pjumps_mu,njumps_mu, pjumps_Q,njumps_Q
end

function op_obtr_diag(p_oprts,n_oprts,pbits,nbits,tdims,wf,twf)
    @inbounds for tmp in p_oprts
        bi = tmp.bi; Npi = tmp.i; Npf = tmp.f
        coeff = ifelse(tmp.phase,-tmp.fac,tmp.fac)
        pbit = pbits[bi]; nbit = nbits[bi]
        idim = tdims[bi]
        l_Nn = length(nbit)
        l_Np = length(pbit)
        tMi = idim - l_Nn + Npi *l_Nn
        tMf = idim - l_Nn + Npf *l_Nn
        @inbounds for nidx = 1:l_Nn
            Mi = tMi + nidx
            Mf = tMf + nidx
            twf[Mf] += coeff .* wf[Mi]
            twf[Mi] += coeff .* wf[Mf]
        end
    end
    for tmp in n_oprts
        bi = tmp.bi
        Nni = tmp.i; Nnf = tmp.f
        phase = tmp.phase
        coeff = ifelse(phase,-tmp.fac,tmp.fac)
        pbit = pbits[bi]; nbit = nbits[bi]
        idim = tdims[bi]
        l_Nn = length(nbit)
        l_Np = length(pbit)
        tMi = idim - l_Nn + Nni
        tMf = idim - l_Nn + Nnf
        @inbounds for pidx = 1:l_Np
            Mi = tMi + pidx * l_Nn
            Mf = tMf + pidx * l_Nn
            twf[Mf] += coeff .* wf[Mi]
            twf[Mi] += coeff .* wf[Mf]
        end
    end
    return nothing
end

function prep_obtr_nondiag(jl2,jr2,p_sps,n_sps,
                           mstates_p,mstates_n,
                           jocc_p,jocc_n,
                           tdims_l,pbits_l,nbits_l,
                           tdims_r,pbits_r,nbits_r,
                           p_bifs,n_bifs,omu,bpar,
                           gfactors,eff_charge)
    lp=length(p_sps); ln=length(n_sps);
    lblock_l = length(nbits_l)
    lblock_r = length(nbits_r)
    ec_p,ec_n = eff_charge
    c_mom_mu = sqrt(3.0/(4.0*pi))
    c_mom_Q = sqrt(5.0/(16.0*pi))

    ## M1 (mu; magnetic dipole moment)
    pjumps_mu = op_M1_nd[ ]
    njumps_mu = op_M1_nd[ ]
    ## E2 (Q; electric quadropole moment)
    pjumps_Q = op_E2_nd[ ]
    njumps_Q = op_E2_nd[ ]

    vec_p_ani = [ false for i=1:length(mstates_p)]
    vec_p_cre = [ false for i=1:length(mstates_p)]
    ex_r2 = 0.0 ; coeff= 0.0
    for i=1:lp
        ni,li,ji,tzi = p_sps[i]
        for f = 1:lp
            nf,lf,jf,tzf = p_sps[f]
            ex_r2 = radius_power(2,ni,li,nf,lf,bpar)
            l6j = wigner6j(Float64,lf,li,1,ji//2,jf//2,1//2)
            s6j = wigner6j(Float64,1//2,1//2,1,ji//2,jf//2,lf)
            Q3j = wigner3j(Float64,ji//2,2,jf//2,1//2,0,-1//2)
            c,mc_s,mc_idxs = possible_mz(p_sps[i],mstates_p)
            a,ma_s,ma_idxs = possible_mz(p_sps[f],mstates_p)
            for mode in ["M1", "E2"]
                if mode == "M1" # slection rule
                    if (ni !=nf || li != lf); continue;end
                    if div(jl2-jr2,2) > 1;continue;end
                elseif mode == "E2"
                    if div(jl2-jr2,2) > 2; continue;end
                end
                @inbounds for (ic,mc) in enumerate(mc_s)
                    vec_p_ani .= false
                    vec_p_ani[mc_idxs[ic]] = true
                    bit_c = bitarr_to_int(vec_p_ani)
                    @inbounds for (ia,ma) in enumerate(ma_s)
                        if div(ma - mc,2) != omu ; continue;end #M conservation
                        vec_p_cre .= false
                        vec_p_cre[ma_idxs[ia]] = true
                        bit_a = bitarr_to_int(vec_p_cre)
                        pfac = sqrt((ji+1)*(jf+1))
                        if mode =="M1"
                            lam = 1
                            sign1 = (-1)^((jf-ma)/2+lf+(ji+3)/2)
                            sign2 = (-1)^((jf-ma)/2+lf+(jf+3)/2)
                            fac1 = sign1 * l6j .* sqrt(lf*(lf+1)*(2*lf+1))
                            fac2 = sign2 * s6j .* sqrt(1.5)
                            m3j = pfac * wigner3j(jf//2,lam,ji//2,-ma//2,omu,mc//2)
                            coeff = m3j * (fac1*gfactors[1]+fac2*gfactors[3]) *c_mom_mu
                        elseif mode =="E2"
                            lam = 2
                            signQ = (-1)^(jf + (1-ma)//2) * (1+(-1)^(li+lf))
                            if signQ == 0; continue;end
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,omu,mc//2)
                            coeff = signQ * m3j *Q3j * pfac * ec_p * c_mom_Q * ex_r2
                        end
                        @inbounds for (bl,br) in p_bifs
                            TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1]
                            @inbounds for (Npi,pPhi) in enumerate(pbits_r[br])
                                TF_connectable_1(pPhi,bit_a,bit_c,TF)
                                if TF[1]==false; continue;end
                                calc_phase_1!(pPhi,bit_a,bit_c,ret)
                                bisearch!(pbits_l[bl],ret[2],ridx)
                                if ridx[1] == 0;continue;end
                                Npf = ridx[1]
                                if mode =="M1"
                                    push!(pjumps_mu,
                                          op_M1_nd(bl,Npf,br,Npi,ret[1]==-1,coeff))
                                elseif mode =="E2"
                                    push!(pjumps_Q,
                                          op_E2_nd(bl,Npf,br,Npi,ret[1]==-1,coeff))
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    vec_n_ani = [ false for i=1:length(mstates_n)]
    vec_n_cre = [ false for i=1:length(mstates_n)]
    for i=1:ln
        ni,li,ji,tzi = n_sps[i]
        for f = 1:ln
            nf,lf,jf,tzf = n_sps[f]
            ex_r2 = radius_power(2,ni,li,nf,lf,bpar)
            l6j = wigner6j(Float64,lf,li,1,ji//2,jf//2,1//2)
            s6j = wigner6j(Float64,1//2,1//2,1,ji//2,jf//2,lf)
            Q3j = wigner3j(Float64,ji//2,2,jf//2,1//2,0,-1//2)
            jc,mc_s,mc_idxs = possible_mz(n_sps[i],mstates_n)
            ja,ma_s,ma_idxs = possible_mz(n_sps[f],mstates_n)
            for mode in ["M1", "E2"]
                if mode == "M1"
                    if (ni !=nf || li != lf);continue;end
                    if div(jl2-jr2,2) > 1;continue;end
                elseif mode == "E2"
                    if div(jl2-jr2,2) > 2; continue;end
                end
                @inbounds for (ic,mc) in enumerate(mc_s)
                    vec_n_ani .= false
                    vec_n_ani[mc_idxs[ic]] = true
                    bit_c = bitarr_to_int(vec_n_ani)
                    if bit_c == 0; continue; end
                    @inbounds for (ia,ma) in enumerate(ma_s)
                        if div(ma -mc,2) != omu; continue;end #M conservation
                        vec_n_cre .= false
                        vec_n_cre[ma_idxs[ia]] = true
                        bit_a = bitarr_to_int(vec_n_cre)
                        if bit_a == 0; continue;end
                        pfac = sqrt((ji+1)*(jf+1))
                        sign1 = (-1)^((jf-ma)/2+lf+(ji+3)/2)
                        sign2 = (-1)^((jf-ma)/2+lf+(jf+3)/2)
                        if mode =="M1"
                            fac1 = sign1 * pfac .* sqrt(lf*(lf+1)*(2*lf+1)) * l6j
                            fac2 = sign2 * pfac * sqrt(1.5) * s6j
                            lam = 1
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,omu,mc//2)
                            coeff = m3j * (fac1*gfactors[2]+fac2*gfactors[4]) * c_mom_mu
                        elseif mode=="E2"
                            lam = 2
                            signQ = (-1)^(jf + (1-ma)//2) * (1+(-1)^(li+lf))
                            if signQ == 0; continue;end
                            m3j =  wigner3j(Float64,jf//2,lam,ji//2,-ma//2,omu,mc//2)
                            coeff = signQ * m3j *Q3j * pfac * ec_n * c_mom_Q * ex_r2
                        end
                        @inbounds for (bl,br) in n_bifs
                            TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1]
                            @inbounds for (Nni,nPhi) in enumerate(nbits_r[br])
                                TF_connectable_1(nPhi,bit_a,bit_c,TF)
                                if TF[1]==false; continue;end
                                calc_phase_1!(nPhi,bit_a,bit_c,ret)
                                bisearch!(nbits_l[bl],ret[2],ridx)
                                if ridx[1] == 0;continue;end
                                Nnf = ridx[1]
                                if mode =="M1"
                                    push!(njumps_mu,
                                          op_M1_nd(bl,Nnf,br,Nni,ret[1]==-1,coeff))
                                elseif mode =="E2"
                                    push!(njumps_Q,
                                          op_E2_nd(bl,Nnf,br,Nni,ret[1]==-1,coeff))
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return pjumps_mu,njumps_mu,pjumps_Q,njumps_Q
end

function op_obtr_nondiag(p_oprts,n_oprts,
                         pbits_l,nbits_l,tdims_l,
                         pbits_r,nbits_r,tdims_r,
                         wf,twf)
    @inbounds for tmp in p_oprts
        bi = tmp.br; Npi = tmp.Nr
        bf = tmp.bl; Npf = tmp.Nl
        coeff = ifelse(tmp.phase,-tmp.fac,tmp.fac)
        nbit_r = nbits_r[bi]; nbit_l = nbits_l[bf]
        idim = tdims_r[bi]
        fdim = tdims_l[bf]
        l_Nn_r = length(nbit_r)
        l_Nn_l = length(nbit_l)
        tMi = idim - l_Nn_r + Npi *l_Nn_r
        tMf = fdim - l_Nn_l + Npf *l_Nn_l
        @inbounds for nidx_r = 1:l_Nn_r
            Mi = tMi + nidx_r
            n_r = nbit_r[nidx_r]
            @inbounds for nidx_l = 1:l_Nn_l
                n_l = nbit_l[nidx_l]
                if n_r != n_l; continue;end
                Mf = tMf + nidx_l
                twf[Mf] += coeff .* wf[Mi]
            end
        end
    end
    @inbounds for tmp in n_oprts
        bi = tmp.br; Nni = tmp.Nr
        bf = tmp.bl; Nnf = tmp.Nl
        coeff = ifelse(tmp.phase,-tmp.fac,tmp.fac)
        pbit_r = pbits_r[bi]; pbit_l = pbits_l[bf]
        nbit_r = nbits_r[bi]; nbit_l = nbits_l[bf]
        idim = tdims_r[bi]
        fdim = tdims_l[bf]
        l_Nn_r = length(nbit_r)
        l_Nn_l = length(nbit_l)
        l_Np_r = length(pbit_r)
        l_Np_l = length(pbit_l)
        tMi = idim - l_Nn_r + Nni
        tMf = fdim - l_Nn_l + Nnf
        @inbounds for pidx_r = 1:l_Np_r
            Mi = tMi + pidx_r * l_Nn_r
            p_r = pbit_r[pidx_r]
            @inbounds for pidx_l = 1:l_Np_l
                p_l = pbit_l[pidx_l]
                if p_r != p_l; continue;end
                Mf = tMf + pidx_l * l_Nn_l
                twf[Mf] += coeff .* wf[Mi]
            end
        end
    end
    return nothing
end

function calc_bifs(Mps_l,Mns_l,Mtot_l,Mps_r,Mns_r,Mtot_r)
    # for right wav
    blocks_r = [ [0,0] ]; deleteat!(blocks_r,1)
    for Mp in Mps_r
        for Mn in Mns_r
            if Mp + Mn == Mtot_r
                push!(blocks_r,[Mp,Mn])
            end
        end
    end
    blocks_l = [ [0,0] ]; deleteat!(blocks_l,1)
    for Mp in Mps_l
        for Mn in Mns_l
            if Mp + Mn == Mtot_l
                push!(blocks_l,[Mp,Mn])
            end
        end
    end
    p_bifs = [ [0,0] ]; deleteat!(p_bifs,1)
    n_bifs = [ [0,0] ]; deleteat!(n_bifs,1)

    for (ir,bl_r) in enumerate(blocks_r)
        Mpr,Mnr = bl_r
        for (il,bl_l) in enumerate(blocks_l)
            Mpl,Mnl = bl_l
            if Mnr == Mnl
                push!(p_bifs,[il,ir])
            end
            if Mpr == Mpl
                push!(n_bifs,[il,ir])
            end
        end
    end
    return p_bifs,n_bifs
end

"""
    transit_main(sntf,target_nuc,jl2,jr2,in_wfs;
                 num_ev_l=100,num_ev_r=100,q=1,is_block=false,is_show=true,
                 calc_EM=true,gfactors=[1.0,0.0,5.586,-3.826],eff_charge=[1.5,0.5])

Calculate the M1&E2 transitions for two wavefunctions

# Arguments
- `sntf`:path to the interaction file
- `target_nuc`:target nucleus in string e.g., "Si28"
- `jl2`:J*2 for the left w.f.
- `jr2`:J*2 for the right w.f.
- `in_wfs`:["path to left wf","path to right wf"]

# Optional arguments
- `num_ev_l=100`:upper limit of the number of eigenvectors for the left w.f.
- `num_ev_r=100`:upper limit of the number of eigenvectors for the right w.f.
- `is_show=true`:to display the TimerOutput
- `gfactors=[1.0,0.0,5.586,-3.826]`:g factors [glp,gln,gsp,gsn]
- `eff_charge=[1.5,0.5]`:effective charges [ep,en]

# Optional arguments (not used, but for future developments)
- `q=1`:block size
- `is_block=false`:to use Block algorithm
"""
function transit_main(sntf,target_nuc,jl2,jr2,in_wfs;
                      num_ev_l=100,num_ev_r=100,
                      q=1,is_block=false,is_show=true,
                      calc_EM=true,
                      gfactors=[1.0,0.0,5.586,-3.826],
                      eff_charge=[1.5,0.5])
    if length(in_wfs) == 0;println("input wfs must be specified");exit();end

    to = TimerOutput()
    Anum = parse(Int64, match(reg,target_nuc).match)
    lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsnt(sntf,Anum)
    hw, bpar = init_ho_by_mass(Anum,1) # mass formula
    if 16 <= Anum <= 40
        hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell
    end

    target_el = replace.(target_nuc, string(Anum)=>"")
    Z,N,vp,vn = getZNA(target_el,Anum,cp,cn)
    mstates_p, mstates_n,mz_p,mz_n = def_mstates(p_sps,n_sps)

    widx = [1,2]
    if jl2 < jr2
        tr = jl2; tl = jr2; jl2 = tl; jr2 = tr; widx=[2,1]
    end

    Mtot_l = jl2; Mtot_r = jr2
    pbits_l,nbits_l,jocc_p,jocc_n,Mps_l,Mns_l,tdims_l = occ(p_sps,mstates_p,mz_p,vp,
                                                            n_sps,mstates_n,mz_n,vn,
                                                            Mtot_l)
    pbits_r,nbits_r,jocc_p,jocc_n,Mps_r,Mns_r,tdims_r = occ(p_sps,mstates_p,mz_p,vp,
                                                            n_sps,mstates_n,mz_n,vn,
                                                            Mtot_r)
    mdim_l = tdims_l[end];   mdim_r = tdims_r[end]
    println("wfl: ",in_wfs[widx[1]])
    if in_wfs[widx[1]] != in_wfs[widx[2]]; println("wfr: ",in_wfs[widx[2]]);end
    println("mdims J2=$jl2 $mdim_l J2=$jr2 $mdim_r Mtots=($Mtot_l, $Mtot_r)")
    wfs_l,jjs_l = read_wav(in_wfs[widx[1]],mdim_l,1;all=true)
    wfs_r,jjs_r = read_wav(in_wfs[widx[2]],mdim_r,1;all=true)
    lblock_l=length(pbits_l); lblock_r=length(pbits_r)
    omu = div(jl2-jr2,2)
    if omu > 2; prinln("Only the M1&E2 are supported now");return nothing;end

    p_bifs,n_bifs = calc_bifs(Mps_l,Mns_l,Mtot_l,Mps_r,Mns_r,Mtot_r)
    # operation of EM transition operators O|Phi>
    @timeit to "prep" begin
        op_p_mu,op_n_mu,op_p_Q,op_n_Q = prep_obtr_nondiag(jl2,jr2,p_sps,n_sps,
                                                          mstates_p,mstates_n,
                                                          jocc_p,jocc_n,
                                                          tdims_l,pbits_l,nbits_l,
                                                          tdims_r,pbits_r,nbits_r,
                                                          p_bifs,n_bifs,omu,bpar,
                                                          gfactors,eff_charge)
    end

    ## eval expectation values
    c_mom_mu = sqrt(4.0*pi/3.0)
    c_mom_Q  = sqrt(16.0*pi/5.0)

    if tri_relation(jl2//2,jr2//2,1)
        @timeit to "M1" begin
        #println("M1 available")
        twf = zeros(Float64,mdim_l)
        lam = 1
        Jf = jl2//2; Ji=jr2//2
        Mf=Mtot_l//2;Mi=Mtot_r//2
        sel = wigner3j(Jf,lam,Ji,-Mf,omu,Mi)
        signmu = (-1)^(Jf-Mf)
        cg = wigner3j(Ji,lam,Jf,-Ji,omu,Jf)
        for (nth_r,wf_r) in enumerate(wfs_r)
            if sel == 0; continue;end
            if nth_r > num_ev_r; continue;end
            twf .= 0.0
            op_obtr_nondiag(op_p_mu,op_n_mu,
                            pbits_l,nbits_l,tdims_l,
                            pbits_r,nbits_r,tdims_r,
                            wf_r,twf)
            if calc_EM
                for (nth_l,wf_l) in enumerate(wfs_l)
                    if nth_l > num_ev_l; continue;end
                    Mred =  dot(wf_l,twf) /sel
                    tx = ""
                    tx *= "Mred.(M1) ($nth_l,$nth_r) " * @sprintf("%9.3f ",Mred) *"  "
                    B_l = Mred^2 / (2*Ji+1); B_r = Mred^2 / (2*Jf+1)
                    tx *= "B(M1) ($nth_l<=$nth_r) " * @sprintf("%9.3f ",B_l) *" B(M1) ($nth_l=>$nth_r) " * @sprintf("%9.3f ",B_r)
                    if in_wfs[1] == in_wfs[2] && nth_l==nth_r
                        gm = Mred * cg *c_mom_mu * signmu
                        tx *= " <gm>= "* @sprintf("%9.3f ",gm)
                    end
                    println(tx)
                end
            else
                nth_l = nth_r
                wf_l= wfs_l[nth_l]
                Mred =  dot(wf_l,twf) /sel
                gm = Mred * cg *c_mom_mu * signmu
                tx = " <gm>= "* @sprintf("%9.3f ",gm)
                println(tx)
            end
        end
        end
    end
    if tri_relation(jl2//2,jr2//2,2)
        #println("E2 available")
        @timeit to "E2" begin
        twf = zeros(Float64,mdim_l)
        lam = 2
        Jf = jl2//2; Mf = Mtot_l//2
        Ji = jr2//2; Mi = Mtot_r//2
        sel = wigner3j(Jf,lam,Ji,-Mf,omu,Mi)
        signQ = (-1)^(Jf-Mf)
        cg = wigner3j(Ji,lam,Jf,-Ji,omu,Jf)
        for (nth_r,wf_r) in enumerate(wfs_r)
            if sel==0;continue;end
            if nth_r > num_ev_r; continue;end
            twf .= 0.0
            op_obtr_nondiag(op_p_Q,op_n_Q,
                            pbits_l,nbits_l,tdims_l,
                            pbits_r,nbits_r,tdims_r,
                            wf_r,twf)
            if calc_EM
                for (nth_l,wf_l) in enumerate(wfs_l)
                    if nth_l > num_ev_l; continue;end
                    Mred =  dot(wf_l,twf) /sel
                    tx = ""
                    tx *= "Mred.(E2) ($nth_l,$nth_r) " * @sprintf("%9.3f ",Mred) *"  "
                    B_l = Mred^2 / (2*Ji+1); B_r = Mred^2 / (2*Jf+1)
                    tx *= "B(E2) ($nth_l<=$nth_r) " * @sprintf("%9.3f ",B_l) *" B(E2) ($nth_l=>$nth_r) " * @sprintf("%9.3f ",B_r)
                    if in_wfs[1] == in_wfs[2] && nth_l==nth_r
                        Q = Mred * cg *c_mom_Q * signQ
                        tx *= " <eQ>= "* @sprintf("%9.3f ",Q)
                    end
                    println(tx)
                end
            else
                wf_l = wfs_l[nth_r]
                Mred =  dot(wf_l,twf) /sel
                Q = Mred * cg *c_mom_Q * signQ
                tx = " <eQ>= "* @sprintf("%9.3f ",Q)
                println(tx)
            end
        end
        end
    end
    if is_show
        show(to, allocations = true,compact = false)
    end
    println("")
    return nothing
end

function eval_moment(Mtot,Rvecs,totJs,
                     p_sps,n_sps,
                     mstates_p,mstates_n,tdims,
                     jocc_p,jocc_n,pbits,nbits,bpar,
                     gfactors,eff_charge)
    mdim = tdims[end]
    c_mom_mu = sqrt(4.0*pi/3.0)
    c_mom_Q = sqrt(16.0*pi/5.0)
    op_p_mu,op_n_mu,op_p_Q,op_n_Q = prep_obtr_diag(p_sps,n_sps,
                                                   mstates_p,mstates_n,
                                                   tdims,jocc_p,jocc_n,
                                                   pbits,nbits,bpar,
                                                   gfactors,eff_charge)
    tx_mom = ""
    twf1 = zeros(Float64,mdim)
    twf2 = zeros(Float64,mdim)
    Mi = Mtot//2; Mf = Mtot//2
    for (nth,wf) in enumerate(Rvecs)
        twf1 .= 0.0 ; twf2 .= 0.0
        # for M1
        p_oprts = op_p_mu;  n_oprts = op_n_mu
        op_obtr_diag(p_oprts,n_oprts,pbits,nbits,tdims,wf,twf1)
        # for E2
        p_oprts = op_p_Q;  n_oprts = op_n_Q
        op_obtr_diag(p_oprts,n_oprts,pbits,nbits,tdims,wf,twf2)

        Ji = Int(2*totJs[nth]) // 2;
        for np = nth:length(Rvecs)
            tx = ""
            wf_l = Rvecs[np]
            Jf = Int(2*totJs[np])//2
            lam = 1
            omu = div(Mf-Mi,2)
            sel = wigner3j(Jf,lam,Ji,-Mf,omu,Mi)
            if sel != 0
                signmu = (-1)^(Jf-Mf)
                cg = wigner3j(Ji,lam,Jf,-Ji,omu,Jf) / sel
                if np ==nth
                    gm = signmu * dot(twf1,wf) *cg  * c_mom_mu
                    tx *= "mu (n=$nth) " * @sprintf("%12.5f ",gm) *"\t"
                else
                    Mred =  dot(wf_l,twf1) /sel
                    B_l = Mred^2 / (2*Ji+1); B_r = Mred^2 / (2*Jf+1)
                    #tx *= "Mred.(M1) ($np,$nth) " * @sprintf("%9.3f ",Mred) *"  "
                    #tx *= "B(M1) ($np<=$nth) " * @sprintf("%9.3f ",B_l)
                    #tx *= " B(M1) ($np=>$nth) " * @sprintf("%9.3f ",B_r)
                end
            end
            lam = 2
            sel =  wigner3j(Jf,lam,Ji,-Mf,omu,Mi)
            if sel != 0
                signQ = ( (-1)^(Jf-Mf))
                cg = wigner3j(Ji,lam,Jf,-Ji,omu,Jf) / sel
                if np == nth
                    Q = dot(twf2,wf) * cg * c_mom_Q * signQ
                    tx *= "Q (n=$nth) " * @sprintf("%12.5f ",Q)
                else
                    Mred = dot(wf_l,twf2) /sel
                    B_l = Mred^2 / (2*Ji+1); B_r = Mred^2 / (2*Jf+1)
                    #tx *= "Mred.(E2) ($np,$nth) " * @sprintf("%9.3f ",Mred) *"  "
                    #tx *= "B(E2) ($np<=$nth) " * @sprintf("%9.3f ",B_l)
                    #tx *= " B(E2) ($np=>$nth) " * @sprintf("%9.3f ",B_r)
                end
            end
            if tx !="";tx_mom *= tx*"\n";end
        end
    end
    return tx_mom
end


text:
using Plots
pyplot()

# struct to hold all parameters for simulation
struct Param
    xmax::Float64
    res::Int64
    dt::Float64
    timesteps::Int64
    dx::Float64
    x::Vector{Float64}
    dk::Float64
    k::Vector{Float64}

    Param() = new(10.0, 512, 0.05, 1000, 2 * 10.0/512,
                  Vector{Float64}(-10.0 + 10.0/512 : 20.0/512 : 10.0),
                   pi / 10.0,
                  Vector{Float64}(vcat(0:512/2 - 1, -512/2 : -1) * pi/10.0))
    Param(xmax::Float64, res::Int64, dt::Float64, timesteps::Int64) = new(
              xmax, res, dt, timesteps,
              2*xmax/res, Vector{Float64}(-xmax+xmax/res:2*xmax/res:xmax),
              pi/xmax, Vector{Float64}(vcat(0:res/2-1, -res/2:-1)*pi/xmax)
          )
end

# struct to hold all operators
mutable struct Operators
    V::Vector{Complex{Float64}}
    PE::Vector{Complex{Float64}}
    KE::Vector{Complex{Float64}}
    wfc::Vector{Complex{Float64}}
end

# Function to initialize the wfc and potential
function init(par::Param, voffset::Float64, wfcoffset::Float64)
    V = 0.5 * (par.x - voffset).^2
    wfc = 3 * exp.(-(par.x - wfcoffset).^2/2)
    PE = exp.(-0.5*im*V*par.dt)
    KE = exp.(-0.5*im*par.k.^2*par.dt)

    opr = Operators(V, PE, KE, wfc)
end

# Function for the split-operator loop
function split_op(par::Param, opr::Operators)

    for i = 1:par.timesteps
        # Half-step in real space
        opr.wfc = opr.wfc .* opr.PE

        # fft to phase space
        opr.wfc = fft(opr.wfc)

        # Full step in phase space
        opr.wfc = opr.wfc .* opr.KE

        # ifft back
        opr.wfc = ifft(opr.wfc)

        # final half-step in real space
        opr.wfc = opr.wfc .* opr.PE

        # plotting density and potential
        density = abs2.(opr.wfc)

        plot([density, real(opr.V)])
        savefig("density" * string(lpad(i, 5, 0)) * ".png")
        println(i)
    end
end

# main function
function main()
    par = Param(10.0, 512, 0.05, 1000)
    opr = init(par, 0.0, 1.0)
    split_op(par, opr)
end

main()

max_stars_repo_name: sanstorik/algorithm-archive

max_issues_repo_name: sanstorik/algorithm-archive

max_forks_repo_name: sanstorik/algorithm-archive


text:
__precompile__()
module RigidBodyRLEnvironments

using ReinforcementLearning
import ReinforcementLearning: step!, reset!, Environment

using POMDPs
using RigidBodyDynamics
using RigidBodyDynamics.OdeIntegrators
using StaticArrays
using RecipesBase

export

    CartpoleSimulator,
    CartpoleEnvironment,
    reset!,
    step!

abstract type RigidBodySimulator end

include("utils.jl")
include("cartpole/cartpole_sim.jl")
include("cartpole/cartpole_env.jl")
include("cartpole/cartpole_recipe.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ e5b53410-63cf-11eb-2532-358a68b59eb9
md"# Variables and basic types in Julia"

# ╔═╡ 16f1e0f0-63d0-11eb-2ee4-8f8be24866e8
md"## Variables"

# ╔═╡ 31972af0-63d0-11eb-3e22-bb3370a9e40d
md"""A variable is simply a name associated to any type of value.

To declare a new variable, we use `=`, also named assignment operator:"""

# ╔═╡ 9d83cca0-63d0-11eb-0346-751060e8cab5
a = 1

# ╔═╡ ffed727d-b072-4797-9667-dc3a1c7e0109
a

# ╔═╡ ba8fdd70-63d0-11eb-0565-f3e77f3f802f
md"Instead of repeating the value, we can then use the variable directly:"

# ╔═╡ a8247100-63d0-11eb-0442-11949579c4ba
b = a + 2

# ╔═╡ f1641eb0-63d0-11eb-2ed0-d5ab1ce695ae
md"""You can name your variables as you want, except for some few reserved words such as `true`, `false`, `if`, `else`, because they are part of the language itself. We will see what they mean later. You can find the list of reserved keywords [here](https://docs.julialang.org/en/v1/base/base/#Keywords).

Julia also allow using **Unicode characters**, which is awesome!

See by yourself:
"""

# ╔═╡ e20a69a0-63d1-11eb-3f1e-159164d5e71d
α = 2

# ╔═╡ 2a77d330-63d2-11eb-16c3-d3be6a6d9975
md"α is written by typing `\alpha` and then `ctrl+enter`.
`ctrl+enter` is used for autocompletion. We can even use emojis:
"

# ╔═╡ f85e81f0-63d1-11eb-2385-0196982bc487
🌴 = 25.0

# ╔═╡ 8ab14f10-63d2-11eb-3874-073f07427a26
md"And it is valid Julia code."

# ╔═╡ 257727a2-63d2-11eb-311b-0b14f3d0f747
α + 🌴

# ╔═╡ b6873b92-63d2-11eb-1a56-1d01d825cf01
md"As a sidenote, it is a convention to name your variable using lower case, and to separate words by underscores if needed, *e.g.*:"

# ╔═╡ e76d15e0-63d2-11eb-111e-8b0637007258
	savings = 5500.0

# ╔═╡ fa35e490-63d2-11eb-210d-0b44aeb5cbfa
bike_price = 5000.0

# ╔═╡ 68fb21b0-63d3-11eb-1792-45baa321e122
remaining_money = savings - bike_price

# ╔═╡ 6f7a56f0-63d3-11eb-113b-b11ac60caaed
md"If I buy this bike, I'll have $(remaining_money) € or $ or whatever on my account."

# ╔═╡ 8813c0d0-63d2-11eb-1723-3b802f8a9f54
md"## Basic types"

# ╔═╡ a067ae80-63d7-11eb-3d1e-39eae17d3434
md"### Listing basic types"

# ╔═╡ f4b09520-63d1-11eb-25d9-7564a91fa7d5
md"""
Values are all defined by a type. For example `1` is called an integer. We can see the type of a variable using `typeof()`

The most basic types are:

- Integers, numbers with no fractional part:
"""

# ╔═╡ ba6332b0-63d1-11eb-17a5-0160686e19e7
cars = 10

# ╔═╡ 6fa8ab30-63d9-11eb-0ce9-0f64d74d9b96
md"And the actual name of the type is:"

# ╔═╡ fe28579e-63d6-11eb-31c6-59b910298cc2
typeof(cars)

# ╔═╡ ccd15a10-63d8-11eb-2cce-9b5b18f2f003
md"Note it is called an `Int64`, not just `Int`. This is because my machine has a 64-bit architecture. If you run on a 32-bits machine, you would get `Int32` by default. This is all technical details, but if you're interested you can find more information [here](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/)"

# ╔═╡ 3279b400-63d6-11eb-093a-4539f1254a83
md"- Floating point numbers, meaning numbers with a fractional part:"

# ╔═╡ 3d9d6b0e-63d6-11eb-2cc1-39a76f290513
max_speed = 250.0

# ╔═╡ 0d5dc660-63d7-11eb-0458-29b554125fce
typeof(max_speed)

# ╔═╡ 53e9e240-63d6-11eb-00dc-098c64e573f1
md"- Characters, just one character defined using `'...'`:"

# ╔═╡ 61ffb030-63d6-11eb-059d-e912b2c828fb
rank = 'A'

# ╔═╡ 1204e040-63d7-11eb-2c61-e3f212329a44
typeof(rank)

# ╔═╡ 72e5b520-63d6-11eb-0463-3bd1d3643d15
md"""- Strings, one or several characters defined using `"..."`:"""

# ╔═╡ 7b9537e0-63d6-11eb-1d93-d3e253614910
brand = "Tesla"

# ╔═╡ 19752c40-63d7-11eb-0500-11c3318c9319
typeof(brand)

# ╔═╡ 48cab050-63dc-11eb-2f06-651fa9acbd54
md"""NB: A single character declared using `"..."` instead of `'...'` is still a String, not a character:"""

# ╔═╡ 69283610-63dc-11eb-1d97-25d73d0ce439
typeof("T")

# ╔═╡ 27efd680-63d7-11eb-0701-67d7c4fc2847
md"- Booleans, logical type, either `true` or `false`:"

# ╔═╡ 5c08ef5e-63d7-11eb-34d0-dfa5f605cc6c
is_awesome = true

# ╔═╡ 6cf7f500-63d7-11eb-285a-4db31a47e386
typeof(is_awesome)

# ╔═╡ 91be39d0-63dc-11eb-3275-f17059662b87
md"NB: `true` and `false` are written in lower case, and are reserved keywords, meaning you can cannot assign values to them."

# ╔═╡ 934d5790-63d7-11eb-3ff1-49bc42858003
md"### Conversion"

# ╔═╡ c4fdd4e0-63d7-11eb-338c-997af2ea699f
md"Some types can be converted one into another. To do so, you can use `convert()`, or the name of the type to convert into as follows:"

# ╔═╡ 34e61ce5-96ad-4a99-a0e3-57ee579834a3
typeof(cars)

# ╔═╡ d40d9330-63d7-11eb-2ca7-95207f940b2e
convert(Float64,cars)

# ╔═╡ b458b410-63d8-11eb-3860-8dca66bf91c7
md"or:"

# ╔═╡ b3875a50-63d8-11eb-25dc-f751d868901d
Float64(cars)

# ╔═╡ ad93f4e0-63d9-11eb-240b-f3a189a9c94c
md"We can also transform it into a string using `string()`:"

# ╔═╡ b8a905f0-63d9-11eb-221d-b11ce15217e9
cars_string = string(cars)

# ╔═╡ 13b8a362-63da-11eb-203d-89ba55e487d0
md"...and back into an integer using `parse()`:"

# ╔═╡ 1c771a40-63da-11eb-10a0-7b33fdaf574b
parse(Int,cars_string)

# ╔═╡ 9f525c00-63d9-11eb-08f0-c1e0b1dd9b91
md"Some types cannot be converted into another without loss of information. For example what if you convert a float like `2.2` into an integer? "

# ╔═╡ 89ce2162-63da-11eb-149a-8d58c6c324ed
Int(2.2)

# ╔═╡ 904c6c40-63da-11eb-3324-43f814d0e102
md"Julia returns an `InexactError` because there is a loss of information here (the `.2` part). But what about `2.0` then?"

# ╔═╡ b16d9020-63da-11eb-23a2-6b67ea2813c0
Int(2.0)

# ╔═╡ b62b3f40-63da-11eb-2206-951745a869ee
md"It works! Why? Because there is no loss of information here. What you should do in the case of `2.2`? Cut out some information:"

# ╔═╡ 9d60cfc0-63da-11eb-093f-6b7f21c7d5fa
rounded = round(2.2)

# ╔═╡ 1db52310-63db-11eb-0c78-eb033ec18343
md"Now we can convert it into an integer safely:"

# ╔═╡ 15146d5e-63db-11eb-303d-0bc6c062359b
Int(rounded)

# ╔═╡ 3562b0e0-63db-11eb-18bc-514ffe3a0d56
md"There is also a shortcut for that:"

# ╔═╡ 40ca4ec0-63db-11eb-0aec-db92ee1848d2
round(Int,2.2)

# ╔═╡ 4602d330-63db-11eb-25f3-69cf04127264
md"You can also use other functions such as `floor()`, `ceil()` or `trunk()`."

# ╔═╡ 10c15500-63dd-11eb-17dc-67783d644ca2
md"### Promotion"

# ╔═╡ 18016940-63dd-11eb-1232-f15f2ce37d58
md"""Promotion is converting multiple values with different types into a single common type with the least loss of information. Not all types have a promotion rule, but Real numbers do:"""

# ╔═╡ 41b008f0-63dd-11eb-0975-7fbab120ad8b
promote(1.0,1)

# ╔═╡ 8c1e4b90-63dd-11eb-34a3-0d118a0eb642
md"Here `promote()` returned both arguments as floating point numbers."

# ╔═╡ Cell order:
# ╟─e5b53410-63cf-11eb-2532-358a68b59eb9
# ╟─16f1e0f0-63d0-11eb-2ee4-8f8be24866e8
# ╟─31972af0-63d0-11eb-3e22-bb3370a9e40d
# ╠═9d83cca0-63d0-11eb-0346-751060e8cab5
# ╠═ffed727d-b072-4797-9667-dc3a1c7e0109
# ╟─ba8fdd70-63d0-11eb-0565-f3e77f3f802f
# ╠═a8247100-63d0-11eb-0442-11949579c4ba
# ╟─f1641eb0-63d0-11eb-2ed0-d5ab1ce695ae
# ╠═e20a69a0-63d1-11eb-3f1e-159164d5e71d
# ╟─2a77d330-63d2-11eb-16c3-d3be6a6d9975
# ╠═f85e81f0-63d1-11eb-2385-0196982bc487
# ╟─8ab14f10-63d2-11eb-3874-073f07427a26
# ╠═257727a2-63d2-11eb-311b-0b14f3d0f747
# ╟─b6873b92-63d2-11eb-1a56-1d01d825cf01
# ╠═e76d15e0-63d2-11eb-111e-8b0637007258
# ╠═fa35e490-63d2-11eb-210d-0b44aeb5cbfa
# ╠═68fb21b0-63d3-11eb-1792-45baa321e122
# ╟─6f7a56f0-63d3-11eb-113b-b11ac60caaed
# ╟─8813c0d0-63d2-11eb-1723-3b802f8a9f54
# ╟─a067ae80-63d7-11eb-3d1e-39eae17d3434
# ╟─f4b09520-63d1-11eb-25d9-7564a91fa7d5
# ╠═ba6332b0-63d1-11eb-17a5-0160686e19e7
# ╟─6fa8ab30-63d9-11eb-0ce9-0f64d74d9b96
# ╠═fe28579e-63d6-11eb-31c6-59b910298cc2
# ╟─ccd15a10-63d8-11eb-2cce-9b5b18f2f003
# ╟─3279b400-63d6-11eb-093a-4539f1254a83
# ╠═3d9d6b0e-63d6-11eb-2cc1-39a76f290513
# ╠═0d5dc660-63d7-11eb-0458-29b554125fce
# ╟─53e9e240-63d6-11eb-00dc-098c64e573f1
# ╠═61ffb030-63d6-11eb-059d-e912b2c828fb
# ╠═1204e040-63d7-11eb-2c61-e3f212329a44
# ╟─72e5b520-63d6-11eb-0463-3bd1d3643d15
# ╠═7b9537e0-63d6-11eb-1d93-d3e253614910
# ╠═19752c40-63d7-11eb-0500-11c3318c9319
# ╟─48cab050-63dc-11eb-2f06-651fa9acbd54
# ╠═69283610-63dc-11eb-1d97-25d73d0ce439
# ╟─27efd680-63d7-11eb-0701-67d7c4fc2847
# ╟─5c08ef5e-63d7-11eb-34d0-dfa5f605cc6c
# ╠═6cf7f500-63d7-11eb-285a-4db31a47e386
# ╟─91be39d0-63dc-11eb-3275-f17059662b87
# ╟─934d5790-63d7-11eb-3ff1-49bc42858003
# ╟─c4fdd4e0-63d7-11eb-338c-997af2ea699f
# ╠═34e61ce5-96ad-4a99-a0e3-57ee579834a3
# ╠═d40d9330-63d7-11eb-2ca7-95207f940b2e
# ╟─b458b410-63d8-11eb-3860-8dca66bf91c7
# ╠═b3875a50-63d8-11eb-25dc-f751d868901d
# ╟─ad93f4e0-63d9-11eb-240b-f3a189a9c94c
# ╠═b8a905f0-63d9-11eb-221d-b11ce15217e9
# ╟─13b8a362-63da-11eb-203d-89ba55e487d0
# ╠═1c771a40-63da-11eb-10a0-7b33fdaf574b
# ╟─9f525c00-63d9-11eb-08f0-c1e0b1dd9b91
# ╠═89ce2162-63da-11eb-149a-8d58c6c324ed
# ╟─904c6c40-63da-11eb-3324-43f814d0e102
# ╠═b16d9020-63da-11eb-23a2-6b67ea2813c0
# ╟─b62b3f40-63da-11eb-2206-951745a869ee
# ╠═9d60cfc0-63da-11eb-093f-6b7f21c7d5fa
# ╟─1db52310-63db-11eb-0c78-eb033ec18343
# ╠═15146d5e-63db-11eb-303d-0bc6c062359b
# ╟─3562b0e0-63db-11eb-18bc-514ffe3a0d56
# ╠═40ca4ec0-63db-11eb-0aec-db92ee1848d2
# ╟─4602d330-63db-11eb-25f3-69cf04127264
# ╟─10c15500-63dd-11eb-17dc-67783d644ca2
# ╟─18016940-63dd-11eb-1232-f15f2ce37d58
# ╠═41b008f0-63dd-11eb-0975-7fbab120ad8b
# ╟─8c1e4b90-63dd-11eb-34a3-0d118a0eb642


text:
module OpenSpiel_jll

include("$(@__DIR__)/../deps/deps.jl")

using CxxWrap
import CxxWrap:argument_overloads
import Base:show, length, getindex, setindex!, keys, values, copy, deepcopy, first, last, step, getfield, setfield!

@wrapmodule(LIB_OPEN_SPIEL)

include("patch.jl")

# export all
for n in names(@__MODULE__(); all=true)
    if Base.isidentifier(n) &&
        !startswith(String(n), "_") &&
        n ∉ (Symbol(@__MODULE__()), :eval, :include)
        @eval export $n
    end
end

function __init__()
    @initcxx
end

end


text:
using SparseArrays
using ViewerGL
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

V = [0.1972649 0.2464999 0.2187192 0.2679542 0.686082 0.5932705 0.8533899 0.7605785 1.1832645 1.3287132 1.4321643 1.5776131 1.0224958 1.2029449 0.9847223 1.1651714 1.3227588 1.353673 0.9990077 1.0299219 1.5318421 1.4893674 1.2650475 1.2225729 0.7126548 0.7569095 0.6814385 0.7256931 0.4524947 0.3425744 0.5812341 0.4713138 1.0960876 1.262034 0.9913078 1.1572542 0.9618413 0.6418082 1.0018243 0.6817912; 1.78830443 1.80930533 1.73800643 1.75900733 1.48466913 1.53732743 1.77955313 1.83221143 1.86718673 2.18004573 1.75147273 2.06433163 1.94003713 1.91927233 1.61178063 1.59101573 1.43677903 1.35742863 1.31064863 1.23129813 1.89724883 1.68246173 1.95000803 1.73522103 1.64723193 1.61188753 1.60814613 1.57280163 1.40879763 1.47291863 1.62949063 1.69361163 2.32095143 2.11285603 2.23739433 2.02929903 1.19415203 1.24807003 1.43147343 1.48539143]
EV = Array{Int64,1}[[1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8], [9, 10], [11, 12], [9, 11], [10, 12], [13, 14], [15, 16], [13, 15], [14, 16], [17, 18], [19, 20], [17, 19], [18, 20], [21, 22], [23, 24], [21, 23], [22, 24], [25, 26], [27, 28], [25, 27], [26, 28], [29, 30], [31, 32], [29, 31], [30, 32], [33, 34], [35, 36], [33, 35], [34, 36], [37, 38], [39, 40], [37, 39], [38, 40]]

VV = [[k] for k=1:size(V,2)]
GL.VIEW( GL.numbering(.05)((V,[VV, EV]),GL.COLORS[1]) );

# arrangement of input data
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV)
V, copEV, copFE = Lar.planar_arrangement(W, cop_EV)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# compute euler characteristic
χ = Lar.euler_characteristic(V, copEV, copFE);
println("χ = $χ ; bicon_comps = $(length(bicon_comps))");

# visualization of component graphs
EW = Lar.cop2lar(copEV)
W = convert(Lar.Points, V')
comps = [ GL.GLLines(W,EW[comp],GL.COLORS[(k-1)%12+1]) for (k,comp) in enumerate(bicon_comps) ];
GL.VIEW(comps);

# visualization of numbered arrangement
VV = [[k] for k=1:size(V,1)]
EV = Lar.cop2lar(copEV);
FE = Lar.cop2lar(copFE);
# FV = cat([[EV[e] for e in face] for face in FE])

# final solid visualization
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,.6));

# polygonal face boundaries
EVs = Lar.FV2EVs(copEV, copFE)
EVs = convert(Array{Array{Array{Int64,1},1},1}, EVs)
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,99,1));


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    tess_set_debug_param(
        inst::TessInst,
        name::AbstractString,
        value::Integer
    )::Bool

Sets an debug integer variable in th Tesseract engine.  Returns `false` if the parameter
was not found.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | inst  |         | The instance to set the variable in
| R | name  |         | The name of the variable to set.
| R | value |         | The value to set.

__Details:__

If the parameter is not a debug setting then the value will not be changed.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_set_debug_param(instance, "classify_learning_debug_level", 0)
true
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_get_param`](@ref)
"""
function tess_set_debug_param(
            inst::TessInst,
            name::AbstractString,
            value::Integer
        )::Bool

    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    local result = ccall(
        (:TessBaseAPISetDebugVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Cstring),
        inst,
        name,
        string(value)
    )

    return result == 1
end

# =========================================================================================
"""
    tess_set_debug_param(
        inst::TessInst,
        name::AbstractString,
        value::Bool
    )::Bool

Sets a debug boolean variable in th Tesseract engine.  Returns `false` if the parameter was
not found.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | inst  |         | The instance to set the variable in
| R | name  |         | The name of the variable to set.
| R | value |         | The value to set.

__Details:__

If the parameter is not a debug setting then the value will not be changed.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_set_debug_param(instance, "textord_debug_tabfind", false)
true
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_get_param`](@ref)
"""
function tess_set_debug_param(
            inst::TessInst,
            name::AbstractString,
            value::Bool
        )::Bool

    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    local result = ccall(
        (:TessBaseAPISetDebugVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Cstring),
        inst,
        name,
        value ? "true" : "false"
    )

    return result == 1
end

# =========================================================================================
"""
    tess_set_debug_param(
        inst::TessInst,
        name::AbstractString,
        value::Float64
    )::Bool

Sets a debug `Float64` variable in th Tesseract engine.  Returns `false` if the parameter
was not found.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | inst  |         | The instance to set the variable in
| R | name  |         | The name of the variable to set.
| R | value |         | The value to set.

__Details:__

This method is implmented for future enchancements.  Currently there are no tesseract debug
parameters.

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_get_param`](@ref)
"""
function tess_set_debug_param(
            inst::TessInst,
            name::AbstractString,
            value::Float64
        )::Bool

    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    local result = ccall(
        (:TessBaseAPISetDebugVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Cstring),
        inst,
        name,
        string(value)
    )

    return result == 1
end

# =========================================================================================
"""
    tess_set_debug_param(
        inst::TessInst,
        name::AbstractString,
        value::AbstractString
    )::Bool

Sets a debug string variable in th Tesseract engine.  Returns `false` if the parameter was
not found.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | inst  |         | The instance to set the variable in
| R | name  |         | The name of the variable to set.
| R | value |         | The value to set.

__Details:__

If the parameter is not a debug setting then the value will not be changed.

__Examples:__

```jldoctest
julia> using Tesseract

julia> download_languages()
true

julia> instance = TessInst()
Allocated Tesseract instance.

julia> tess_set_debug_param(instance, "debug_file", "debug.log")
true
```

See also: [`tess_params`](@ref), [`tess_params_parsed`](@ref), [`tess_get_param`](@ref)
"""
function tess_set_debug_param(
            inst::TessInst,
            name::AbstractString,
            value::AbstractString
        )::Bool

    if is_valid(inst) == false
        @error "Instance has been freed."
        return false
    end

    local result = ccall(
        (:TessBaseAPISetDebugVariable, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Cstring),
        inst,
        name,
        value
    )

    return result == 1
end


text:
# Von-Mises Fisher: a multivariate distribution useful in directional statistics

# Useful notes:
# http://www.mitsuba-renderer.org/~wenzel/vmf.pdf
# Some of the code adapted from http://www.unc.edu/~sungkyu/manifolds/randvonMisesFisherm.m
# as well as the movMF R package.

immutable VonMisesFisher <: ContinuousMultivariateDistribution
    mu::Vector{Float64}
    kappa::Float64

    function VonMisesFisher{T <: Real}(mu::Vector{T}, kappa::Float64)
        mu = mu ./ norm(mu)
        if kappa < 0
            throw(ArgumentError("kappa must be a nonnegative real number."))
        end
        new(float64(mu), kappa)
    end
end

length(d::VonMisesFisher) = length(d.mu)
mean(d::VonMisesFisher) = d.mu
scale(d::VonMisesFisher) = d.kappa

insupport{T<:Real}(d::VonMisesFisher, x::AbstractVector{T}) = abs(sum(x) - 1.) < 1e-8

function _logpdf{T<:Real}(d::VonMisesFisher, x::DenseVector{T}; stable=true)
    if abs(d.kappa - 0.0) < eps()
        return 0.25 / pi
    end
    if stable
        # As suggested by Wenzel Jakob: http://www.mitsuba-renderer.org/~wenzel/vmf.pdf
        return d.kappa * dot(d.mu, x) - d.kappa + log(d.kappa) - log(2*pi) - log(1-exp(-2*d.kappa))
    else
        # As described on Wikipedia
        p = length(d)
        logCpk = 0.0
        if p == 3
            logCpk = log(d.kappa) - log(2 * pi * (exp(kappa) - exp(-kappa)))
        else
            logCpk = (p/2 - 1) * log(d.kappa) - (p/2) * log(2*pi) - log(besselj(p/2-1, d.kappa))
        end
        return d.kappa * dot(d.mu, x) + logCpk
    end
end

# sampling (TODO: make it consistent with the common API)

function rand(d::VonMisesFisher, n::Int)
    randvonMisesFisher(n, d.kappa, d.mu)
end

function randvonMisesFisher(n, kappa, mu)
    m = length(mu)
    w = rW(n, kappa, m)
    v = rand(MvNormal(zeros(m-1), eye(m-1)), n)

    # normalize each column of v
    for j = 1:n
        s = 0.
        vj = view(v,:,j)
        for i = 1:size(v,1)
            s += abs2(vj[i])
        end
        s = sqrt(s)
        for i = 1:size(v,1)
            vj[i] /= s
        end
    end
    v = v'

    r = sqrt(1.0 .- w .^ 2)
    for j = 1:size(v,2) v[:,j] = v[:,j] .* r; end
    x = hcat(v, w)
    mu = mu / norm(mu)
    return rotMat(mu)'*x'
end

# Randomly sample W
function rW(n, kappa, m)
    y = zeros(n)
    l = kappa;
    d = m - 1;
    b = (- 2. * l + sqrt(4. * l * l + d * d)) / d;
    x = (1. - b) / (1. + b);
    c = l * x + d * log(1. - x * x);
    w = 0
    for i=1:n
        done = false
        while !done
            z = rand(Beta(d / 2., d / 2.))
            w = (1. - (1. + b) * z) / (1. - (1. - b) * z);
            u = rand()
            if l * w + d * log(1. - x * w) - c >= log(u)
                done = true
            end
        end
        y[i] = w
    end
    return y
end

# Rotation helper function
function rotMat(b)
    d = length(b)
    b= b/norm(b)
    a = [zeros(d-1,1); 1]
    alpha = acos(a'*b)[1]
    c = b - a * (a'*b); c = c / norm(c)
    A = a*c' - c*a'
    return eye(d) + sin(alpha)*A + (cos(alpha) - 1)*(a*a' +c*c')
end

# Each row of x assumed to be ~ VonMisesFisher(mu, kappa)
# MLE notes from: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.186.1887&rep=rep1&type=pdf

function fit_mle(::Type{VonMisesFisher}, x::Matrix{Float64})
    (n,p) = size(x)
    sx = sum(x, 1)
    mu = sx[:] / norm(sx)
    rbar = norm(sx) / n
    kappa0 = rbar * (p-rbar^2) / (1-rbar^2)  # Eqn. 4
    # TODO: Include a few Newton steps to get a better approximation.
    # A(p,kappa) = besselj(p/2, kappa) / besselj(p/2-1, kappa)
    # apk0 = A(p,kappa0)
    # kappa1 = kappa0 + (apk0 - rbar) / (1 - apk0^2 - (p-1)*apk0/kappa0)
    # apk1 = A(p,kappa1)
    # kappa2 = kappa1 + (apk1 - rbar) / (1 - apk1^2 - (p-1)*apk1/kappa1)
    return VonMisesFisher(mu, kappa0)#, kappa1, kappa2)
end


text:
# # Social networks with Graphs.jl

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../schoolyard.mp4" type="video/mp4">
# </video>
# ```

# Many ABM frameworks provide graph infrastructure for analysing network properties of agents.
# Agents.jl is no different in that aspect, we have [`GraphSpace`](@ref) for when spatial structure
# is not important, but connections are.

# What if you wish to model something a little more complex? Perhaps a school yard full of students
# running around (in space), interacting via some social network. This is precisely the scenario that
# the [MASON](https://cs.gmu.edu/~eclab/projects/mason/) ABM framework uses as an introductory example
# in their [documentation](https://cs.gmu.edu/~eclab/projects/mason/manual.pdf).

# Rather than implementing an Agents.jl⸺specific graph structure, we can interface with
# [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl): a high class library for managing
# and implementing graphs, which can be re-used to establish social networks within existing spaces.

# To begin, we load in some dependencies

using Agents
using SimpleWeightedGraphs: SimpleWeightedDiGraph # will make social network
using SparseArrays: findnz                        # for social network connections
using Random: MersenneTwister                     # reproducibility

# And create a very simple agent without any extra properties
mutable struct Student <: AbstractAgent
    id::Int
    pos::Tuple{Float64,Float64}
end

# ## Rules of the schoolyard

# It's lunchtime, and the students are going out to play.
# We assume the school building is in the centre of our space, with some fences around the building.
# A teacher monitors the students, and makes sure they don't stray too far towards the fence.
# We use a `teacher_attractor` force to simulate a teacher's attentiveness.
# Students head out to the schoolyard in random directions, but adhere to some social norms.

# Each student has one *friend* and one *foe*. These are chosen at random in our model, so it's
# possible that for any pair of students, one likes the other but this feeling is not reciprocated.
# The bond between pairs is chosen at random between 0 and 1, with a bond of 1 being the strongest.
# If the bond is *friendly*, agents wish above all else to be near their *friend*.
# Bonds that are *unfriendly* see students moving as far away as possible from their *foe*.

# ## Initialising the model

function schoolyard(;
    numStudents = 50,
    teacher_attractor = 0.15,
    noise = 0.1,
    max_force = 1.7,
    spacing = 4.0,
    seed = 6998,
)
    model = ABM(
        Student,
        ContinuousSpace((100, 100), spacing; periodic = false);
        properties = Dict(
            :teacher_attractor => teacher_attractor,
            :noise => noise,
            :buddies => SimpleWeightedDiGraph(numStudents),
            :max_force => max_force,
        ),
        rng = MersenneTwister(seed)
    )
    for student in 1:numStudents
        ## Students begin near the school building
        add_agent!(model.space.extent .* 0.5 .+ Tuple(rand(model.rng, 2)) .- 0.5, model)

        ## Add one friend and one foe to the social network
        friend = rand(model.rng, filter(s -> s != student, 1:numStudents))
        add_edge!(model.buddies, student, friend, rand(model.rng))
        foe = rand(model.rng, filter(s -> s != student, 1:numStudents))
        add_edge!(model.buddies, student, foe, -rand(model.rng))
    end
    model
end

# Our model contains the `buddies` property, which is our Graphs.jl directed, weighted graph.
# As we can see in the loop, we choose one `friend` and one `foe` at random for each `student` and
# assign their relationship as a weighted edge on the graph.

# ## Movement dynamics

distance(pos) = sqrt(pos[1]^2 + pos[2]^2)
scale(L, force) = (L / distance(force)) .* force

function agent_step!(student, model)
    ## place a teacher in the center of the yard, so we don’t go too far away
    teacher = (model.space.extent .* 0.5 .- student.pos) .* model.teacher_attractor

    ## add a bit of randomness
    noise = model.noise .* (Tuple(rand(model.rng, 2)) .- 0.5)

    ## Adhere to the social network
    network = model.buddies.weights[student.id, :]
    tidxs, tweights = findnz(network)
    network_force = (0.0, 0.0)
    for (widx, tidx) in enumerate(tidxs)
        buddiness = tweights[widx]
        force = (student.pos .- model[tidx].pos) .* buddiness
        if buddiness >= 0
            ## The further I am from them, the more I want to go to them
            if distance(force) > model.max_force # I'm far enough away
                force = scale(model.max_force, force)
            end
        else
            ## The further I am away from them, the better
            if distance(force) > model.max_force # I'm far enough away
                force = (0.0, 0.0)
            else
                L = model.max_force - distance(force)
                force = scale(L, force)
            end
        end
        network_force = network_force .+ force
    end

    ## Add all forces together to assign the students next position
    new_pos = student.pos .+ noise .+ teacher .+ network_force
    move_agent!(student, new_pos, model)
end

# Applying the rules for movement is relatively simple. For the network specifically,
# we find the student's `network` and figure out how far apart they are. We scale this
# by the `buddiness` factor (how much force we should apply), then figure out if
# that force should be in a positive or negative direction (*friend* or *foe*?).

# The `findnz` function is something that may require some further explanation.
# Graphs uses sparse vectors internally to efficiently represent data.
# When we find the `network` of our `student`, we want to convert the result to
# a dense representation by **find**ing the **n**on-**z**ero (`findnz`) elements.

# ## Visualising the system

# Now, we can watch the dynamics of the social system unfold:
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

model = schoolyard()

function static_preplot!(ax, model)
    obj = CairoMakie.scatter!([50 50]; color = :red) # Show position of teacher
    CairoMakie.hidedecorations!(ax) # hide tick labels etc.
    CairoMakie.translate!(obj, 0, 0, 5) # be sure that the teacher will be above students
end

abm_video(
    "schoolyard.mp4", model, agent_step!, dummystep;
    framerate = 15, frames = 40,
    title = "Playgound dynamics",
    static_preplot!,
)

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../schoolyard.mp4" type="video/mp4">
# </video>
# ```