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

using Markdown
using InteractiveUtils

# ╔═╡ 5a98d3d0-1d31-11eb-0caf-5dca48313597
using Pkg, DrWatson

# ╔═╡ 61f78d42-1d31-11eb-31f1-63eb86dabf9f
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics
end

# ╔═╡ 29db8358-1d24-11eb-166a-e1b1c886af45
md"## Listing 1.2"

# ╔═╡ af985f1c-1d23-11eb-335b-c3abd97e66ea
begin
	array1 = [(2n+1)^2 for n in 1:5]
	array2 = [sqrt(i) for i in array1]
	typeof(1:5), typeof(array1), typeof(array2)
end

# ╔═╡ d4e0542a-1d23-11eb-1ef2-3d3b86cd30cb
1:5, array1, array2

# ╔═╡ fe1c39ce-1d23-11eb-0790-7d996206bd99
md"##### Open `Live docs` and start typing or click on `sqrt`."

# ╔═╡ d3cdd15a-1d23-11eb-3e0c-91a6135d3407
sqrt

# ╔═╡ 06a876ce-1d2f-11eb-2128-e1e8c67779cb
md"## Listing 1.3"

# ╔═╡ 293b0710-1d2f-11eb-2da1-2d516ce4c85f
begin
	Text(@time begin
			data1 = Float64[]
				for _ in 1:10^6
					group = Float64[]
						for _ in 1:5*10^2
							push!(group,rand())
						end
					push!(data1,mean(group))
				end
			"98% of the means lie in the estimated range: ",
                        (quantile(data1,0.01), quantile(data1,0.99))
		end
	)
end

# ╔═╡ e1c7f0e8-1d3b-11eb-2545-53bb0fe2569a
md"##### Note that the result of the `@time ...` macro is written to the console (where you started Pluto)."

# ╔═╡ 0ade9930-1d2f-11eb-2329-334dc91bc5a9
md"## Listing 1.4"

# ╔═╡ b283773a-1d3b-11eb-21b8-89af728c6fed
@time begin
	data2 = [mean(rand(5*10^2)) for _ in 1:10^6]
        "98% of the means lie in the estimated range: ",
                        (quantile(data2,0.01),quantile(data2,0.99))
end

# ╔═╡ a95250f0-1d3b-11eb-0492-5576bc284913
md"## Listing 1.5"

# ╔═╡ f2bc632e-1d3c-11eb-0fab-b741d7d2d04f
begin
	data3 = [1,2,3]
	s = 0
	beta, gamma = 2, 1

	for i in 1:length(data3)
		global s += beta*data3[i]
		data3[i] *= -1
	end

	"Sum of data in external scope: ", s, data3
end

# ╔═╡ 36795a2c-1d3d-11eb-22fe-7f694a7a816c
begin

	function sumData(beta)
		s = 0
		for i in 1:length(data3)
			s += data3[i] + gamma
		end
		return s
	end

	"Sum of data in a function: ", s, sumData(beta/2)
end

# ╔═╡ e0d70c94-1d3d-11eb-1caf-4dc5c83c9482
md"##### As mentioned before, `begin` and `end` is a very important construct in Pluto notebooks!"

# ╔═╡ 0a6f9b48-1d3e-11eb-336f-fd577cf395ef
begin
	"This will be printed!"
end

# ╔═╡ 338544f6-1d3e-11eb-1ca7-c994c82458f5
md"##### Suppress printing of the final result by a `;` after `end`."

# ╔═╡ 20635ec6-1d3e-11eb-3ad1-b99b1fe6b4c2
begin
	"This will not be printed!"
end;

# ╔═╡ Cell order:
# ╠═29db8358-1d24-11eb-166a-e1b1c886af45
# ╠═5a98d3d0-1d31-11eb-0caf-5dca48313597
# ╠═61f78d42-1d31-11eb-31f1-63eb86dabf9f
# ╠═af985f1c-1d23-11eb-335b-c3abd97e66ea
# ╠═d4e0542a-1d23-11eb-1ef2-3d3b86cd30cb
# ╟─fe1c39ce-1d23-11eb-0790-7d996206bd99
# ╠═d3cdd15a-1d23-11eb-3e0c-91a6135d3407
# ╟─06a876ce-1d2f-11eb-2128-e1e8c67779cb
# ╠═293b0710-1d2f-11eb-2da1-2d516ce4c85f
# ╟─e1c7f0e8-1d3b-11eb-2545-53bb0fe2569a
# ╟─0ade9930-1d2f-11eb-2329-334dc91bc5a9
# ╠═b283773a-1d3b-11eb-21b8-89af728c6fed
# ╟─a95250f0-1d3b-11eb-0492-5576bc284913
# ╠═f2bc632e-1d3c-11eb-0fab-b741d7d2d04f
# ╠═36795a2c-1d3d-11eb-22fe-7f694a7a816c
# ╟─e0d70c94-1d3d-11eb-1caf-4dc5c83c9482
# ╠═0a6f9b48-1d3e-11eb-336f-fd577cf395ef
# ╟─338544f6-1d3e-11eb-1ca7-c994c82458f5
# ╠═20635ec6-1d3e-11eb-3ad1-b99b1fe6b4c2


text:
import TPT

using Base.Test

using DataFrames
using Optim
using Plots; pyplot()

println("--- TM Energy ---")

# Read elemental parameters (while checking if we are in the correct place)
p = readtable(joinpath("data", "parameters", "tm_optim.csv"))

#
# Prepare the directory to store the results
#
!isdir("results") && mkdir("results")

resdir = joinpath("results", "tm_energy")
!isdir(resdir) && mkdir(resdir)

# Physical constants
kB = TPT.kB

# number of elements available
M = size(p, 1)

# Various thermodynamic quantities
enr = [:F, :U, :U_nfe, :U_tb, :U_bond, :U_es, :S, :S_gas, :S_conf, :S_ref, :S_nfe, :S_tb]

for B in 3:M
  A = 5
  # (A,B) = ind2sub((M,M), k)

  X₁ = p[:X][A]
  X₂ = p[:X][B]

  println("$(X₁)-$(X₂)")

  T = float(p[:T][A])

  for x in [:σ, :ρ, :m, :zs, :rc, :a, :zd, :rd, :Ed, :Wd, :R0]
    v = [ p[x][A], p[x][B] ]
    @eval ($x) = $v
  end

  for E in enr
    ΔE = Symbol(:Δ, E)
    @eval ($E) = Vector{Float64}(11)
    @eval ($ΔE) = Vector{Float64}(11)
  end

  res = Vector{TPT.TPTSystem}(11)

  Threads.@threads for i in 1:11
    x₂ = (i-1) / 10

    ρ₀ = (1-x₂)*ρ[1] + x₂*ρ[2]
    c = [1-x₂, x₂]

    ahs = TPT.AHS(ρ = ρ₀, σ = σ, c = c, approx = :RFA)
    wca = TPT.LWCA(ahs, T, struct = :full)

    pse = TPT.BretonnetSilbert(zs, rc, a)
    nfe = TPT.NFE(ahs, pse)
    tb = TPT.HHTB(c, zd, Ed, Wd, R0)
    nfetb = TPT.NFETB(nfe, tb)

    sys = TPT.TPTSystem(wca, nfetb, m = m)
    res[i] = sys

    U_nfe[i] = TPT.internal_pair(sys.pert.nfe, sys.ref)
    U_tb[i] = TPT.internal(sys.pert.tb, sys.ref)
    U_bond[i] = TPT.internal(sys.ref, tb, TPT.pairpotential_bond(tb))
    U_es[i] = TPT.internal_es(sys.pert.nfe, sys.ref)

    U[i] = U_nfe[i] + U_tb[i] + U_es[i]

    S_gas[i] = TPT.entropy_gas(sys)
    S_conf[i] = TPT.entropy_conf(sys)
    # S_ref[i] = TPT.entropy(sys.ref)
    S_ref[i] = 0.
    S_nfe[i] = TPT.entropy(nfetb.nfe, sys.ref, T)
    # S_tb[i] = TPT.entropy(nfetb.tb, sys.ref, T)
    S_tb[i] = 0.

    S[i] = S_gas[i] + S_conf[i] + S_ref[i] + S_nfe[i] + S_tb[i]

    K = TPT.kinetic(sys) # this can be omitted

    F[i] = K + U[i] - kB*T*S[i]
  end

  for i in 1:11
    x₂ = (i-1) / 10
    for E in enr
      ΔE = Symbol(:Δ, E)
      @eval ($ΔE)[$i] = ($E)[$i] - (1-$x₂)*($E)[1] - $x₂*($E)[11]
    end
  end

  #
  # Plot various quantities of mixing
  #
  # settings for plots
  default(xlabel="x($X₂)", ylabel="Energy of mixing (kJ/mol)", ylims=())

  plot() # initialize the plot

  # free energy (total) and internal energy
  for E in [:F, :U, :U_nfe, :U_tb, :U_bond, :U_es]
    ΔE = Symbol(:Δ, E)
    @eval plot!(0:0.1:1.0, 2625.5*($ΔE), label=$(string(E)))
  end

  # entropic contribution
  for E in [:S, :S_gas, :S_conf, :S_ref, :S_nfe, :S_tb]
    ΔE = Symbol(:Δ, E)
    @eval plot!(0:0.1:1.0, 2625.5*$(T*kB)*($ΔE), label=$(string("T", E)))
  end

  # Save the plot as png
  png(joinpath(resdir, "$(A)-$(B)_$(X₁)-$(X₂)_F"))

  # Effective hard-sphere diameter
  σ₁ = [ TPT.hsdiameter(res[i].ref)[1,1] for i in 1:11 ]
  σ₂ = [ TPT.hsdiameter(res[i].ref)[2,2] for i in 1:11 ]
  σmin = min(minimum(σ₁), minimum(σ₂))
  σmax = max(maximum(σ₁), maximum(σ₂))
  default(ylabel = "Effective HS diameter (a.u.)", ylims=(0.99σmin, 1.01σmax))
  plot(0:0.1:1, σ₁, label=X₁)
  plot!(0:0.1:1, σ₂, label=X₂)
  png(joinpath(resdir, "$(A)-$(B)_$(X₁)-$(X₂)_σ"))
end

@testset "TM Energy" begin
end

ext: jl


text:
#=
Utilities for testing QuantEcon

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-26
=#
using HDF5, JLD, MAT

const test_path = dirname(@__FILE__())
const data_path = joinpath(test_path, "data")

const quad_data_file_name = joinpath(data_path, "matlab_quad.mat")
const ml_quad_data_url = "https://github.com/spencerlyon2/QuantEcon.jl/releases/download/v0.0.1/matlab_quad.mat"

if !(isfile(quad_data_file_name))
    try
        download(ml_quad_data_url, quad_data_file_name)
    catch
        m = """
        Could not download data for quad tests. They will not run properly
        right now. Try again when you have an internet connection
        """

        warn(m)
    end
end


text:
module DenoiseDNA
using FASTX, BioSequences,  BioAlignments # BioJulia
import BioAlignments.pairalign
using StatsFuns, Loess,LazyStack, StatsBase,LoopVectorization, CodecZlib,DataFrames
using CairoMakie,UnicodePlots, Poppler_jll # visualisation
using RCall

include.(["plotting.jl",
        "fileIO.jl",
        "kmers.jl",
        "DADA.jl"])

# Write your package code here.
export quality_plot,
        dada,
        custering,
        trim,
        dereplicate,
        get_files,
        learnErrors
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "USRPHardwareDriver"
version = v"4.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/EttusResearch/uhd.git", "6bd0be9cda5db97081e4f3ee3127c45eed21239c")
]

dependencies = [
    Dependency("boost_jll"; compat="=1.76.0"), # max gcc7
    Dependency("libusb_jll", compat="1.0.24"),
]

# Bash recipe for building across all platforms
script = raw"""
cd uhd/host

#code generators, only needed at build time
apk add py3-mako py3-ruamel.yaml

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_EXAMPLES=OFF \
    -DENABLE_TESTS=OFF \
    -DENABLE_PYTHON_API=OFF \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# TODO: Windows has several issues with boost threads. There is a WIP branch:
# https://github.com/JuliaTelecom/uhd/tree/juliatelecom/patch-v4.1.0.1
platforms = expand_cxxstring_abis(filter!(p -> !Sys.iswindows(p) && !in(arch(p),("armv7l","armv6l")), supported_platforms(;experimental=true)))

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("uhd_cal_rx_iq_balance", :uhd_cal_rx_iq_balance),
    ExecutableProduct("uhd_cal_tx_iq_balance", :uhd_cal_tx_iq_balance),
    ExecutableProduct("uhd_find_devices", :uhd_find_devices),
    ExecutableProduct("uhd_adc_self_cal", :uhd_adc_self_cal),
    ExecutableProduct("uhd_cal_tx_dc_offset", :uhd_cal_tx_dc_offset),
    ExecutableProduct("uhd_config_info", :uhd_config_info),
    ExecutableProduct("uhd_image_loader", :uhd_image_loader),
    ExecutableProduct("uhd_usrp_probe", :uhd_usrp_probe),
    LibraryProduct("libuhd", :libuhd)
]

# Build the tarballs, and possibly a `build.jl` as well.
# gcc7 constraint from boost
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module util
using Base.Threads: @threads
using ..Formatting, ..Distributions, ..StatsBase
using Sobol
using Random: GLOBAL_RNG, MersenneTwister

export unsqueeze

################################################################################
##                                                                            ##
##                                Generic                                     ##
##                                                                            ##
################################################################################

e_k(T, n, k) = begin; out = zeros(T, n); out[k] = 1; out; end
e_k(n, k) = e_k(Float32, n, k)
e1(T, n) = e_k(T, n, 1)
e1(n) = e_k(Float32, n, 1)

unsqueeze(xs, dim) = reshape(xs, (size(xs)[1:dim-1]..., 1, size(xs)[dim:end]...));
vflatten(x) = reduce(vcat, x)
collapsedims1_2(X::AbstractArray) = reshape(X, prod(size(X)[1:2]), size(X)[3])
"""
    get_strtype_wo_params(x)

Useful for highly-parameterized models. The print (or typeof) for these structs
result in the overall type, plus a load of parameters, or details in
parentheses. This function just returns the first part: the overall type,
without all the other junk. This uses the Base.show methods, but captures the
output via IOBuffer.
"""
function get_strtype_wo_params(x)
    _io = IOBuffer()
    print(_io, x)    # (does not print to stdout)
    tstr = String(take!(_io))
    close(_io)
    ixparen = findfirst("(", tstr)
    ixparen === nothing && return tstr
    return tstr[1:first(ixparen)-1]
end

function sq_diff_matrix(X, Y)
    #=
    Constructs $(x_i - y_j)^T (x_i - y_j)$ matrix where
    X = Array(n_x, d)
    Y = Array(n_y, d)
    return: Array(n_x, n_y)
    =#
    normsq_x = sum(a -> a^2, X, dims=2)
    normsq_y = sum(a -> a^2, Y, dims=2)
    @assert size(normsq_x, 2) == 1 && size(normsq_y, 2) == 1
    out = normsq_x .+ normsq_y'    # outer
    @assert size(out) == (size(normsq_x, 1), size(normsq_y, 1))
    out .-= 2*X * Y'
    return out
end

sq_diff_matrix(X::MaskedArray, Y) = _sq_diff_matrix_masked(X, Y)
sq_diff_matrix(X, Y::MaskedArray) = _sq_diff_matrix_masked(X, Y)
sq_diff_matrix(X::MaskedArray, Y::MaskedArray) = _sq_diff_matrix_masked(X, Y)
function _sq_diff_matrix_masked(X, Y)
    # This is over an order of magnitude slower than the above
    n_x, n_y = size(X, 1), size(Y, 1)
    [sum(x->x^2, rm_nan(X[i,:]-Y[j,:])) for i in 1:n_x, j in 1:n_y]
end

function _softmax_lse!(out::AbstractVecOrMat{T}, xs::AbstractVecOrMat{T}) where T<:AbstractFloat
    #=
    Simultaneous softmax and logsumexp. Useful for calculating an estimate of logp(x) out of
    lots of MC samples, as well as calculating the approximate posterior.
    We essentially get the Logsumexp for free here.
    (adapted from Mike I's code in NNlib.)
    =#
    lse = Array{T, 1}(undef, size(out, 2))

    @threads for j = 1:size(xs, 2)
        @inbounds begin
            # out[end, :] .= maximum(xs, 1)
            m = xs[1, j]
            for i = 2:size(xs, 1)
                m = max(m, xs[i, j])
            end
            # out .= exp(xs .- out[end, :])
            for i = 1:size(out, 1)
                out[i, j] = exp(xs[i, j] - m)
            end
            # out ./= sum(out, 1)
            s = zero(eltype(out))
            for i = 1:size(out, 1)
                s += out[i, j]
            end
            for i = 1:size(out, 1)
                out[i, j] /= s
            end
            lse[j] = log(s) + m
        end
    end
    return out, lse
end

function softmax_lse!(out, xs; dims=1)
    if dims == 1
        _softmax_lse!(out, xs)
    elseif dims == 2
        xsT, outT = Matrix(xs'), out'  # Matrix = faster to avoid strides for larger matrix
        s, l = _softmax_lse!(outT, xsT)
        s', l
    else
        error("Unsupported dims argument (expecting (1,2)): $dims")
    end
end

softmax_lse!(xs; dims=1) = softmax_lse!(xs, xs; dims=dims)
softmax_lse(xs; dims=1) = softmax_lse!(similar(xs), xs; dims=dims)

################################################################################
##                                                                            ##
##                             Random samples                                 ##
##    These are all copied from my AxUtil package to avoid further imports    ##
##                                                                            ##
################################################################################

sobol_gaussian(n, d) = sobol_gaussian(GLOBAL_RNG, n, d)
"""
    sobol_gaussian([rng::MersenneTwister,] n, d)

Randomized Quasi Monte Carlo samples from a `d` dimensional unit Gaussian. This
function returns `n` such samples in a ``n × d`` matrix. The QMC sequence used
here is the Sobol sequence.
"""
function sobol_gaussian(rng::MersenneTwister, n, d)
    (n == 0) && return zeros(d,0)'
    s = SobolSeq(d)
    p = reduce(hcat, [next!(s) for i = 1:n])
    ϵ = rand(rng, d)
    prand = [(p[j,:] .+ ϵ[j]) .% 1.0 for j in 1:d]
    p = reduce(vcat, [quantile.(Normal(), prand[j])' for j in 1:d])'
    return p
end

randomised_sobol(n, d) = randomised_sobol(GLOBAL_RNG, n, d)
"""
    randomised_sobol([rng::MersenneTwister,] n, d)
Generate `n` samples from `d`-dimensional Sobol sequence, and applying a random
affine mapping to the collection on the [0,1]^d torus.
"""
function randomised_sobol(rng::MersenneTwister, n, d)
    s = SobolSeq(d)
    p = reduce(hcat, [next!(s) for i = 1:n])
    ϵ = rand(rng, d)
    prand = [(p[j,:] .+ ϵ[j]) .% 1.0 for j in 1:d]
    return hcat(prand...)
end

# randomised Sobol within rectangle
uniform_rand_sobol(n, lims...) = uniform_rand_sobol(GLOBAL_RNG, n, lims...)
"""
    uniform_rand_sobol([rng::MersenneTwister,] n, lims)

Randomised sobol sequence (see `randomised_sobol`), mapped to be a quasi random
sample of a uniform distribution on a hyper-cuboid with limits specified by
`lims`. For example:

    uniform_rand_sobol(N, [x_l, x_u], [y_l, y_u])

generates a QMC sample on the rectangle ``[x_l, x_u] × [y_l, y_u]```, with
the first coordinate drawn at random in ``[x_l, x_u]`` and so on.
"""
function uniform_rand_sobol(rng::MersenneTwister, n, lims...)
    d = length(lims)
    rsob = randomised_sobol(rng, n, d)
    for (i, interval) in enumerate(lims)
        @assert (length(interval) == 2) format("interval {:d} does not have length 2", i)
        rsob[:,i] .*= diff(interval)
        rsob[:,i] .+= interval[1]
    end
    return rsob
end

function categorical_sampler(p::AbstractVector, n::Int)
    """
        categorical_sampler(p::AbstractVector, n::Int)

    Categorical sampling using a linear scan strategy, which is surprisingly
    efficient for small `n`. Unlike `rand(Categorical(p), n)`, this function
    accepts Float32 as well as Float64, and is robust in case of `sum(p) != 1`.
    """
    m = length(p)
    x = zeros(Int32, n)

    function linearsearch(p::AbstractArray, m::Int64, rn)
        cs = 0.0
        for ii in 1:m
            @inbounds cs += p[ii]
            if cs > rn
                return ii
            end
        end
        return m
    end

    sump = sum(p)
    for i in 1:n
        rn = rand()*sump
        @inbounds x[i] = linearsearch(p, m, rn)
    end
    return x
end

################################################################################
##                                                                            ##
##                    Preprocessing: StandardScaler                           ##
##                                                                            ##
################################################################################
# This was originally from my Mocap project, but works ~ like Sklearn's version

"""
    fit(MyStandardScaler, X, dims)
Fit a standardisation to a matrix `X` s.t. that the rows/columns have mean zero and standard deviation 1.
This operation calculates the mean and standard deviation and outputs a MyStandardScaler object `s` which
allows this same standardisation to be fit to any matrix using `transform(s, Y)` for some matrix `Y`. Note
that the input matrix `X` is *not* transformed by this operation. Instead use the above `transform` syntax
on `X`.
Note a couple of addendums:
1. Any columns/rows with constant values will result in a standard deviation of 1.0, not 0. This is to avoid NaN errors from the transformation (and it is a natural choice).
2. If only a subset of the rows/columns should be standardised, an additional argument of the indices may be given as:
    `fit(MyStandardScaler, X, operate_on, dims)`
    The subset to operate upon will be maintained through all `transform` and `invert` operations.
"""
mutable struct MyStandardScaler{T}
    μ::AbstractArray{T}
    σ::AbstractArray{T}
    dims
end

Base.copy(s::MyStandardScaler) = MyStandardScaler(copy(s.μ), copy(s.σ), copy(s.operate_on), copy(s.dims))

StatsBase.fit(::Type{MyStandardScaler}, X::AbstractArray, dims=1) =
    MyStandardScaler(mean(MaskedArray(X), dims=dims),
    std(MaskedArray(X), dims=dims), dims) |> post_fit

function post_fit(s::MyStandardScaler)
    bad_ixs = (s.σ .== 0)
    s.σ[bad_ixs] .= 1
    s
end

function scale_transform(s::MyStandardScaler, X::AbstractArray, dims=s.dims)
    (dims != s.dims) && @warn "dim specified in transform is different to specification during fit."
    mu, sigma = s.μ, s.σ
    (X .- mu) ./ sigma
end

function invert(s::MyStandardScaler, X::AbstractArray, dims=s.dims)
    (dims != s.dims) && @warn "dim specified in inversion is different to specification during fit."
    mu, sigma = s.μ, s.σ
    (X .* sigma) .+ mu
end

################################################################################
##                                                                            ##
##               Handling missing (NaN) values: MaskedArray                   ##
##                                                                            ##
################################################################################

"""
    MaskedArray(x::AbstractArray)

I've been dissatisfied for some time about how AD engines etc tend to handle missing values (historically, they basically couldn't, and required a hack... the hack is not terribly hard, but it requires additional admin for coding and data structures etc.) This was an attempt to show that this could be easily handled using a new type of array. But... turns out that a lot more work is needed and I need to understand more about Julia internals to make this work properly (especially Broadcasting), and would need to spend quite some time integrating with Flux's Tracker which is now ~ defunct anyway. So I gave up halfway! Nevertheless, I found the resulting class moderately useful.

The idea is that missing values (represented by NaNs¹) can be handled by a (`.mask`) matrix which determines where the missing values go. The array is then represented with a (`.data`) matrix **only** with numbers in ℝ -- where the missing values might be replaced with any value, but often zero. To recover the original matrix with missing values, the `as_nan` function creates a copy of the `.data` and places the `NaN`s back in their correct place from the mask.

This allows a couple of useful properties for writing AD-enabled code:

1. The missing value mask is carried around *with* the data, rather than requiring separate book-keeping.
1. The AD code can just operate on the `.data`, but then the mask is used to zero out all the gradients which should be ignored.
1. We can natively extract the dimensions of the missing entries e.g. to use as an input mask for a network using `.mask`.
1. Useful versions of `sum`, `mean`, `std` and `var` can be (and have been) defined which ignores the `NaN`s as usually required rather than needing to use some alternative library like NaNMath.

**Broadcasting Issues**:

THIS DOES NOT ALLOW BROADCASTING, and getting to grips with how broadcasting works is necessary to sort this. Simply subtyping from AbstractArray prob doesn't work. The major problem right now is that by default broadcasting iterates over the entire object *one element at a time*. This results in using single indices (scalar values) from the object, which the `.getindex` method implicitly doesn't support :(. To get this to work I'll need to catch when none of the indices are slices and force the result back into an array prior to --> MaskedArray. But even this won't help in the end, as it will result in an Array of 1x1 MaskedArrays. There's no inspiration I can get from TrackedArrays or DualNumbers etc. Utimately I'll need to rtfm here: https://docs.julialang.org/en/latest/manual/interfaces/#man-interfaces-broadcasting-1

... BUT even then I can't natively plug these arrays into Flux, since I'll need to write lots of custom methods to interact with TrackedArrays, and this could easily go wrong. It's almost certainly best just to use MaskedArrays as a useful container and apply the data and mask explicitly during the Flux forward pass... 😂

**Footnote:**

¹I know Julia introduced `missing` with very good support and apparently performance some time ago, but this idea predates this being well supported. I don't know how well it is supported in AD these days. This whole data structure might be entirely unnecessary now.
"""
struct MaskedArray{T, N} <: AbstractArray{T, N}
    data::AbstractArray{T, N}
    mask::AbstractArray{Bool, N}
end
MaskedArray(T::DataType, x::AbstractArray, b::AbstractArray) = MaskedArray(T.(x), T.(b))

function make_masked_array(x::AbstractArray)
    mask = isnan.(x)
    x = copy(x)
    x[mask] .= 0
    return MaskedArray(x, .!mask)
end
MaskedArray(T::DataType, x::AbstractArray) = make_masked_array(T.(x))
MaskedArray(x::AbstractArray) = make_masked_array(x)

data(x::MaskedArray) = x.data
mask(x::MaskedArray) = x.mask
data(x::AbstractArray) = x
mask(x::AbstractArray) = fill(true, size(x))

# A key problem is here => this works fine provided >=1 of the indices are
Base.getindex(x::MaskedArray, I...) = MaskedArray(x.data[I...], x.mask[I...])
Base.size(x::MaskedArray) = size(data(x))
Base.size(x::MaskedArray, j::Int) = size(data(x), j)
Base.length(x::MaskedArray) = length(data(x))
Base.show(io::IO, ::MIME"text/plain", x::MaskedArray) = (println("Masked array with data:");
    flush(stdout); display(as_nan(x)))
Base.display(io::IO, x::MaskedArray) = (println("Masked array with data:");
    flush(stdout);  display(as_nan(x)))

function _as_missing(x::MaskedArray{T,N}) where {T,N}
    y::AbstractArray{Union{Missing, T}} = copy(x.data)
    y[.!x.mask] .= missing
    y
end

# (Possibly) useful ops re applying the mask
as_nan(x::MaskedArray{T,N}) where {T,N} = replace(_as_missing(x), missing=>T(NaN))
rm_nan(x::MaskedArray{T,N}) where {T,N} = collect(skipmissing(_as_missing(x)))
zero_maskvals(x::MaskedArray, y::AbstractArray) = y .* mask(x)
zero_maskvals(x::MaskedArray) = zero_maskvals(x, data(x))
zero_maskvals(x::AbstractArray) = x
apply_mask(x::MaskedArray, y::AbstractArray) = MaskedArray(zero_maskvals(x, y), mask(x))

stacked(x::MaskedArray) = vcat(data(x), mask(x))
stacked(x::MaskedArray, I...) = vcat(data(x)[I...], mask(x)[I...])
hstacked(x::MaskedArray) = hcat(data(x), mask(x))
hstacked(x::MaskedArray, I...) = hcat(data(x)[I...], mask(x)[I...])
const vstacked = stacked

# # For vcat in encode method
# Base.vcat(x::MaskedArray, y::MaskedArray) = MaskedArray(vcat(data(x), data(y)),
#     vcat(mask(x), mask(y)))
Base.reshape(x::MaskedArray, sz::Int...) = MaskedArray(reshape(data(x), sz...),
    reshape(mask(x), sz...))
Base.reshape(x::MaskedArray, sz::Tuple{Vararg{Int,N}}) where N = MaskedArray(
    reshape(data(x), sz...), reshape(mask(x), sz...))

Base.repeat(x::MaskedArray, a...; kws...) = MaskedArray(repeat(data(x), a...; kws...),
    repeat(mask(x), a...; kws...))

# For comparison and use as Dict keys
Base.hash(x::MaskedArray, h::UInt) = hash((data(x), mask(x)), h)
Base.:(==)(x::MaskedArray, y::MaskedArray) = hash(x) == hash(y)

# Standard operators
Base.Array(x::MaskedArray) = as_nan(x)
Base.:+(x::MaskedArray, y::AbstractArray) = MaskedArray(data(x) + y, mask(x))
Base.:+(x::AbstractArray, y::MaskedArray) = MaskedArray(x + data(y), mask(y))
Base.:-(x::MaskedArray, y::AbstractArray) = MaskedArray(data(x) - y, mask(x))
Base.:-(x::AbstractArray, y::MaskedArray) = MaskedArray(x - data(y), mask(y))
Base.:*(x::MaskedArray, y::AbstractArray) = MaskedArray(data(x) * y, mask(x))
Base.:*(x::AbstractArray, y::MaskedArray) = MaskedArray(x * data(y), mask(y))
Base.:/(x::MaskedArray, y::AbstractArray) = MaskedArray(data(x) / y, mask(x))
Base.:/(x::AbstractArray, y::MaskedArray) = MaskedArray(x / data(y), mask(y))
Base.adjoint(x::MaskedArray) = MaskedArray(adjoint(data(x)), adjoint(mask(x)))

Base.sum(x::MaskedArray) = sum(zero_maskvals(x))
Base.sum(x::MaskedArray; dims) = sum(zero_maskvals(x), dims=dims)
StatsBase.mean(x::MaskedArray; dims) = sum(x, dims=dims) ./ sum(mask(x), dims=dims)
function StatsBase.var(x::MaskedArray; dims)
    mu = mean(x, dims=dims)
    return sum(z->z^2, (data(x) .- mu) .* mask(x), dims=dims) ./ (sum(mask(x), dims=dims) .- 1)
end
StatsBase.std(x::MaskedArray; dims) = sqrt.(var(x; dims=dims))

Base.selectdim(A::MaskedArray, d::Integer, i) = MaskedArray(Array(selectdim(data(A), d, i)),
    Array(selectdim(mask(A), d, i)))
end


text:
function isPalindrome(str) :: Bool
    len = length(str)
    head = if isodd(len)
        str[1:Int(ceil(len/2) - 1)]
    else
        str[1:Int(ceil(len/2))]
    end

    str[Int(ceil(len/2))+1:len] == reverse(head)
end

function shortestPalindrome(str) :: String
    str == "" && return ""
    len = length(str)
    reverseStr = reverse(str)
    for ind in 1:len
        isPalindrome(reverseStr[ind:len]) &&
            return join([reverseStr[1:ind-1], str])
    end
end

using Test
@test shortestPalindrome("aacecaaa") == "aaacecaaa"
@test shortestPalindrome("abcd") == "dcbabcd"
@test shortestPalindrome("") == ""
@test shortestPalindrome("aba") == "aba"
@test shortestPalindrome("abbacd") == "dcabbacd"
@test shortestPalindrome("abb") == "bbabb"
@test shortestPalindrome("aaaaa") == "aaaaa"
@test shortestPalindrome("aabba") == "abbaabba"
@test shortestPalindrome("a") == "a"
@test shortestPalindrome("babbbabbaba") == "ababbabbbabbaba"
@test shortestPalindrome("aaaabbaa") == "aabbaaaabbaa"

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
#=
  Instances for assignment.jl

=#

assignment_instances = Dict(

# Data from
# Winston "Operations Research", Assignment Problems, page 393f
# added the fifth column
# See http://www.hakank.org/minizinc/assignment.mzn
:1 => Dict(:op => :minimize,
           :m => [[14,  5, 8,  7, 15],
                       [ 2, 12, 6,  5,  3],
                       [ 7,  8, 3,  9,  7],
                       [ 2,  4, 6, 10,  1]]),

#
# Winston "Operations Research", page 398, swimming team example
# (original version]
# See http://www.hakank.org/minizinc/assignment2.mzn
#
:2 => Dict(:op => :minimize,
           :m => [[54, 54, 51, 53],
                     [51, 57, 52, 52],
                     [50, 53, 54, 56],
                     [56, 54, 55, 53]]),

#
# Winston "Operations Research", page 398, swimming team example
# See http://www.hakank.org/minizinc/assignment2_2.mzn
# expanded version
#
:3 => Dict(:op => :minimize,
            :m => [[54, 54, 51, 53,   50,60,70,80,90,100],
                     [51, 57, 52, 52,   40,50,60,70,80, 90],
                     [50, 53, 54, 56,   40,50,60,80,93, 69],
                     [56, 54, 55, 53,   60,80,40,60,50,100]]),

#
# Winston "Operations Research", page 399
#
# """
# Tom Cruise, Freddy Prinze Jr, Harrison Ford, and Matt LeBlanc
# are marooned on a desert island with Jennifer Anniston,
# Courtney Cos, Gwynneth Paltrow, and Julia Roberts.
# The 'compatibility matrix' in Table 52 indicate how much happiness
# each couple would experience if the spend all their time toghether.
# The happiness earned by a couple is proportional to the fraction
# of time the spend toghether.
# ...
# The optimal solution requires that that each person send all their
# time with one person of the opposite sex, so this result is often
# referred to as the Marriage Theorem.
# """
#
# See http://www.hakank.org/minizinc/assignment3.mzn

# males:
# 1 "Tom Cruise"
# 2 "Freddie Prinz Jr"
# 3 "Harrison Ford"
# 4 "Mark LeBlanc"
#
# females:
# 1 "Jennifer Anniston"
# 2 "Courtney Cox"
# 3 "Gwynneth Paltrow"
# 4 "Julia Roberts"
:4 => Dict(:op => :maximize,
            :m =>
                       [[7, 5, 8, 2],
                        [7, 8, 9, 4],
                        [3, 5, 7, 9],
                        [5, 5, 6, 7]]),

# From
#  "SAS OR 9.1 User's Guide Mathematical Programming"
# """
# Consider assigning five programmers to five programming jobs. Each
# programmer prefers specific programming job over others. [...]
# Suppose you ask each programmer to rank the jobs according to preference
# (using 1 for the most preferred job and 5 for the least preffered job].
# PROC ASSIGN maximizes the total preference of the group by minimizing the
# sum of the preferences.
#
#    PROGRAMMER     JOB1 JOB2 JOB3 JOB4 JOB5
#    PROGRAMMER1    4    1    3    5    2
#              2    2    1    3    4    5
#              3    3    2    4    1    5
#              4    2    3    4    5    1
#              5    4    2    3    1    5
#
# """
#
# See http://www.hakank.org/minizinc/assignment5.mzn
#
:5 => Dict(:op => :minimize,
            :m => [[4, 1, 3, 5, 2],
                     [2, 1, 3, 4, 5],
                     [3, 2, 4, 1, 5],
                     [2, 3, 4, 5, 1],
                     [4, 2, 3, 1, 5]]),

#
# From GLPK:s example assign.mod:
# """
# The assignment problem is one of the fundamental combinatorial
# optimization problems.
#
# In its most general form, the problem is as follows:
#
# There are a number of agents and a number of tasks. Any agent can be
# assigned to perform any task, incurring some cost that may vary
# depending on the agent-task assignment. It is required to perform all
# tasks by assigning exactly one agent to each task in such a way that
# the total cost of the assignment is minimized.
#
# (From Wikipedia, the free encyclopedia.]
# """
#
# """
# These data correspond to an example from [Christofides].
# """
#
# See http://www.hakank.org/minizinc/assignment6.mzn
#
:6 => Dict(:op => :minimize,
               :m => [[13, 21, 20, 12,  8, 26, 22, 11],
                    [12, 36, 25, 41, 40, 11,  4,  8],
                    [35, 32, 13, 36, 26, 21, 13, 37],
                    [34, 54,  7,  8, 12, 22, 11, 40],
                    [21,  6, 45, 18, 24, 34, 12, 48],
                    [42, 19, 39, 15, 14, 16, 28, 46],
                    [16, 34, 38,  3, 34, 40, 22, 24],
                    [26, 20,  5, 17, 45, 31, 37, 43]])

)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Tests from Flux, probably not the optimal testset organisation!

@testset "CUDA" begin
  x = randn(5, 5)
  cx = cu(x)
  @test cx isa CuArray

  @test_skip onecold(cu([1.0, 2.0, 3.0])) == 3  # passes with CuArray with Julia 1.6, but fails with JLArray

  x = onehotbatch([1, 2, 3], 1:3)
  cx = cu(x)
  @test cx isa OneHotMatrix && cx.indices isa CuArray
  @test (cx .+ 1) isa CuArray

  xs = rand(5, 5)
  ys = onehotbatch(1:5,1:5)
  @test collect(cu(xs) .+ cu(ys)) ≈ collect(xs .+ ys)
end

@testset "onehot gpu" begin
  y = onehotbatch(ones(3), 1:2) |> cu;
  @test (repr("text/plain", y); true)

  gA = rand(3, 2) |> cu;
  @test_broken gradient(A -> sum(A * y), gA)[1] isa CuArray  # fails with JLArray, bug in Zygote?
end

@testset "onecold gpu" begin
  y = onehotbatch(ones(3), 1:10) |> cu;
  l = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j']
  @test onecold(y) isa CuArray
  @test y[3,:] isa CuArray
  @test onecold(y, l) == ['a', 'a', 'a']
end

@testset "onehot forward map to broadcast" begin
  oa = OneHotArray(rand(1:10, 5, 5), 10) |> cu
  @test all(map(identity, oa) .== oa)
  @test all(map(x -> 2 * x, oa) .== 2 .* oa)
end

@testset "show gpu" begin
  x = onehotbatch([1, 2, 3], 1:3)
  cx = cu(x)
  # 3-arg show
  @test contains(repr("text/plain", cx), "1  ⋅  ⋅")
  @test contains(repr("text/plain", cx), string(typeof(cx.indices)))
  # 2-arg show, https://github.com/FluxML/Flux.jl/issues/1905
  @test repr(cx) == "Bool[1 0 0; 0 1 0; 0 0 1]"
end


text:
############################################################
# Softmax Layer
############################################################
@defstruct SoftmaxLayer Layer (
  name :: AbstractString = "softmax",
  (dim :: Int = -2, dim != 0),
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops))
)
@characterize_layer(SoftmaxLayer,
  can_do_bp => true,
)

type SoftmaxLayerState <: LayerState
  layer      :: SoftmaxLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}

  dims  :: Vector{Int}
  etc   :: Any
end

function setup_etc(backend::CPUBackend, layer::SoftmaxLayer, dims::Vector{Int}, data_type, inputs)
  nothing
end

function setup(backend::Backend, layer::SoftmaxLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type  = eltype(inputs[1])
  blobs      = Blob[make_blob(backend, data_type, size(input)) for input in inputs]
  blobs_diff = Blob[make_blob(backend, data_type, size(input)) for input in inputs]
  dims = map(inputs) do input
    total_dim = ndims(input)
    dim = layer.dim < 0 ? layer.dim + total_dim + 1 : layer.dim
    @assert 1 <= dim <= total_dim
    @assert dim != total_dim # should not operate on the mini-batch dimension
    dim
  end

  etc = setup_etc(backend, layer, dims, data_type, inputs)

  state = SoftmaxLayerState(layer, blobs, blobs_diff, dims, etc)
  return state
end
function shutdown(backend::CPUBackend, state::SoftmaxLayerState)
  map(destroy, state.blobs)
end

function forward(backend::CPUBackend, state::SoftmaxLayerState, inputs::Vector{Blob})
  for ii = 1:length(inputs)
    input  = inputs[ii].data
    output = state.blobs[ii].data
    op_dim = state.dims[ii]

    dim_pre, dim_prob, dim_post = split_dims(input, op_dim)

    for i = 0:dim_pre-1
      for j = 0:dim_post-1
        idx = Int[i + dim_pre*(k + dim_prob*j) for k=0:dim_prob-1] + 1

        maxval = -Inf
        for k in idx
          @inbounds maxval = max(maxval, input[k])
        end
        for k in idx
          @inbounds output[k] = exp(input[k]-maxval)
        end
        the_sum = 0.0
        for k in idx
          @inbounds the_sum += output[k]
        end
        for k in idx
          @inbounds output[k] /= the_sum
        end
      end
    end
  end
end

function backward(backend::CPUBackend, state::SoftmaxLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  for ii = 1:length(inputs)
    diff = diffs[ii]
    if isa(diff, CPUBlob)
      top_diff = state.blobs_diff[ii].data
      output = state.blobs[ii].data

      dim_pre, dim_prob, dim_post = split_dims(output, state.dims[ii])
      for i = 0:dim_pre-1
        for j = 0:dim_post-1
          idx = Int[i + dim_pre*(k + dim_prob*j) for k=0:dim_prob-1] + 1
          dot_prod = 0.0
          for k in idx
            @inbounds dot_prod += top_diff[k] * output[k]
          end

          for k in idx
            @inbounds diff.data[k] = (top_diff[k]-dot_prod)*output[k]
          end
        end
      end
    end
  end
end


text:
function f=GradGrad(OrdPoly)
DG.OrdPoly=OrdPoly;
[DG.w,DG.xw]=GaussLobattoQuad(OrdPoly);
[DG.DW,DG.DS,DG.DV,DG.DVT]=DerivativeMatrix(DG);
c=zeros(OrdPoly+1,1);
for i=1:OrdPoly+1
  c(i)=(i-OrdPoly/2-1)^2;
end
f=zeros(OrdPoly+1,1);
for i=1:OrdPoly+1
  for l=1:OrdPoly+1
    grad_phi_i=DG.DS(l,i);
    w=DG.w(l);
    t=0;
    for j=1:OrdPoly+1
      t=t+c(j)*DG.DS(l,j);
    end
    f(i)=f(i)+w*t*grad_phi_i;
  end
  f(i)=f(i)/DG.w(i);
end
f2=DG.DS*c;
f2=DG.DW*f2;
end


text:
function sqOut(A::Assoc)
    AtA = deepcopy(A)
    if ! isa(A.val[1], Number)
        AtA = logical(A)
    end
    AA = Adj(AtA)
    AAtAA = AA * AA';
    AtA.A = AAtAA;
    AtA.col = AtA.row
    return AtA
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
# polygon benchmarks
using Luxor, BenchmarkTools, TimerOutputs

const to = TimerOutput()

function testpoly()
    pgon = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    @timeit to "polycentroid" begin
        for i in 1:10
            pc = polycentroid(pgon)
        end
    end

    @timeit to "polyportion" begin
        for i in 1:10
            ps = polyportion(pgon, i / 10)
        end
    end

    @timeit to "polyremainder" begin
        for i in 1:10
            ps = polyremainder(pgon, i / 10)
        end
    end

    @timeit to "polyarea" begin
        for i in 1:10
            pa = polyarea(pgon)
        end
    end

    @timeit to "polyperimeter" begin
        for i in 1:10
            pa = polyperimeter(pgon)
        end
    end

    @timeit to "polydistances" begin
        for i in 1:10
            pa = polydistances(pgon)
        end
    end

    @timeit to "polysample" begin
        for i in 1:10
            pa = polysample(pgon, 200)
        end
    end

    @timeit to "polyrotate" begin
        for i in 1:10
            polyrotate!(pgon, π / 5, center = O + (20, 20))
        end
    end

    pgon1 = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    @timeit to "polyreflect" begin
        for i in 1:10
            polyreflect!(pgon1, O, O + (49.9, 51.1))
        end
    end

    pgon1 = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    @timeit to "polymove" begin
        for i in 1:10
            polymove!(pgon1, O, O + (2000.1, 1999.9))
        end
    end

    pgon1 = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    @timeit to "polyscale" begin
        for i in 1:10
            polyscale!(pgon1, 2, 0.5, center = O + (200.1, 199.9))
        end
    end

    @timeit to "polysortbyangle" begin
        for i in 1:10
            polysortbyangle(pgon)
        end
    end

    @timeit to "polysortbydistance" begin
        for i in 1:10
            polysortbydistance(pgon, pgon[1])
        end
    end

    @timeit to "polyorientation" begin
        for i in 1:10
            polyorientation(pgon)
        end
    end

    @timeit to "polysplit" begin
        for i in 1:10
            polysplit(pgon, O - (0, 100), O - (0, 100))
        end
    end

    @timeit to "polyfit" begin
        for i in 1:10
            polyfit(pgon, 200)
        end
    end

    @timeit to "polysmooth" begin
        for i in 1:10
            polysmooth(pgon, 5, :none)
        end
    end

    @timeit to "offsetpoly" begin
        for i in 1:10
            offsetpoly(pgon, 5)
        end
    end

    @timeit to "polyremovecollinearpoints" begin
        for i in 1:10
            polyremovecollinearpoints(pgon)
        end
    end

    pgon1 = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    pgon2 = star(O, 201.9, 23, 0.8, π / 5, vertices = true)
    @timeit to "polyintersect" begin
        for i in 1:10
            polyintersect(pgon1, pgon2)
        end
    end

    pgon1 = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
    pgon2 = star(O, 201.9, 23, 0.8, π / 5, vertices = true)
    @timeit to "polyintersections" begin
        for i in 1:10
            polyintersections(pgon1, pgon2)
        end
    end

    @timeit to "polytriangulate" begin
    pgons = []
    pgon = star(O, 201.9, 23, 0.8, π / 2, vertices = true)
        for i in 1:10
            pgons = polytriangulate(pgon)
        end
    end

    return pgon1, pgon2, pgons

end

# polyfit

function benchone()
    pgon = star(O, 201.9, 23, 0.8, π / 2, :stroke)
    sethue("red")
    p1 = p2 = p3 = Point[]
    for i in 1:10
        p1, p2, p3 = testpoly()
    end
    return p1, p2, p3
end

function drawtest()
    @png begin
        reset_timer!(to)
        p1, p2, p3 = benchone()
        print_timer(to)
        println()
        poly(p1, :stroke)
        scale(0.6)
        poly(p2, :stroke)
        sethue("green")
        scale(0.7)
        for i in p3
            randomhue()
            poly(i, :fill)
        end
    end 500 500 "/tmp/polybenchmark.png"
end

drawtest()

#=
────────────────────────────────────────────────────────────────────────────────────
                                             Time                   Allocations
                                     ──────────────────────   ───────────────────────
          Tot / % measured:               132ms / 99.4%            109MiB / 100%

 Section                     ncalls     time   %tot     avg     alloc   %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────
 polytriangulate                 10   81.9ms  62.6%  8.19ms   87.3MiB  80.6%  8.73MiB
 polyintersect                   10   17.5ms  13.4%  1.75ms   1.26MiB  1.17%   129KiB
 polyintersections               10   16.6ms  12.7%  1.66ms   1.20MiB  1.11%   123KiB
 polysmooth                      10   4.95ms  3.78%   495μs   2.67MiB  2.46%   273KiB
 polyfit                         10   3.14ms  2.40%   314μs   9.50MiB  8.77%   973KiB
 polysortbydistance              10   2.53ms  1.94%   253μs   1.93MiB  1.78%   198KiB
 polysample                      10   1.00ms  0.77%   100μs   1.21MiB  1.11%   124KiB
 polyorientation                 10    681μs  0.52%  68.1μs   1.26MiB  1.17%   129KiB
 polyportion                     10    578μs  0.44%  57.8μs    302KiB  0.27%  30.2KiB
 offsetpoly                      10    491μs  0.38%  49.1μs   79.7KiB  0.07%  7.97KiB
 polyremovecollinearpoints       10    352μs  0.27%  35.2μs    546KiB  0.49%  54.6KiB
 polyremainder                   10    294μs  0.23%  29.4μs    299KiB  0.27%  29.9KiB
 polysortbyangle                 10    273μs  0.21%  27.3μs    177KiB  0.16%  17.7KiB
 polysplit                       10    265μs  0.20%  26.5μs    314KiB  0.28%  31.4KiB
 polyreflect                     10   80.5μs  0.06%  8.05μs     0.00B  0.00%    0.00B
 polydistances                   10   51.3μs  0.04%  5.13μs    128KiB  0.12%  12.8KiB
 polyperimeter                   10   50.8μs  0.04%  5.08μs    128KiB  0.12%  12.8KiB
 polyarea                        10   24.0μs  0.02%  2.40μs     0.00B  0.00%    0.00B
 polyrotate                      10   10.6μs  0.01%  1.06μs     0.00B  0.00%    0.00B
 polycentroid                    10   9.86μs  0.01%   986ns     0.00B  0.00%    0.00B
 polymove                        10   5.53μs  0.00%   553ns     0.00B  0.00%    0.00B
 polyscale                       10   4.26μs  0.00%   426ns     0.00B  0.00%    0.00B
 ────────────────────────────────────────────────────────────────────────────────────
=#


text:
"""
Example: a AR(2) time series with 1000 observations. The variance quadruples after 500 observations.
"""

# Packages
using Plots
using Distributions

# Generate an AR(2) time series with a doubling in standard deviation after half the observations observations:
# y(t) = -1.0 + 0.5 * y(t-1) - 0.5 * y(t- 2) + Σ
α = -1.0
β = [0.5, -0.5]
b = length(β)
T = 1_000
y = zeros(T, 1);
Σ = ones(T)
for t in 1:T
    if t <= b
        y[t] = α + rand(Normal(0.0, 1.0))
    else
        if t <= T/2
            # Initial volatility (1.0)
            y[t] = α + y[t-1:-1:t-b]' * β + rand(Normal(0.0, 1.0))
        else
            # Double volatility (2.0)
            Σ[t] = 4.0
            y[t] = α + y[t-1:-1:t-b]' * β + rand(Normal(0.0, sqrt(Σ[t])))
        end
    end
end

# Construct priors
priors = TVVAR(y, fill(0.0, b+1, 1), Matrix(1000.0I, b+1, b+1), Matrix(1.0I, 1, 1), 0.95, 0.95)

# Estimate
est = estimate(priors)

# Plot simulated data and estimated means
pl = scatter(est.y, color = :blue, markeralpha = 0.5, label = "observed", title = "Time series: AR(2)", legend = :topleft)
pl = plot!(est.μ, color = :blue, linewidth = 2; label = "predicted")

savefig(pl,"./example/timeseries_uni.png")

# Plot the estimated vs. true coefficients
cols = [:blue :red :green :yellow]
pl = plot(est.m[:, :, 1], ylim = [-2, 2], color = cols, label = "", title = "Coefficients")
pl = plot!(fill(α, T), label = "a", color = cols[1], linestyle = :dash)
pl = plot!(fill(β[2], T), label = "b(1)", color = cols[2], linestyle = :dash)
pl = plot!(fill(β[1], T), label = "b(2)", color = cols[3], linestyle = :dash)

savefig(pl,"./example/coefficients_uni.png")

# Plot true and estimated variances
pl = plot(est.Σ[:, 1, 1], color = :blue, linewidth = 1, label = "", ylim = [0.0, 10.0], title = "Variance")
pl = plot!(Σ, color = :blue, linestyle = :dash, label = "")

savefig(pl,"./example/variance_uni.png")


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations_euler = CompressibleEulerEquations2D(gamma)

initial_condition = initial_condition_sedov_self_gravity
boundary_conditions = boundary_condition_sedov_self_gravity

surface_flux = flux_hll
volume_flux  = flux_chandrashekar
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations_euler, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver_euler = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-4, -4)
coordinates_max = ( 4,  4)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=100_000,
                periodicity=false)

semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, solver_euler,
                                          boundary_conditions=boundary_conditions)

###############################################################################
# semidiscretization of the hyperbolic diffusion equations
equations_gravity = HyperbolicDiffusionEquations2D()

# TODO: Taal, define initial/boundary conditions here for gravity?

solver_gravity = DGSEM(polydeg, flux_lax_friedrichs)

semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_condition, solver_gravity,
                                            boundary_conditions=boundary_conditions,
                                            source_terms=source_terms_harmonic)

###############################################################################
# combining both semidiscretizations for Euler + self-gravity
parameters = ParametersEulerGravity(background_density=0.0, # aka rho0
                                    gravitational_constant=6.674e-8, # aka G
                                    cfl=2.4,
                                    resid_tol=1.0e-4,
                                    n_iterations_max=100,
                                    timestep_gravity=timestep_gravity_erk52_3Sstar!)

semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters)

###############################################################################
# ODE solvers, callbacks etc.
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()

amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=1.0,
                                          alpha_min=0.0,
                                          alpha_smooth=false,
                                          variable=density_pressure)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=2,
                                      max_level =8, max_threshold=0.0003)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=1,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.0)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

analysis_interval = 100
alive_callback = AliveCallback(analysis_interval=analysis_interval)

analysis_callback = AnalysisCallback(semi_euler, interval=analysis_interval,
                                     save_analysis=true,
                                     extra_analysis_integrals=(energy_total, energy_kinetic, energy_internal))

# TODO: Taal decide, first AMR or save solution etc.
callbacks = CallbackSet(summary_callback, amr_callback, stepsize_callback,
                        save_restart, save_solution,
                        analysis_callback, alive_callback)

###############################################################################
# run the simulation
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary
println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout)


text:
using LinearAlgebra
using KrylovKit
using OMEinsum
using Random
using Zygote

"""
tensor order graph: from left to right, top to bottom.
```
a ────┬──── c    a──────┬──────c     a─────b
│     b     │    │      │      │     │     │
├─ d ─┼─ e ─┤    │      b      │     ├──c──┤
│     g     │    │      │      │     │     │
f ────┴──── h    d──────┴──────e     d─────e
```
"""

"""
    i, j = ktoij(k,Ni,Nj)
    LinearIndices -> CartesianIndices
"""
function ktoij(k,Ni,Nj)
    Cart = CartesianIndices((1:Ni,1:Nj))
    Index = Cart[k]
    Index[1],Index[2]
end

safesign(x::Number) = iszero(x) ? one(x) : sign(x)

"""
    qrpos(A)

Returns a QR decomposition, i.e. an isometric `Q` and upper triangular `R` matrix, where `R`
is guaranteed to have positive diagonal elements.
"""
qrpos(A) = qrpos!(copy(A))
function qrpos!(A)
    mattype = _mattype(A)
    F = qr!(mattype(A))
    Q = mattype(F.Q)
    R = F.R
    phases = safesign.(diag(R))
    rmul!(Q, Diagonal(phases))
    lmul!(Diagonal(conj!(phases)), R)
    return Q, R
end

"""
    lqpos(A)

Returns a LQ decomposition, i.e. a lower triangular `L` and isometric `Q` matrix, where `L`
is guaranteed to have positive diagonal elements.
"""
lqpos(A) = lqpos!(copy(A))
function lqpos!(A)
    mattype = _mattype(A)
    F = qr!(mattype(A'))
    Q = mattype(mattype(F.Q)')
    L = mattype(F.R')
    phases = safesign.(diag(L))
    lmul!(Diagonal(phases), Q)
    rmul!(L, Diagonal(conj!(phases)))
    return L, Q
end

function cellones(A)
    Ni, Nj = size(A)
    D = size(A[1,1],1)
    Cell = Array{_arraytype(A[1,1]){ComplexF64,2},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        Cell[i,j] = _mattype(A){ComplexF64}(I, D, D)
    end
    return Cell
end

function ρmap(ρ,Ai,J)
    Nj = size(Ai,1)
    for j = 1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        ρ = ein"(dc,csb),dsa -> ab"(ρ,Ai[jr],conj(Ai[jr]))
    end
    return ρ
end

function initialA(M, D)
    Ni, Nj = size(M)
    arraytype = _arraytype(M[1,1])
    A = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        d = size(M[i,j], 4)
        A[i,j] = arraytype(rand(ComplexF64, D, d, D))
    end
    return A
end
"""
    getL!(A,L; kwargs...)

````
┌ A1─A2─    ┌      L ─
ρ │  │    = ρ   =  │
┕ A1─A2─    ┕      L'─
````

ρ=L'*L, return L, where `L`is guaranteed to have positive diagonal elements.
L = cholesky!(ρ).U
If ρ is not exactly positive definite, cholesky will fail
"""
function getL!(A,L; kwargs...)
    Ni,Nj = size(A)
    for j = 1:Nj, i = 1:Ni
        λ,ρs,info = eigsolve(ρ->ρmap(ρ,A[i,:],j), L[i,j]'*L[i,j], 1, :LM; ishermitian = false, maxiter = 1, kwargs...)
        @debug "getL eigsolv" λ info sort(abs.(λ))
        info.converged == 0 && @warn "getL not converged"
        ρ = ρs[1] + ρs[1]'
        ρ ./= tr(ρ)
        F = svd!(ρ)
        Lo = lmul!(Diagonal(sqrt.(F.S)), F.Vt)
        _, L[i,j] = qrpos!(Lo)
    end
    return L
end

"""
    getAL(A,L)

Given an MPS tensor `A` and `L` ，return a left-canonical MPS tensor `AL`, a gauge transform `R` and
a scalar factor `λ` such that ``λ AR R = L A``
"""
function getAL(A,L)
    Ni,Nj = size(A)
    arraytype = _arraytype(A[1,1])
    AL = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    Le = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    λ = zeros(Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        D, d, = size(A[i,j])
        Q, R = qrpos!(reshape(L[i,j]*reshape(A[i,j], D, d*D), D*d, D))
        AL[i,j] = reshape(Q, D, d, D)
        λ[i,j] = norm(R)
        Le[i,j] = rmul!(R, 1/λ[i,j])
    end
    return AL, Le, λ
end

function getLsped(Le, A, AL; kwargs...)
    Ni,Nj = size(A)
    L = Array{_arraytype(A[1,1]){ComplexF64,2},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        λ , Ls, info = eigsolve(X -> ein"(dc,csb),dsa -> ab"(X,A[i,j],conj(AL[i,j])), Le[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "getLsped eigsolve" λ info sort(abs.(λ))
        info.converged == 0 && @warn "getLsped not converged"
        _, L[i,j] = qrpos!(Ls[1])
    end
    return L
end

"""
    leftorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)

Given an MPS tensor `A`, return a left-canonical MPS tensor `AL`, a gauge transform `L` and
a scalar factor `λ` such that ``λ AL L = L A``, where an initial guess for `L` can be
provided.
"""
function leftorth(A,L=cellones(A); tol = 1e-12, maxiter = 100, kwargs...)
    L = getL!(A,L; kwargs...)
    AL, Le, λ = getAL(A,L;kwargs...)
    numiter = 1
    while norm(L.-Le) > tol && numiter < maxiter
        L = getLsped(Le, A, AL; kwargs...)
        AL, Le, λ = getAL(A, L; kwargs...)
        numiter += 1
    end
    L = Le
    return AL, L, λ
end

"""
    rightorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)

Given an MPS tensor `A`, return a gauge transform R, a right-canonical MPS tensor `AR`, and
a scalar factor `λ` such that ``λ R AR^s = A^s R``, where an initial guess for `R` can be
provided.
"""
function rightorth(A,L=cellones(A); tol = 1e-12, maxiter = 100, kwargs...)
    Ni,Nj = size(A)
    arraytype = _arraytype(A[1,1])
    Ar = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    Lr = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        Ar[i,j] = permutedims(A[i,j],(3,2,1))
        Lr[i,j] = permutedims(L[i,j],(2,1))
    end
    AL, L, λ = leftorth(Ar,Lr; tol = tol, maxiter = maxiter, kwargs...)
    R = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    AR = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        R[i,j] = permutedims(L[i,j],(2,1))
        AR[i,j] = permutedims(AL[i,j],(3,2,1))
    end
    return R, AR, λ
end

"""
    LRtoC(L,R)

```
 ── Cᵢⱼ ──  =  ── Lᵢⱼ ── Rᵢⱼ₊₁ ──
```
"""
function LRtoC(L, R)
    Ni, Nj = size(L)
    arraytype = _arraytype(L[1,1])
    C = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    for j in 1:Nj,i in 1:Ni
        jr = j + 1 - (j + 1 > Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
    end
    return C
end

"""
    FLm = FLmap(ALi, ALip, Mi, FL, J)

ALip means ALᵢ₊₁
```
  ┌──        ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁   ──   ...          a ────┬──── c
  │          │     │        │                         │     b     │
 FLm   =   FLᵢⱼ ─ Mᵢⱼ   ── Mᵢⱼ₊₁    ──   ...          ├─ d ─┼─ e ─┤
  │          │     │        │                         │     g     │
  ┕──        ┕──  ALᵢ₊₁ⱼ ─ ALᵢ₊₁ⱼ₊₁ ──   ...          f ────┴──── h
```
"""
function FLmap(ALi, ALip, Mi, FL, J)
    Nj = size(ALi,1)
    FLm = copy(FL)
    for j=1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        FLm = ein"((adf,abc),dgeb),fgh -> ceh"(FLm,ALi[jr],Mi[jr],ALip[jr])
    end
    return FLm
end

"""
    FRm = FRmap(ARi, ARip, Mi, FR, J)

ARip means ARᵢ₊₁
```
 ──┐       ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐            a ────┬──── c
   │                │          │      │            │     b     │
──FRm  =   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ  ──FRᵢⱼ          ├─ d ─┼─ e ─┤
   │                │          │      │            │     g     │
 ──┘       ... ─ ARᵢ₊₁ⱼ₋₁ ─ ARᵢ₊₁ⱼ  ──┘            f ────┴──── h
```
"""
function FRmap(ARi, ARip, Mi, FR, J)
    Nj = size(ARi,1)
    FRm = copy(FR)
    for j=1:Nj
        jr = J-(j-1) + (J-(j-1) < 1)*Nj
        FRm = ein"((ceh,abc),dgeb),fgh -> adf"(FRm,ARi[jr],Mi[jr],ARip[jr])
    end
    return FRm
end

function FLint(AL, M)
    Ni,Nj = size(AL)
    arraytype = _arraytype(AL[1,1])
    FL = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        D = size(AL[i,j],1)
        dL = size(M[i,j],1)
        FL[i,j] = arraytype(rand(ComplexF64, D, dL, D))
    end
    return FL
end

function FRint(AR, M)
    Ni,Nj = size(AR)
    arraytype = _arraytype(AR[1,1])
    FR = Array{arraytype{ComplexF64,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        D = size(AR[i,j],1)
        dR = size(M[i,j],3)
        FR[i,j] = arraytype(rand(ComplexF64, D, dR, D))
    end
    return FR
end

"""
    λL, FL = leftenv(ALu, ALd, M, FL = FLint(ALu,M); kwargs...)

Compute the left environment tensor for MPS A and MPO M, by finding the left fixed point
of ALu - M - ALd contracted along the physical dimension.
```
 ┌──  ALuᵢⱼ  ── ALuᵢⱼ₊₁   ──   ...         ┌──
 │     │        │                          │
FLᵢⱼ ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...  = λLᵢⱼ FLᵢⱼ
 │     │        │                          │
 ┕──  ALdᵢᵣⱼ  ─ ALdᵢᵣⱼ₊₁  ──   ...         ┕──
```
"""
leftenv(ALu, ALd, M, FL = FLint(ALu,M); kwargs...) = leftenv!(ALu, ALd, M, copy(FL); kwargs...)
function leftenv!(ALu, ALd, M, FL; kwargs...)
    Ni,Nj = size(ALu)
    λL = zeros(eltype(FL[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        λLs, FL1s, info= eigsolve(X->FLmap(ALu[i,:], conj(ALd[ir,:]), M[i,:], X, j), FL[i,j], 1, :LM; maxiter=1000 , ishermitian = false, kwargs...)
        @debug "leftenv! eigsolv" λLs info sort(abs.(λLs))
        info.converged == 0 && @warn "leftenv not converged"
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FL[i,j] = FL1s[1]
                λL[i,j] = λLs[1]
            else
                FL[i,j] = FL1s[2]
                λL[i,j] = λLs[2]
            end
        else
            FL[i,j] = FL1s[1]
            λL[i,j] = λLs[1]
        end
    end
    return λL, FL
end

"""
    λR, FR = rightenv(ARu, ARd, M, FR = FRint(ARu,M); kwargs...)

Compute the right environment tensor for MPS A and MPO M, by finding the left fixed point
of AR - M - conj(AR) contracted along the physical dimension.
```
   ... ─── ARuᵢⱼ₋₁ ── ARuᵢⱼ  ──┐          ──┐
            │          │       │            │
   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ   ──FRᵢⱼ  = λRᵢⱼ──FRᵢⱼ
            │          │       │            │
   ... ─   ARdᵢᵣⱼ₋₁ ─ ARdᵢᵣⱼ ──┘          ──┘
```
"""
rightenv(ARu, ARd, M, FR = FRint(ARu,M); kwargs...) = rightenv!(ARu, ARd, M, copy(FR); kwargs...)
function rightenv!(ARu, ARd, M, FR; kwargs...)
    Ni,Nj = size(ARu)
    λR = zeros(eltype(FR[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        λRs, FR1s, info= eigsolve(X->FRmap(ARu[i,:], conj(ARd[ir,:]), M[i,:], X, j), FR[i,j], 1, :LM;maxiter=1000 , ishermitian = false, kwargs...)
        @debug "rightenv! eigsolv" λRs info sort(abs.(λRs))
        info.converged == 0 && @warn "rightenv not converged"
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                FR[i,j] = FR1s[1]
                λR[i,j] = λRs[1]
            else
                FR[i,j] = FR1s[2]
                λR[i,j] = λRs[2]
            end
        else
            FR[i,j] = FR1s[1]
            λR[i,j] = λRs[1]
        end
    end
    return λR, FR
end

"""
    ACm = ACmap(ACij, FLj, FRj, Mj, II)

```
                                ┌─────── ACᵢⱼ ─────┐
                                │        │         │             a ────┬──── c
┌─────── ACm  ─────┐      =     FLᵢⱼ ─── Mᵢⱼ ───── FRᵢⱼ          │     b     │
│        │         │            │        │         │             ├─ d ─┼─ e ─┤
                                FLᵢ₊₁ⱼ ─ Mᵢ₊₁ⱼ ──  FRᵢ₊₁ⱼ        │     g     │
                                │        │         │             f ────┴──── h
                                .        .         .
                                .        .         .
                                .        .         .
```
"""
function ACmap(ACij, FLj, FRj, Mj, II)
    Ni = size(FLj,1)
    ACm = copy(ACij)
    for i=1:Ni
        ir = II+i-1 - (II+i-1 > Ni)*Ni
        ACm = ein"((adf,abc),dgeb),ceh -> fgh"(FLj[ir],ACm,Mj[ir],FRj[ir])
    end
    return ACm
end

"""
    Cmap(Cij, FLjp, FRj, II)

```
                    ┌────Cᵢⱼ ───┐
                    │           │            a ─── b
┌──── Cm ───┐   =   FLᵢⱼ₊₁ ──── FRᵢⱼ         │     │
│           │       │           │            ├─ c ─┤
                    FLᵢ₊₁ⱼ₊₁ ── FRᵢ₊₁ⱼ       │     │
                    │           │            d ─── e
                    .           .
                    .           .
                    .           .
```
"""
function Cmap(Cij, FLjp, FRj, II)
    Ni = size(FLjp,1)
    Cm = copy(Cij)
    for i=1:Ni
        ir = II+i-1 - (II+i-1 > Ni)*Ni
        Cm = ein"(acd,ab),bce -> de"(FLjp[ir],Cm,FRj[ir])
    end
    return Cm
end

"""
    ACenv(AC, FL, M, FR;kwargs...)

Compute the up environment tensor for MPS `FL`,`FR` and MPO `M`, by finding the up fixed point
        of `FL - M - FR` contracted along the physical dimension.
```
┌─────── ACᵢⱼ ─────┐
│        │         │
FLᵢⱼ ─── Mᵢⱼ ───── FRᵢⱼ
│        │         │
FLᵢ₊₁ⱼ ─ Mᵢ₊₁ⱼ ──  FRᵢ₊₁ⱼ  =  λACᵢⱼ ┌──── ACᵢⱼ ───┐
│        │         │                │      │      │
.        .         .
.        .         .
.        .         .
```
"""
ACenv(AC, FL, M, FR; kwargs...) = ACenv!(copy(AC), FL, M, FR; kwargs...)
function ACenv!(AC, FL, M, FR; kwargs...)
    Ni,Nj = size(AC)
    λAC = zeros(eltype(AC[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        λACs, ACs, info = eigsolve(X->ACmap(X, FL[:,j], FR[:,j], M[:,j], i), AC[i,j], 1, :LM; maxiter=1000 ,ishermitian = false, kwargs...)
        @debug "ACenv! eigsolve" λACs info sort(abs.(λACs))
        info.converged == 0 && @warn "ACenv Not converged"
        if length(λACs) > 1 && norm(abs(λACs[1]) - abs(λACs[2])) < 1e-12
            @show λACs
            if real(λACs[1]) > 0
                AC[i,j] = ACs[1]
                λAC[i,j] = λACs[1]
            else
                AC[i,j] = ACs[2]
                λAC[i,j] = λACs[2]
            end
        else
            AC[i,j] = ACs[1]
            λAC[i,j] = λACs[1]
        end
    end
    return λAC, AC
end

"""
    Cenv(C, FL, FR;kwargs...)

Compute the up environment tensor for MPS `FL` and `FR`, by finding the up fixed point
    of `FL - FR` contracted along the physical dimension.
```
┌────Cᵢⱼ ───┐
│           │
FLᵢⱼ₊₁ ──── FRᵢⱼ
│           │
FLᵢ₊₁ⱼ₊₁ ── FRᵢ₊₁ⱼ   =  λCᵢⱼ ┌──Cᵢⱼ ─┐
│           │                │       │
.           .
.           .
.           .
```
"""
Cenv(C, FL, FR; kwargs...) = Cenv!(copy(C), FL, FR; kwargs...)
function Cenv!(C, FL, FR; kwargs...)
    Ni,Nj = size(C)
    λC = zeros(eltype(C[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        jr = j + 1 - (j==Nj) * Nj
        λCs, Cs, info = eigsolve(X->Cmap(X, FL[:,jr], FR[:,j], i), C[i,j], 1, :LM; maxiter=1000 ,ishermitian = false, kwargs...)
        @debug "Cenv! eigsolve" λCs info sort(abs.(λCs))
        info.converged == 0 && @warn "Cenv Not converged"
        if length(λCs) > 1 && norm(abs(λCs[1]) - abs(λCs[2])) < 1e-12
            @show λCs
            if real(λCs[1]) > 0
                C[i,j] = Cs[1]
                λC[i,j] = λCs[1]
            else
                C[i,j] = Cs[2]
                λC[i,j] = λCs[2]
            end
        else
            C[i,j] = Cs[1]
            λC[i,j] = λCs[1]
        end
    end
    return λC, C
end

function ACCtoAL(ACij,Cij)
    D,d, = size(ACij)
    QAC, RAC = qrpos(reshape(ACij,(D*d, D)))
    QC, RC = qrpos(Cij)
    errL = norm(RAC-RC)
    # @show errL
    reshape(QAC*QC', (D, d, D)), errL
end

function ACCtoAR(ACij,Cijr)
    D,d, = size(ACij)
    LAC, QAC = lqpos(reshape(ACij,(D, d*D)))
    LC, QC = lqpos(Cijr)
    errR = norm(LAC-LC)
    # @show errR
    reshape(QC'*QAC, (D, d, D)), errR
end

"""
    itoir(i,Ni,Nj)

````
i -> (i,j) -> (i,jr) -> ir
````
"""
function itoir(i,Ni,Nj)
    Liner = LinearIndices((1:Ni,1:Nj))
    Cart = CartesianIndices((1:Ni,1:Nj))
    Index = Cart[i]
    i,j = Index[1],Index[2]
    jr = j - 1 + (j==1)*Nj
    Liner[i,jr]
end

function ALCtoAC(AL,C)
    Ni,Nj = size(AL)
    ACij = [ein"asc,cb -> asb"(AL[i],C[i]) for i=1:Ni*Nj]
    reshape(ACij,Ni,Nj)
end

"""
    AL, AR = ACCtoALAR(AC, C)

QR factorization to get `AL` and `AR` from `AC` and `C`

````
──ALᵢⱼ──Cᵢⱼ──  =  ──ACᵢⱼ──  = ──Cᵢ₋₁ⱼ ──ARᵢⱼ──
  │                  │                  │
````
"""
function ACCtoALAR(AC, C)
    Ni,Nj = size(AC)
    ALijerrL = [ACCtoAL(AC[i],C[i]) for i=1:Ni*Nj]
    AL = reshape([ALijerrL[i][1] for i=1:Ni*Nj],Ni,Nj)
    errL = Zygote.@ignore sum([ALijerrL[i][2] for i=1:Ni*Nj])
    ARijerrR = [ACCtoAR(AC[i],C[itoir(i,Ni,Nj)]) for i=1:Ni*Nj]
    AR = reshape([ARijerrR[i][1] for i=1:Ni*Nj],Ni,Nj)
    errR = Zygote.@ignore sum([ARijerrR[i][2] for i=1:Ni*Nj])
    return AL, AR, errL, errR
end

"""
Compute the error through all environment `AL,C,AR,FL,M,FR`

````
        ┌── AC──┐
        │   │   │           ┌── AC──┐
MAC1 =  FL─ M ──FR  =  λAC  │   │   │
        │   │   │

        ┌── AC──┐
        │   │   │           ┌──C──┐
MAC2 =  FL─ M ──FR  =  λAC  │     │
        │   │   │
         ── AR──┘

── MAC1 ──    ≈    ── MAC2 ── AR ──
    │                         │
````
"""
function error(AL,C,AR,FL,M,FR)
    Ni,Nj = size(AL)
    AC = ALCtoAC(AL, C)
    err = 0
    for j = 1:Nj, i = 1:Ni
        MAC = ACmap(AC[i,j], FL[:,j], FR[:,j], M[:,j], i)
        MAC -= ein"(apc,dpc),dsb -> asb"(MAC,conj(AR[i,j]),AR[i,j])
        err += norm(MAC)
    end
    return err
end

"""
    λL, FL = obs_FL(ALu, ALd, M, FL = FLint(AL,M); kwargs...)

Compute the observable left environment tensor for MPS A and MPO M, by finding the left fixed point
of AL - M - conj(AL) contracted along the physical dimension.
```
 ┌──  ALuᵢⱼ  ── ALuᵢⱼ₊₁   ──   ...         ┌──
 │     │        │                          │
FLᵢⱼ ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...  = λLᵢⱼ FLᵢⱼ
 │     │        │                          │
 ┕──  ALdᵢᵣⱼ  ─ ALdᵢᵣⱼ₊₁  ──   ...         ┕──
```
"""
obs_FL(ALu, ALd, M, FL = FLint(ALu,M); kwargs...) = obs_FL!(ALu, ALd, M, copy(FL); kwargs...)
function obs_FL!(ALu, ALd, M, FL; kwargs...)
    Ni,Nj = size(ALu)
    λL = zeros(eltype(FL[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = Ni + 1 - i
        λLs, FL1s, info= eigsolve(X->FLmap(ALu[i,:], ALd[ir,:], M[i,:], X, j), FL[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "obs_FL eigsolve" λLs info sort(abs.(λLs))
        info.converged == 0 && @warn "obs_FL Not converged"
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FL[i,j] = FL1s[1]
                λL[i,j] = λLs[1]
            else
                FL[i,j] = FL1s[2]
                λL[i,j] = λLs[2]
            end
        else
            FL[i,j] = FL1s[1]
            λL[i,j] = λLs[1]
        end
    end
    return λL, FL
end

"""
    λR, FR = obs_FR(ARu, ARd, M, FR = FRint(ARu,M); kwargs...)

Compute the observable right environment tensor for MPS A and MPO M, by finding the left fixed point
of AR - M - conj(AR) contracted along the physical dimension.
```
   ... ─── ARuᵢⱼ₋₁ ── ARuᵢⱼ  ──┐          ──┐
            │          │       │            │
   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ   ──FRᵢⱼ  = λRᵢⱼ──FRᵢⱼ
            │          │       │            │
   ... ─   ARdᵢᵣⱼ₋₁ ─ ARdᵢᵣⱼ ──┘          ──┘
```
"""
obs_FR(ARu, ARd, M, FR = FRint(ARu,M); kwargs...) = obs_FR!(ARu, ARd, M, copy(FR); kwargs...)
function obs_FR!(ARu, ARd, M, FR; kwargs...)
    Ni,Nj = size(ARu)
    λR = zeros(eltype(FR[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = Ni + 1 - i
        λRs, FR1s, info= eigsolve(X->FRmap(ARu[i,:], ARd[ir,:], M[i,:], X, j), FR[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "obs_FR! eigsolve" λRs info sort(abs.(λRs))
        info.converged == 0 && @warn "obs_FR! Not converged"
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                FR[i,j] = FR1s[1]
                λR[i,j] = λRs[1]
            else
                FR[i,j] = FR1s[2]
                λR[i,j] = λRs[2]
            end
        else
            FR[i,j] = FR1s[1]
            λR[i,j] = λRs[1]
        end
    end
    return λR, FR
end

function BgFLint(AL, M)
    Ni,Nj = size(AL)
    arraytype = _arraytype(AL[1,1])
    BgFL = Array{arraytype{ComplexF64,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        D1 = size(AL[i,j],1)
        D2 = size(AL[irr,j],1)
        dL1 = size(M[i,j],1)
        dL2 = size(M[ir,j],1)
        BgFL[i,j] = arraytype(rand(ComplexF64, D1, dL1, dL2, D2))
    end
    return BgFL
end

"""
    BgFLm = BgFLmap(ALi, ALip, Mi, Mip, BgFLij, J)

ALip means ALᵢ₊₂
Mip means Mᵢ₊₁
```
  ┌──        ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁    ──   ...       a ────┬──── c
  │          │     │        │                       │     b     │
  │          │  ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...       ├─ d ─┼─ e ─┤
BgFLm   =  BgFLᵢⱼ  │        │                       │     f     │
  │          │  ─ Mᵢ₊₁ⱼ ── Mᵢ₊₁ⱼ₊₁   ──   ...       ├─ g ─┼─ h ─┤
  │          │     │        │                       │     j     │
  ┕──        ┕──  ALᵢ₊₂ⱼ ─ ALᵢ₊₂ⱼ₊₁  ──   ...       i ────┴──── k
```
"""
function BgFLmap(ALi, ALip, Mi, Mip, BgFLij, J)
    Nj = size(ALi,1)
    BgFLm = copy(BgFLij)
    for j=1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        BgFLm = ein"(((adgi,abc),dfeb),gjhf),ijk -> cehk"(BgFLm,ALi[jr],Mi[jr],Mip[jr],ALip[jr])
    end
    return BgFLm
end

"""
    λL, BgFL = bigleftenv(ALu, ALd, M, BgFL = BgFLint(AL,M); kwargs...)

Compute the up and down left environment tensor for MPS A and MPO M, by finding the left fixed point
of AL - M - M - conj(AL) contracted along the physical dimension.
```
   ┌──  ALuᵢⱼ ── ALuᵢⱼ₊₁   ──   ...           ┌──
   │     │        │                           │
   │  ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...           │
 BgFLᵢⱼ  │        │                   = λLᵢⱼ BgFLᵢⱼ
   │  ─ Mᵢ₊₁ⱼ ── Mᵢ₊₁ⱼ₊₁   ──   ...           │
   │     │        │                           │
   ┕── ALdᵢ₊₂ⱼ ─ ALdᵢ₊₂ⱼ₊₁ ──   ...           ┕──
```
"""
bigleftenv(ALu, ALd, M, BgFL = BgFLint(ALu,M); kwargs...) = bigleftenv!(ALu, ALd, M, copy(BgFL); kwargs...)
function bigleftenv!(ALu, ALd, M, BgFL; kwargs...)
    Ni,Nj = size(ALu)
    λL = zeros(eltype(BgFL[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        λLs, BgFL1s, info= eigsolve(X->BgFLmap(ALu[i,:], ALd[irr,:], M[i,:], M[ir,:], X, j), BgFL[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "bigleftenv! eigsolve" λLs info sort(abs.(λLs))
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                BgFL[i,j] = BgFL1s[1]
                λL[i,j] = λLs[1]
            else
                BgFL[i,j] = BgFL1s[2]
                λL[i,j] = λLs[2]
            end
        else
            BgFL[i,j] = BgFL1s[1]
            λL[i,j] = λLs[1]
        end
    end
    return λL, BgFL
end

function BgFRint(AR, M)
    Ni,Nj = size(AR)
    arraytype = _arraytype(AR[1,1])
    BgFR = Array{arraytype{ComplexF64,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        D1 = size(AR[i,j],3)
        D2 = size(AR[irr,j],3)
        dR1 = size(M[i,j],3)
        dR2 = size(M[ir,j],3)
        BgFR[i,j] = arraytype(rand(ComplexF64, D1, dR1, dR2, D2))
    end
    return BgFR
end

"""
    FRm = FRmap(ARi, ARip, Mi, FR, J)

ARip means ARᵢ₊₁
```
 ──┐          ...  ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐        a ────┬──── c
   │                    │          │      │        │     b     │
 ──│          ... ──── Mᵢⱼ₋₁  ──  Mᵢⱼ   ──│        ├─ d ─┼─ e ─┤
  BgFRm   =             │          │     BgFRm     │     f     │
 ──│          ... ──── Mᵢ₊₁ⱼ₋₁ ── Mᵢ₊₁ⱼ ──│        ├─ g ─┼─ h ─┤
   │                    │          │      │        │     j     │
 ──┘          ...  ─ ARᵢ₊₂ⱼ₋₁ ─── ARᵢ₊₂ⱼ──┘        i ────┴──── k
```
"""
function BgFRmap(ARi, ARip, Mi, Mip, BgFR, J)
    Nj = size(ARi,1)
    BgFRm = copy(BgFR)
    for j=1:Nj
        jr = J-(j-1) + (J-(j-1) < 1)*Nj
        BgFRm = ein"(((cehk,abc),dfeb),gjhf),ijk -> adgi"(BgFRm,ARi[jr],Mi[jr],Mip[jr],ARip[jr])
    end
    return BgFRm
end

"""
    λR, BgFR = bigrightenv(ARu, ARd, M, BgFR = BgFRint(ARu,M); kwargs...)

Compute the up and down right environment tensor for MPS A and MPO M, by finding the left fixed point
of AR - M - M - conj(AR) contracted along the physical dimension.
```
     ──┐          ...  ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐
       │                    │          │      │
     ──│          ... ──── Mᵢⱼ₋₁  ──  Mᵢⱼ   ──│
λRᵢⱼ BgFRᵢⱼ   =             │          │     BgFRᵢⱼ
     ──│          ... ──── Mᵢ₊₁ⱼ₋₁ ── Mᵢ₊₁ⱼ ──│
       │                    │          │      │
     ──┘          ...  ─ ARᵢ₊₂ⱼ₋₁ ─── ARᵢ₊₂ⱼ──┘
```
"""
bigrightenv(ARu, ARd, M, BgFR = BgFRint(ARu,M); kwargs...) = bigrightenv!(ARu, ARd, M, copy(BgFR); kwargs...)
function bigrightenv!(ARu, ARd, M, BgFR; kwargs...)
    Ni,Nj = size(ARu)
    λR = zeros(eltype(BgFR[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        λRs, BgFR1s, info= eigsolve(X->BgFRmap(ARu[i,:], ARd[irr,:], M[i,:], M[ir,:], X, j), BgFR[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "bigrightenv eigsolve" λRs info sort(abs.(λRs))
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                BgFR[i,j] = BgFR1s[1]
                λR[i,j] = λRs[1]
            else
                BgFR[i,j] = BgFR1s[2]
                λR[i,j] = λRs[2]
            end
        else
            BgFR[i,j] = BgFR1s[1]
            λR[i,j] = λRs[1]
        end
    end
    return λR, BgFR
end

function norm_FLint(AL)
    Ni,Nj = size(AL)
    arraytype = _arraytype(AL[1,1])
    norm_FL = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        D = size(AL[i,j],1)
        norm_FL[i,j] = arraytype(rand(ComplexF64, D, D))
    end
    return norm_FL
end

"""
    FL_normm = norm_FLmap(ALui, ALdi, FL_norm, J)

```
   ┌──        ┌──  ALuᵢⱼ ── ALuᵢⱼ₊₁ ──  ...     a───┬───c
  FLm   =    FLm     │        │                 │   b   │
   ┕──        ┕──  ALdᵢⱼ ── ALdᵢⱼ₊₁ ──  ...     d───┴───e
```
"""
function norm_FLmap(ALui, ALdi, FL_norm, J)
    Nj = size(ALui,1)
    FL_normm = copy(FL_norm)
    for j=1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        FL_normm = ein"(ad,abc),dbe -> ce"(FL_normm,ALui[jr],ALdi[jr])
    end
    return FL_normm
end

norm_FL(ALu, ALd, FL_norm = norm_FLint(ALu); kwargs...) = norm_FL!(ALu, ALd, FL_norm; kwargs...)
function norm_FL!(ALu, ALd, FL_norm; kwargs...)
    Ni,Nj = size(ALu)
    λL = zeros(eltype(FL_norm[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        λLs, FL_norms, info= eigsolve(X->norm_FLmap(ALu[i,:], ALd[i,:], X, j), FL_norm[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "norm_FL eigsolve" λLs info sort(abs.(λLs))
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FL_norm[i,j] = FL_norms[1]
                λL[i,j] = λLs[1]
            else
                FL_norm[i,j] = FL_norms[2]
                λL[i,j] = λLs[2]
            end
        else
            FL_norm[i,j] = FL_norms[1]
            λL[i,j] = λLs[1]
        end
    end
    return λL, FL_norm
end

function norm_FRint(AR)
    Ni,Nj = size(AR)
    arraytype = _arraytype(AR[1,1])
    norm_FR = Array{arraytype{ComplexF64,2},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        D = size(AR[i,j],1)
        norm_FR[i,j] = arraytype(rand(ComplexF64, D, D))
    end
    return norm_FR
end

"""
    FR_normm = norm_FRmap(ARui, ARdi, FR_norm, J)

```
──┐       ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐      a───┬───c
 FRm  =            │          │     FRm     │   b   │
──┘       ... ─── ARᵢⱼ₋₁ ─── ARᵢⱼ  ──┘      d───┴───e
```
"""
function norm_FRmap(ARui, ARdi, FR_norm, J)
    Nj = size(ARui,1)
    FR_normm = copy(FR_norm)
    for j=1:Nj
        jr = J-(j-1) + (J-(j-1) < 1)*Nj
        FR_normm = ein"(ce,abc),dbe -> ad"(FR_normm,ARui[jr],ARdi[jr])
    end
    return FR_normm
end

norm_FR(ARu, ARd, FR_norm = norm_FRint(ARu); kwargs...) = norm_FR!(ARu, ARd, FR_norm; kwargs...)
function norm_FR!(ARu, ARd, FR_norm; kwargs...)
    Ni,Nj = size(ARu)
    λL = zeros(eltype(FR_norm[1,1]),Ni,Nj)
    for j = 1:Nj, i = 1:Ni
        λLs, FR_norms, info= eigsolve(X->norm_FRmap(ARu[i,:], ARd[i,:], X, j), FR_norm[i,j], 1, :LM; ishermitian = false, kwargs...)
        @debug "norm_FR! eigsolve" λLs info sort(abs.(λLs))
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FR_norm[i,j] = FR_norms[1]
                λL[i,j] = λLs[1]
            else
                FR_norm[i,j] = FR_norms[2]
                λL[i,j] = λLs[2]
            end
        else
            FR_norm[i,j] = FR_norms[1]
            λL[i,j] = λLs[1]
        end
    end
    return λL, FR_norm
end


text:
using AtomEnergyLevels
using PyPlot
pygui(true)

function main(z, rgrid; config = atom[z], xc = xc_func(), n = 500)
    x, en, HOMO = log.(rgrid), similar(rgrid), similar(rgrid)
    for (i, xmax) in enumerate(x)
        results = lda(z, range(-35, xmax, length = n), conf = config, xc = xc)
        en[i] = results.energy.total
        HOMO[i] = maximum(map(λ -> λ[1], values(results.orbitals)))
    end
    return en, HOMO
end

#=
Jorge Garza, Rubicelia Vargas, and Alberto Vela.
Numerical self-consistent-field method to solve the Kohn-Sham equations
in confined many-electron atoms // Phys. Rev. E. 58 (1998) 3949-3954
DOI:10.1103/PhysRevE.58.3949
=#

r = range(0.90, 6, length = 50)
E, HOMO = main(2, r, config = c"[He]", xc = xc_func(xc_type_lda, false, LDA_X, ρ -> 0))

begin
    title("Total energy of the helium atom for several confinements")
    xlabel("Rc (a.u.)")
    ylabel("Total energy (a.u.)")
    grid("on")

    plot(r, E)

    # TABLE III.
    plot(1:6, [1.354, -2.384, -2.682, -2.718, -2.723, -2.724], linestyle="None", marker="o",
         label="Jorge Garza et al // Phys. Rev. E. 58 (1998) 3949")
    legend(loc = "upper right", fancybox = "true")
end

#=
Jan C. A. Boeyens. Ionization Radii of Compressed Atoms.
J. CHEM. SOC. FARADAY TRANS., 1994, 90(22), 3377-3381

NB. Hartree-Fock-Slater method is simply Xα functional, where α = 1

Data taken from Fig. 4 One-electron compression curve for the valence level of oxygen.
=#

r = [2.2059, 2.2752, 2.3446, 2.4002, 2.4696, 2.5526,
     2.6494, 2.7598, 2.8839, 3.0218, 3.1733, 3.3522,
     3.5722, 3.8195, 4.0805, 4.3688, 4.6980, 5.0272,
     5.3564, 5.6992, 6.0283, 6.3711, 6.7139, 7.0429,
     7.3582, 7.7010, 8.0437, 8.3728, 8.7019, 9.0446,
     9.3874]

ϵₕ = [0.95604,  0.81394,  0.66592,  0.51786,  0.37576,
      0.23370,  0.09166, -0.03848, -0.17452, -0.31052,
     -0.44649, -0.57054, -0.68264, -0.78281, -0.87110,
     -0.94154, -0.97038, -0.99330, -1.01622, -1.02725,
     -1.03239, -1.03749, -1.04259, -1.03587, -1.04104,
     -1.04021, -1.03939, -1.04452, -1.04373, -1.04291,
     -1.04208] # Ry

E, HOMO = main(8, r, config = c"[He] 2s2 2p4", xc = xc_func(xc_type_lda, false, ρ -> LDA_X(ρ, α = 1), ρ -> 0))

begin
    ax1 = PyPlot.axes()
    ax1.set_xlim([ 2.0, 10])
    ax1.set_ylim([-1.5,  1])
    plot(r, ϵₕ)
    plot(r, 2HOMO)
    grid("on")
end


text:
# Callbacks functions

function fbs_trace(t, θ, loss, reg)
    @printf("iter %5d : loss = %8.4f, reg = %8.4f, total = %8.4f, ",
            t, loss, reg, loss+reg)
    @printf("non zero = %d\n", sum(θ .!= 0.))
end

function fista_trace(t, θ, loss, reg, α)
    @printf("iter %5d : loss = %8.4f, reg = %8.4f, total = %8.4f, ",
            t, loss, reg, loss+reg)
    @printf("α = %4.2f, non zero = %d\n",
            α, sum(θ .!= 0.))
end

function fasta_trace(t, θ, loss, reg, λ, α, bt)
    @printf("iter %5d : loss = %8.4f, reg = %8.4f, total = %8.4f, ",
            t, loss, reg, loss+reg)
    @printf("λ = %8.5f, α = %4.2f, backtracks = %2d, non zero = %d\n",
            λ, α, bt, sum(θ .!= 0.))
end


text:
## Leer datos
using Random

## Leer datos
include("../../utils/preprocesado_datos.jl")
include("../../utils/validation.jl")
include("../../algoritmo-enfriamiento/enfriamiento-simulado.jl")

file_path= "src/Instancias_APC/"
struct DataFile
    route::String
    class_atributte::String
end

files = [
    DataFile(file_path*"ionosphere.arff", "class"),
    DataFile(file_path*"parkinsons.arff", "Class"),
    DataFile(file_path*"spectf-heart.arff", "Class"),
]
# ¡Tiene que acabar en /
csv_file_path  = "src/resultados/Enfriamiento-Simulado/"
prefijo = "ES-"
process_name = [
    DataFile(csv_file_path*prefijo*"ionosphere.result.csv", "Datos Iosfera"),
    DataFile(csv_file_path*prefijo*"parkinsons.result.csv", "Datos parkinson "),
    DataFile(csv_file_path*prefijo*"spectf-heart.result.csv", "Datos ataques corazón")
]
Random.seed!(0)
println("Procedemos a calcular los datos con $(Threads.nthreads()) hebras disponibles")

for i in 1:length(files)
    # Seleccionamos parkinson
    file = files[i]
    data , labels = DataLabelArff(file.route, file.class_atributte)
    ## Entrenar para obtener el peso
    println("==============================================")
    println("$(process_name[i].class_atributte)")

    ## Evaluar resultado con cross validation
    VerboseCrossValidation(data, labels, 5, EnfriamientoSimulado_LearnerOneNN, process_name[i].route)
end


text:
using Jecco.AdS5_3_1

grid = SpecCartGrid3D(
    x_min            = -20.0,
    x_max            =  20.0,
    x_nodes          =  48,
    y_min            = -20.0,
    y_max            =  20.0,
    y_nodes          =  48,
    u_outer_min      =  0.1,
    u_outer_max      =  1.005,
    u_outer_domains  =  3,
    u_outer_nodes    =  32,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

potential = Phi8Potential(
    oophiM2 = -1.0,
    oophiQ  = 0.1,
)

id = BlackBranePert(
    energy_dens = 0.817,
    phi0        = 1.0,
    phi2        = 0.29819,
    oophiM2     = -1.0,
    a4_ampx     = -5.e-2,
    a4_ampy     = -5.e-2,
    a4_kx       = 1,
    a4_ky       = 1,
    AH_pos      = 1.0,
    xmax        = grid.x_max,
    xmin        = grid.x_min,
    ymin        = grid.y_min,
    ymax        = grid.y_max,
)

evoleq = AffineNull(
    phi0           = 1.0,
    potential      = potential,
    gaugecondition = ConstantAH(u_AH = 1.0),
)

diag = DiagAH(
    find_AH_every_t    = 1.0,
)

io = InOut(
    out_boundary_every_t        = 0.2,
    out_bulk_every_t            = 5.0,
    out_gauge_every_t           = 0.2,
    out_bulkconstrained_every_t = 5.0,
    checkpoint_every_walltime_hours = 1,
)

integration = Integration(
    tmax            = 500.0,
    ODE_method      = AdS5_3_1.VCABM3(),
    adaptive        = true,
    filter_poststep = true,
)

run_model(grid, id, evoleq, diag, integration, io)


text:
using Plots
using RecipesBase

function get_width_bounds(mid_point_x, interval_width)
    half_width = 0.5 * interval_width
    return mid_point_x .+ [-half_width, half_width]
end

function get_interval_vertices(width_bounds, length_bounds)
    lower_bound_width, upper_bound_width = width_bounds
    lower_bound_length, upper_bound_length = length_bounds
    vertices_x = [
        lower_bound_width,
        lower_bound_width,
        upper_bound_width,
        upper_bound_width
    ]
    vertices_y = [
        lower_bound_length,
        upper_bound_length,
        upper_bound_length,
        lower_bound_length
    ]
    return vertices_x, vertices_y
end

@shorthands confidenceinterval

@recipe function f(
        ::Type{Val{:confidenceinterval}},
        args...;
        interval_width=0.75,
        interval_color=:skyblue,
        line_color=:black
)
    interval_data = args[2]
    mid_point_x = plotattributes[:series_plotindex]
    width_bounds = get_width_bounds(mid_point_x, interval_width)

    legend := false
    @series begin
        # Confidence Interval Bar
        seriestype := :shape
        seriesalpha --> 0.5
        fillcolor := interval_color
        length_bounds = [interval_data[1], interval_data[3]]
        vertices = get_interval_vertices(width_bounds, length_bounds)
        x := vertices[1]
        y := vertices[2]
        ()
    end

    @series begin
        seriestype := :shape
        seriescolor := line_color
        line_value = interval_data[2]
        x := width_bounds
        y := [line_value, line_value]
        ()
    end
end


text:
"""

Increases size of arrays in dc structure if needed

"""
function increase_size!(dc::CutoffDistances, maxdim::Int)
  dc.maxdim[1] = maxdim
  resize!(dc.d,dc.maxdim[1])
  resize!(dc.iat,dc.maxdim[1])
  resize!(dc.jat,dc.maxdim[1])
  resize!(dc.imol,dc.maxdim[1])
  resize!(dc.jmol,dc.maxdim[1])
  @. dc.d[dc.nd[1]+1:dc.maxdim[1]] = 0.
  @. dc.iat[dc.nd[1]+1:dc.maxdim[1]] = 0
  @. dc.jat[dc.nd[1]+1:dc.maxdim[1]] = 0
  @. dc.imol[dc.nd[1]+1:dc.maxdim[1]] = 0
  @. dc.jmol[dc.nd[1]+1:dc.maxdim[1]] = 0
  nothing
end

# Increase size by 50% of current size

function increase_size!(dc::CutoffDistances)
  maxdim = round(Int,1.5*dc.maxdim[1])
  increase_size!(dc, maxdim)
  nothing
end


text:
module ModStokesFastDirect

using ModifiedStokesSolver

# Dummy module. Fast direct solver not included in code.

end


text:
using ChaosTools
using Test

"""
    saturation_point(x, y; threshold = 0.01, dxi::Int = 1, tol = 0.2)
Decompose the curve `y(x)` into linear regions using `linear_regions(x, y; dxi, tol)`
and then attempt to find a saturation point where the the first slope
of the linear regions become `< threshold`.

Return the `x` value of the saturation point.
"""
function saturation_point(Ds, E1s; threshold = 0.01, kwargs...)
    lrs, slops = ChaosTools.linear_regions(Ds, E1s; kwargs...)
    i = findfirst(x -> x < threshold, slops)
    return i == 0 ? Ds[end] : Ds[lrs[i]]
end

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

@testset "Estimate Delay" begin

    ds = Systems.henon()
    data = trajectory(ds,1000;dt=1)
    x = data[:,1]
    @test estimate_delay(x,"first_zero") <= 2
    @test estimate_delay(x,"first_min")  <= 2
    @test estimate_delay(x,"exp_decay")  <= 2
    # @test 3 <= estimate_delay(x,"mutual_inf"; k=1) <= 10
    # @test 3 <= estimate_delay(x,"mutual_inf"; k=10) <= 10

    ds = Systems.roessler(ones(3))
    dt = 0.01
    data = trajectory(ds,200,dt=dt)
    x = data[:,1]
    @test 1.3 <= estimate_delay(x,"first_zero")*dt <= 1.7
    @test 2.6 <= estimate_delay(x,"first_min")*dt  <= 3.4

    dt = 0.1
    data = trajectory(ds,2000,dt=dt)
    x = data[:,1]
    @test 1.3 <= estimate_delay(x,"first_zero")*dt <= 1.7
    @test 2.6 <= estimate_delay(x,"first_min")*dt  <= 3.4
    @test 1.3 <= estimate_delay(x,"mutual_inf")*dt <= 1.8
    @test 1.3 <= estimate_delay(x,"mutual_inf"; maxtau=150, k = 2)*dt <= 1.7

    ds = Systems.lorenz()
    dt = 0.05
    data = trajectory(ds,1000;dt=dt)
    x = data[500:end,1]
    @test 0.1 <= estimate_delay(x,"mutual_inf")*dt <= 0.4
    @test 0.1 <= estimate_delay(x,"exp_decay")*dt  <= 0.4
    #
    dt = 0.1
    data = trajectory(ds,2000;dt=dt)
    x = data[:,1]
    @test 0.1 <= estimate_delay(x,"exp_decay")*dt  <= 0.4
end

@testset "Estimate Dimension" begin
    s = sin.(0:0.1:1000)
    τ = 15
    Ds = 1:5
    E1s = estimate_dimension(s, τ, Ds)
    E2s = stochastic_indicator(s, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.01) == 1

    # @test minimum(E2s) < 0.1 # THIS TEST FAILS

    ds = Systems.roessler();τ=15; dt=0.1
    data = trajectory(ds,1000;dt=dt)
    s = data[:,1]
    Ds = 1:5
    E1s = estimate_dimension(s, τ, Ds)
    E2s = stochastic_indicator(s, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]
    @test minimum(E2s) < 0.3

    ds = Systems.lorenz();τ=5; dt=0.01
    data = trajectory(ds,500;dt=dt)
    s = data[:,1]
    Ds = 1:5
    E1s = estimate_dimension(s, τ, Ds)
    E2s = stochastic_indicator(s, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]

    # @test minimum(E2s) < 0.1 # THIS TEST FAILS

    #Test against random signal
    E2s = stochastic_indicator(rand(100000), 1, 1:5)
    @test minimum(E2s) > 0.9

end


text:
include("../src/ttv.jl")
include("/Users/ericagol/Computer/Julia/regress.jl")

n = 8
t0 = 7257.93115525
#h  = 0.12
#h  = 0.075
h  = 0.05
tmax = 600.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+2
end
println("ntt: ",ntt)
tt1 = zeros(n,maximum(ntt))
tt2 = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
@time ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
@time ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
# Write out every 10th step to a file:
@time ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar,"test_out.txt",10)
# Now call with half the timestep:
count2 = zeros(Int64,n)
ttv_elements!(n,t0,h/10.,tmax,elements,tt2,count2,0.0,0,0,rstar)
println("Timing error: ",maximum(abs.(tt1-tt2))*24.*3600.," sec")

# Check type stabilty:

dtdq0 = zeros(n,maximum(ntt),7,n)
dtdelements0 = @code_warntype ttv_elements!(n,t0,h,tmax,elements,tt1,count1,dtdq0,rstar)

using PyPlot

# Make a plot of some TTVs:

fig,axes = subplots(4,2)

for i=2:8
  ax = axes[i-1]
  fn = zeros(Float64,2,count1[i])
  sig = ones(count1[i])
  tti1 = tt1[i,1:count1[i]]
  tti2 = tt2[i,1:count2[i]]
  for j=1:count1[i]
    fn[1,j] = 1.0
    fn[2,j] = round(Int64,(tti1[j]-elements[i,3])/elements[i,2])
  end
  coeff,cov = regress(fn,tti1,sig)
  tt_ref1 = coeff[1]+coeff[2]*fn[2,:]
  ttv1 = (tti1-tt_ref1)*24.*60.
#  coeff,cov = regress(fn,tti2,sig)
  tt_ref2 = coeff[1]+coeff[2]*fn[2,:]
  ttv2 = (tti2-tt_ref2)*24.*60.
  ax[:plot](tti1,ttv1)
  ax[:plot](tti2,ttv2)
#  ax[:plot](tti2,((ttv1-ttv2)-mean(ttv1-ttv2)))
  ax[:plot](tti2,ttv1-ttv2)
  println(i," ",coeff," ",elements[i,2:3]," ",coeff[1]-elements[i,3]," ",coeff[2]-elements[i,2])
#  println(i," ",maximum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec ", minimum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec" )
  println(i," ",maximum(ttv1-ttv2)*60.," sec ", minimum(ttv1-ttv2)*60.," sec")
end