text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

#

#

# In a conjunction the left hand side vessel is called *parent* vessel and
# the right hand side vessel is called *daughter* vessel.
#
# The conjunction is solved by imposing the conservation of mass and
# __*total*__
# pressure at the interface node $j$. Two additional relations are obtained
# by extrapolating the outgoing characteristics from the two vessels. From
# parent vessel outlet we have
# $$
#   W_1 = u_1 + 4c_1,
# $$
# and for the daughter vessel inlet we have
# $$
#   W_2 = u_2 - 4c_2.
# $$
# The mass conservation reads $A_1u_1 - A_2u_2 = 0$, and the total pressure
# conservation requires $P_{t1} = P_{t2}$. By defining the unknown vector $U$
# as
# $$
#   U = \{U_i\} = \left\{ \begin{array}{c}
#                           u_1 \\
#                           u_2 \\
#                           A_1^{1/4} \\
#                           A_2^{1/4}
#                        \end{array} \right\} , \quad i =1, ..., 4,
# $$
# the four relations read
# $$
#   F = \left\{f_i \right\} =
#   \begin{cases}
#     U_1 + 4k_1U_3 - W_1^* = 0, \\
#     U_2 - 4k_2U_4 - W_2^* = 0, \\
#     U_1U_3^4 - U_2U_4^4 = 0, \\
#     \beta_1 \left(\tfrac{U_3^2}{A_{01}^{1/2}} -1 \right) +
#     \frac{1}{2}\rho U_1^2 -
#     \beta_2 \left(\tfrac{U_4^2}{A_{02}^{1/2}} -1 \right) -
#     \frac{1}{2}\rho U_2^2 = 0,
#   \end{cases}
# $$
# where $c_i = k_i A_i^{1/4}$, $k_i = \sqrt{3/2 \gamma_i}$. $F(U)=0$ is solved
# iteratively with Newton's method [^1]
# $$
#   \begin{cases}
#     J \cdot \delta U = - F(U), \\
#     U^{new} = U + \delta U,
#   \end{cases}
# $$
# where $J$ is the Jacobian
# $$
#   J = \left[ \begin{array}{cccc}
#         1 & 0 & 4k_1 & 0 \\
#         0 & 1 & 0 & -4k_2 \\
#         U_3^4 & -U_4^4 & 4U_1 U_3^3 & -4 U_2 U_4^3 \\
#         \rho U_1 & -\rho U_2 & 2 \beta_1 U_3/A_{01}^{1/2} &
#           -2\beta_2 U_4/A_{02}^{1/2}
#       \end{array} \right].
# $$
# The entire process is handled by
# [`solveConjunction`](conjunctions.html#solveConjunction) function.

# *function* __`solveConjunction`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------------- --------------------------------------------------------
# `b`                 `::Blood` model matrix row.
#
# `v1`                `::Vessel` data structure for the parent vessel.
#
# `v2`                `::Vessel` data structure for the daughter vessel.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Firstly, the unknown vector $U$ and parameters vector $k$ are initialised
# with data from parent and daughter vessels.
# ----------------------------------------------------------------------------
#
function solveConjunction(b :: Blood, v1 :: Vessel, v2 :: Vessel)

  U = [v1.u[end],
       v2.u[ 1 ],
       sqrt(sqrt(v1.A[end])),
       sqrt(sqrt(v2.A[ 1 ]))]

  k1 = sqrt(0.5*3*v1.gamma[end])
  k2 = sqrt(0.5*3*v2.gamma[ 1 ])
  k  = [k1, k2]
  # Riemann invariants are initialised by
  # [`calculateWstarConj`](conjunctions.html#calculateWstarConj), the
  # Jacobian is calculated by
  # [`calculateJacobianConj`](conjunctions.html#calculateJacobianConj), and
  # the $F$ vector is calculated by
  # [`calculateFofUconj`](conjunctions.html#calculateFofUconj), and
  W = calculateWstarConj(U, k)
  J = calculateJacobianConj(b, v1, v2, U, k)
  F = calculateFofUconj(b, v1, v2, U, k, W)
  # Newton-Raphson method iterates until a tolerance on the error is met. The
  # tolerance is defined on $U$ and $F$.
  nr_toll_U = 1.e-5
  nr_toll_F = 1.e-5
  # At the beginning of each iteration, $\delta U$ is calculated by means of
  # LU decomposition by julia's command `\`. $U$ is updated with an iterative
  # step of 0.01.
  while true
    dU = J\(-F)
    # U_new = U + 0.01*dU
    U_new = U + dU
    # In the case the solution diverges, $F$ will contain `nan` elements. In
    # this case `openBF` exits the simulation and returns an error.
    if any(isnan(dot(F,F)))
      println(F)
      break
    end
    # When all the elements in $\delta U$ and $F$ meet the set tolerances, the
    # solver stops the iterative loop.
    u_ok = 0
    f_ok = 0
    for i in 1:length(dU)
      if abs(dU[i]) <= nr_toll_U || abs(F[i]) <= nr_toll_F
        u_ok += 1
        f_ok += 1
      end
    end

    if u_ok == length(dU) || f_ok == length(dU)
      U = U_new
      break
    # If tolerances are not met, new $U$, $W$, and $F$ vectors are computed
    # and a new iteration starts.
    else
      U = U_new
      W = calculateWstarConj(U, k)
      F = calculateFofUconj(b, v1, v2, U, k, W)
    end
  end

  # Once the solution is found, quantities at both sides of the interface are
  # updated.
  v1.u[end] = U[1]
  v2.u[ 1 ] = U[2]

  v1.A[end] = U[3]*U[3]*U[3]*U[3]
  v1.Q[end] = v1.u[end]*v1.A[end]

  v2.A[ 1 ] = U[4]^4
  v2.Q[ 1 ] = v2.u[1]*v2.A[1]

  v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
  v2.P[ 1 ] = pressure(v2.A[ 1 ], v2.A0[ 1 ], v2.beta[ 1 ], v2.Pext)

  v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
  v2.c[ 1 ] = waveSpeed(v2.A[ 1 ], v2.gamma[ 1 ])
end

# *function* __`calculateWstarConj`__ $\rightarrow$
# `W::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Riemann invariants are computed starting from unknown and parameters
# vectors as
# $$
#   W_1 = u_1 + 4c_1, \quad W_2 = u_2 - 4c_2.
# $$
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `W`         `::Array` containing outgoing characteristics at both sides
#             of the conjunction.
# ----------------------------------------------------------------------------
#
function calculateWstarConj(U :: Array, k :: Array)

  W1 = U[1] + 4*k[1] * U[3]
  W2 = U[2] - 4*k[2] * U[4]

  return [W1, W2]
end

# *function* __`calculateFofUconj`__ $\rightarrow$ `F::Array`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `b`            `::Blood` data structure.
#
# `v1`           `::Vessel` parent vessel data structure.
#
# `v2`           `::Vessel` daughter vessel data structure.
#
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
#
# `W`            `::Array` outgoing characteristics array.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# $F$ array is computed by imposing the conservation of mass and total
# pressure. $F$ reads
# $$
#   F = \left\{f_i \right\} =
#   \begin{cases}
#     U_1 + 4k_1U_3 - W_1^* = 0, \\
#     U_2 - 4k_2U_4 - W_2^* = 0, \\
#     U_1U_3^4 - U_2U_4^4 = 0, \\
#     \beta_1 \left(\tfrac{U_3^2}{A_{01}^{1/2}} -1 \right) +
#     \frac{1}{2}\rho U_1^2 -
#     \beta_2 \left(\tfrac{U_4^2}{A_{02}^{1/2}} -1 \right) -
#     \frac{1}{2}\rho U_2^2 = 0,
#   \end{cases}
# $$
# ----------------------------------------------------------------------------
# Static pressure conservation can be imposed by replacing `f4` with
#
#     f4 = v1.beta*( U[3] ^2 /sqrt(v1.A0) - 1) -
#        ( v2.beta*((U[4])^2 /sqrt(v2.A0) - 1) )
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `F`         `::Array` Newton's method relations.
# ----------------------------------------------------------------------------
#
function calculateFofUconj(b :: Blood, v1 :: Vessel, v2 :: Vessel,
                           U :: Array, k :: Array, W :: Array)

  f1 = U[1] + 4*k[1]*U[3] - W[1]

  f2 = U[2] - 4*k[2]*U[4] - W[2]

  f3 = U[1]*(U[3]*U[3]*U[3]*U[3]) - U[2]*(U[4]*U[4]*U[4]*U[4])

  f4 = 0.5*b.rho*U[1]*U[1] + v1.beta[end]*(U[3]*U[3]/sqrt(v1.A0[end]) - 1) -
     ( 0.5*b.rho*U[2]*U[2] + v2.beta[ 1 ]*(U[4]*U[4]/sqrt(v2.A0[ 1 ]) - 1) )

  return [f1, f2, f3, f4]
end

# *function* __`calculateJacobianConj`__ $\rightarrow$
# `W::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `b`            `::Blood` data structure.
#
# `v1`           `::Vessel` parent vessel data structure.
#
# `v2`           `::Vessel` daughter vessel data structure.
#
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The Jacobian is computed as
# $$
#   J = \left[ \begin{array}{cccc}
#         1 & 0 & 4k_1 & 0 \\
#         0 & 1 & 0 & -4k_2 \\
#         U_3^4 & -U_4^4 & 4U_1 U_3^4 & -4 U_2 U_4^4 \\
#         \rho U_1 & -\rho U_2 & 2 \beta_1 U_3/A_{01}^{1/2} &
#           -2\beta_2 U_4/A_{02}^{1/2}
#       \end{array} \right].
# $$
# ----------------------------------------------------------------------------
#
# The conservation of static pressure is imposed by setting to zero elements
# $J_{4,1}$ and $J_{4,2}$.
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `J`         `::Array{Float, 2}` Jacobian matrix.
# ----------------------------------------------------------------------------
#
function calculateJacobianConj(b :: Blood, v1 :: Vessel, v2 :: Vessel,
                               U :: Array, k :: Array)

  J = eye(4)

  J[1,3] =  4*k[1]
  J[2,4] = -4*k[2]

  J[3,1] =  U[3]*U[3]*U[3]*U[3]
  J[3,2] = -U[4]*U[4]*U[4]*U[4]
  J[3,3] =  4*U[1]*(U[3]*U[3]*U[3]^3)
  J[3,4] = -4*U[2]*(U[4]*U[4]*U[4]^3)

  J[4,3] =  2*v1.beta[end]*U[3]/sqrt(v1.A0[end])
  J[4,4] = -2*v2.beta[ 1 ]*U[4]/sqrt(v2.A0[ 1 ])

  J[4,1] =  b.rho*U[1]
  J[4,2] = -b.rho*U[2]

  return J
end

# ### References
#
# [^1]: Press, William H. Numerical recipes 3rd edition:
# The art of scientific computing. Cambridge university press, 2007.


text:
module TabularMakie

export tplot, lplot

using AbstractPlotting
using AbstractPlotting: MakieLayout
using AbstractPlotting.MakieLayout: Optional, LegendEntry, EntryGroup
# import some internal functions of AbstractPlotting
# TODO: test behavior of these functions

using AbstractPlotting: Automatic
include("categorical.jl")
using .CategoricalConversions: categorical_labels, categorical_range, categorical_positions, categorical_trait, Categorical, HasRefPool, Continuous

using DataFrames
using Chain: @chain
using Underscores: @_
using StructArrays
using UnPack: @unpack
#using Statistics
using NamedTupleTools: delete

using InlineTest: @testset, @test

include("mini-language.jl")
include("attribute-dicts.jl")
include("ticks.jl")
include("plot.jl")
include("layout.jl")
include("legend-entries.jl")
include("legend-layout.jl")
include("compose.jl")

end

ext: jl

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
arg_dict = Dict{Any, Any}(
"run_type" => 1,
"jac_type" => 3,
"newton_globalize_euler" => true,
"order" => 1,
"real_time" => false,
"variable_type" => :conservative,
"smb_name" => "src/mesh_files/2D_airfoilnp4.smb",
"dmg_name" => "src/mesh_files/2D_airfoil.dmg",
"IC_name" => "ICFreeStream",
"numBC" => 2,
"BC1" => [8],
"BC1_name" => "FreeStreamBC",
"BC2" => [5],
"BC2_name" => "noPenetrationBC",
"aoa" => 0.,
"Ma" => 0.15,
"use_DG" => true,
"Flux_name" => "RoeFlux",
"t_max" => 500000.0,
"delta_t" => 0.0001,
"res_tol" => 1e-10,
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"step_tol" => 1e-10,
"itermax" => 5000000000, # For RK4
"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => true,
"output_freq" => 100000,
"solve" => true,
"do_postproc" => true,
"exact_soln_func" => "ICFreeStream"
)


text:
# simple wrapper to give indices a custom wrapping behaviour
struct SparseArray{T,N} <: AbstractArray{T,N}
    data::Dict{CartesianIndex{N},T}
    dims::NTuple{N,Int64}
    function SparseArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
        return new{T,N}(Dict{CartesianIndex{N},T}(), dims)
    end
    function SparseArray(a::SparseArray{T,N}) where {T,N}
        new{T,N}(copy(a.data), a.dims)
    end
end
SparseArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} =
    SparseArray{T,N}(undef, dims)
SparseArray{T}(::UndefInitializer, dims...) where {T} = SparseArray{T}(undef, dims)

nonzero_pairs(a::SparseArray) = pairs(a.data)
nonzero_keys(a::SparseArray) = keys(a.data)
nonzero_values(a::SparseArray) = values(a.data)
nonzero_length(a::SparseArray) = length(a.data)

_zero!(x::SparseArray) = empty!(x.data)
_sizehint!(x::SparseArray, n) = sizehint!(x.data, n)

# elementary getindex and setindex!
@inline function Base.getindex(a::SparseArray{T,N}, I::CartesianIndex{N}) where {T,N}
    @boundscheck checkbounds(a, I)
    return get(a.data, I, zero(T))
end
Base.@propagate_inbounds Base.getindex(a::SparseArray{T,N}, I::Vararg{Int,N}) where {T,N} =
                                        getindex(a, CartesianIndex(I))

@inline function Base.setindex!(a::SparseArray{T,N}, v, I::CartesianIndex{N}) where {T,N}
    @boundscheck checkbounds(a, I)
    if v != zero(v)
        a.data[I] = v
    else
        delete!(a.data, I) # does not do anything if there was no key corresponding to I
    end
    return v
end
Base.@propagate_inbounds Base.setindex!(a::SparseArray{T,N},
                                        v, I::Vararg{Int,N}) where {T,N} =
                                            setindex!(a, v, CartesianIndex(I))

@inline function increaseindex!(a::SparseArray{T,N}, v, I::CartesianIndex{N}) where {T,N}
    @boundscheck checkbounds(a, I)
    iszero(v) && return
    h = a.data
    index = Base.ht_keyindex2!(h, I)
    @inbounds begin
        if index > 0
            currentv = h.vals[index]
            newv = currentv + convert(T, v)
            if iszero(newv)
                Base._delete!(h, index)
            else
                h.age += 1
                h.keys[index] = I
                h.vals[index] = newv
            end
        else
            newv = convert(T, v)
            Base._setindex!(h, newv, I, -index)
        end
    end
    return newv
end

# following code is used to index with ranges etc
_newindex(i::Int, range::Int) = i == range ? () : nothing
function _newindex(i::Int, range::AbstractVector{Int})
    k = findfirst(==(i), range)
    k === nothing ? nothing : (k,)
end
_newindices(I::Tuple{}, indices::Tuple{}) = ()
function _newindices(I::Tuple, indices::Tuple)
    i = _newindex(I[1], indices[1])
    Itail = _newindices(Base.tail(I), Base.tail(indices))
    (i === nothing || Itail === nothing) && return nothing
    return (i..., Itail...)
end

_findfirstvalue(v, r) = findfirst(==(v), r)
# slicing should produce SparseArray
function Base._unsafe_getindex(::IndexCartesian, a::SparseArray{T,N},
                                I::Vararg{<:Union{Int,AbstractVector{Int}},N}) where {T,N}
    @boundscheck checkbounds(a, I...)
    indices = Base.to_indices(a, I)
    b = SparseArray{T}(undef, length.(Base.index_shape(indices...)))
    for (k, v) in a.data
        newI = _newindices(k.I, indices)
        if newI !== nothing
            b[newI...] = v
        end
    end
    return b
end

Base.Array(a::SparseArray{T,N}) where {T,N} = Array{T,N}(a)
function Base.Array{T,N}(a::SparseArray) where {T,N}
    d = fill(zero(T), size(a))
    for (I, v) in a.data
        d[I] = v
    end
    d
end

SparseArray(a::AbstractArray{T,N}) where {T,N} = SparseArray{T,N}(a)
SparseArray{T}(a::AbstractArray{<:Any,N}) where {T,N} = SparseArray{T,N}(a)
function SparseArray{T,N}(a::AbstractArray{<:Any,N}) where {T,N}
    d = SparseArray{T,N}(undef, size(a))
    for I in CartesianIndices(a)
        iszero(a[I]) && continue
        d[I] = a[I]
    end
    return d
end
Base.convert(::Type{S}, a::S) where {S<:SparseArray} = a
Base.convert(S::Type{<:SparseArray}, a::AbstractArray) = S(a)

function SparseArray(A::Adjoint{T,<:SparseArray{T,2}}) where T
    B = SparseArray{T}(undef, size(A))
    for (I, v) in parent(A).data
        B[I[2], I[1]] = conj(v)
    end
    return B
end
function SparseArray(A::Transpose{T,<:SparseArray{T,2}}) where T
    B = SparseArray{T}(undef, size(A))
    for (I, v) in parent(A).data
        B[I[2], I[1]] = v
    end
    return B
end

Base.copy(a::SparseArray) = SparseArray(a)

Base.size(a::SparseArray) = a.dims

function Base.copy!(dst::SparseArray, src::SparseArray)
    axes(dst) == axes(src) || throw(ArgumentError(
        "arrays must have the same axes for copy! (consider using `copyto!`)"))

    _zero!(dst)
    for (I, v) in nonzero_pairs(src)
        dst[I] = v
    end
    return dst
end

Base.similar(a::SparseArray, ::Type{S}, dims::Dims{N}) where {S,N} =
    SparseArray{S}(undef, dims)


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Profile
using HDF5
using Random
using PyCall
using Arpack

@testset "tpsci" begin
    atoms = []

    r = 1
    a = 1
    push!(atoms,Atom(1,"H", [0, 0*a, 0*r]))
    push!(atoms,Atom(2,"H", [0, 0*a, 1*r]))
    push!(atoms,Atom(3,"H", [0, 1*a, 2*r]))
    push!(atoms,Atom(4,"H", [0, 1*a, 3*r]))
    push!(atoms,Atom(5,"H", [0, 2*a, 4*r]))
    push!(atoms,Atom(6,"H", [0, 2*a, 5*r]))
    push!(atoms,Atom(7,"H", [0, 3*a, 6*r]))
    push!(atoms,Atom(8,"H", [0, 3*a, 7*r]))
    push!(atoms,Atom(9,"H", [0, 4*a, 8*r]))
    push!(atoms,Atom(10,"H",[0, 4*a, 9*r]))
    push!(atoms,Atom(11,"H",[0, 5*a, 10*r]))
    push!(atoms,Atom(12,"H",[0, 5*a, 11*r]))

    clusters    = [(1:2),(3:4),(5:6),(7:8),(9:12)]
    init_fspace = [(1,1),(1,1),(1,1),(1,1),(2,2)]
    clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(2,2),(2,2),(2,2)]
    na = 6
    nb = 6

    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    nroots = 1

    # get integrals
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));

    @printf(" Do FCI\n")
    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    fci = pyimport("pyscf.fci")
    cisolver = pyscf.fci.direct_spin1.FCI()
    cisolver.max_cycle = 200
    cisolver.conv_tol = 1e-8
    nelec = na + nb
    norb = size(ints.h1,1)
    #e_fci, v_fci = cisolver.kernel(ints.h1, ints.h2, norb, nelec, ecore=0, nroots =nroots)

    #e_fci = [-18.33022092,
    #         -18.05457644]
    e_fci  = [-18.33022092,
              -18.05457645,
              -18.02913047,
              -17.99661027
             ]

    for i in 1:length(e_fci)
        @printf(" %4i %12.8f %12.8f\n", i, e_fci[i], e_fci[i]+ints.h0)
    end

    # localize orbitals
    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    #FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Rotate Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    #
    # define clusters
    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))

    e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1, rdm1,
                                       max_iter_oo=40, verbose=0, gconv=1e-6, method="bfgs")
    #FermiCG.pyscf_write_molden(mol,Cl*U,filename="cmf.molden")
    ints = FermiCG.orbital_rotation(ints,U)

    e_ref = e_cmf - ints.h0

    max_roots = 100
    # build Hamiltonian, cluster_basis and cluster ops
    #display(Da)
    #cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=2, max_roots=max_roots)
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots,
                                                       init_fspace=init_fspace, rdm1a=Da, rdm1b=Db)
    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)
    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    FermiCG.add_cmf_operators!(cluster_ops, cluster_bases, ints, Da, Db);

    ref_fock = FermiCG.FockConfig(init_fspace)

    if true

        ci_vector = FermiCG.ClusteredState(clusters, ref_fock, R=nroots)

        @time e0, v0 = FermiCG.tpsci_ci(ci_vector, cluster_ops, clustered_ham, incremental=true,
                                  thresh_cipsi=1e-3, thresh_foi=1e-9, thresh_asci=1e-4, conv_thresh=1e-5, matvec=1);

        ref = [-18.32973618]

        @test isapprox(abs.(ref), abs.(e0), atol=1e-8)
    end

    nroots = 4

    ci_vector = FermiCG.ClusteredState(clusters, ref_fock, R=nroots)

    #1 excitons
    ci_vector[ref_fock][ClusterConfig([2,1,1])] = [0,1,0,0]
    ci_vector[ref_fock][ClusterConfig([1,2,1])] = [0,0,1,0]
    ci_vector[ref_fock][ClusterConfig([1,1,2])] = [0,0,0,1]

    #e0, v0 = FermiCG.tpsci_ci(ci_vector, cluster_ops, clustered_ham, incremental=false,
    #                          thresh_cipsi=1e-2, thresh_foi=1e-4, thresh_asci=1e-2, conv_thresh=1e-4);
    e0, v0 = FermiCG.tpsci_ci(ci_vector, cluster_ops, clustered_ham, incremental=true,
                              thresh_cipsi=1e-2, thresh_foi=1e-4, thresh_asci=1e-2, conv_thresh=1e-4);

    ci_vector4 = FermiCG.extract_roots(v0,[4])
    # todo now add root following for tpsci

    e2, v1 = FermiCG.compute_pt2(v0, cluster_ops, clustered_ham, thresh_foi=1e-8, matvec=3)

    ref = [-18.32932467
           -18.05349474
           -18.02775313
           -17.99514933
          ]
    @test isapprox(abs.(ref), abs.(e0+e2), atol=1e-7)

    rotations = FermiCG.hosvd(v0, cluster_ops)
    for ci in clusters
        FermiCG.rotate!(cluster_ops[ci.idx], rotations[ci.idx])
        FermiCG.rotate!(cluster_bases[ci.idx], rotations[ci.idx])
        FermiCG.check_basis_orthogonality(cluster_bases[ci.idx])
    end

    #cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);
    #FermiCG.add_cmf_operators!(cluster_ops, cluster_bases, ints, Da, Db);

    e0a, v0a = FermiCG.tpsci_ci(ci_vector, cluster_ops, clustered_ham, incremental=true,
                                thresh_cipsi=1e-2, thresh_foi=1e-4, thresh_asci=1e-2);

    H = FermiCG.build_full_H(v0a, cluster_ops, clustered_ham)
    sig1 = H*FermiCG.get_vectors(v0a)
    sig2 = FermiCG.tps_ci_matvec(v0a, cluster_ops, clustered_ham)

    @test isapprox(norm(sig1-sig2), 0.0, atol=1e-12)

    guess = deepcopy(v0a)
    FermiCG.rand!(guess)
    e0b, v0b = FermiCG.tps_ci_direct(guess, cluster_ops, clustered_ham, conv_thresh=1e-9);
    e0c, v0c = FermiCG.tps_ci_davidson(guess, cluster_ops, clustered_ham, conv_thresh=1e-9);

    @test isapprox(abs.(e0a), abs.(e0b), atol=1e-9)
    @test isapprox(abs.(e0a), abs.(e0c), atol=1e-9)

    println(" Now test pt2 correction")

    ref = [-18.32916288
           -18.05357935
           -18.02800015
           -17.99499973]

    e2a, v1a = FermiCG.compute_pt2(v0a, cluster_ops, clustered_ham, thresh_foi=1e-8, matvec=3)
    @test isapprox(abs.(ref), abs.(e0a+e2a), atol=1e-7)

    e2b = FermiCG.compute_batched_pt2(v0a, cluster_ops, clustered_ham, thresh_foi=1e-8)
    @test isapprox(abs.(ref), abs.(e0a+e2b), atol=1e-7)

    ci_vector = FermiCG.ClusteredState(clusters, ref_fock, R=nroots)
    ci_vector[ref_fock][ClusterConfig([2,1,1])] = [0,1,0,0]
    ci_vector[ref_fock][ClusterConfig([1,2,1])] = [0,0,1,0]
    ci_vector[ref_fock][ClusterConfig([1,1,2])] = [0,0,0,1]

    sig1 = FermiCG.open_matvec_serial2(ci_vector, cluster_ops, clustered_ham, nbody=4, thresh=1e-8)
    sig2 = FermiCG.open_matvec_thread(ci_vector, cluster_ops, clustered_ham, nbody=4, thresh=1e-8)
    sig3 = FermiCG.open_matvec_thread2(ci_vector, cluster_ops, clustered_ham, nbody=4, thresh=1e-8)

    @test isapprox(norm(sig1), norm(sig2), atol=1e-16)
    @test isapprox(norm(sig1), norm(sig3), atol=1e-16)

    @test isapprox(norm(FermiCG.get_vectors(sig1) - FermiCG.get_vectors(sig2)), 0.0, atol=1e-16)
    @test isapprox(norm(FermiCG.get_vectors(sig1) - FermiCG.get_vectors(sig3)), 0.0, atol=1e-16)

end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
module FwiFlow
    using Conda
    using PyCall
    using Reexport
    @reexport using ADCME
    using LinearAlgebra
    using PyPlot
    using Random
    using JSON
    using DataStructures
    using Dierckx
    using Parameters
    using MAT

    export DATADIR
    DATADIR = "$(@__DIR__)/../docs/data"

    function ADCME.:Session(args...;kwargs...)
        config = tf.ConfigProto(
            device_count = Dict("GPU"=> 0) # do not use any GPU devices for all ops except FWI
        )
        sess = tf.Session(config=config)
    end
    include("Core.jl")
    include("Utils.jl")
    include("FWI.jl")

end


text:
module MultiplicationPartialBase

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export multiplication

function multiplication(m)

    #tot_cover_of_i += x_axis_covered*y_axis_covered

    @expression(m, static_new_variable_1[i=1:length(boxes), j=1:length(boxes); i!=j], 0.5*(m[:x_axis_covered][i, j] + m[:y_axis_covered][i, j]))
    @expression(m, static_new_variable_2[i=1:length(boxes), j=1:length(boxes); i!=j], 0.5*(m[:x_axis_covered][i, j] - m[:y_axis_covered][i, j]))

    # Piecewise linearisation start, need two, one for new each variable
    @variable(m, piecewise1[i=1:length(boxes), j=1:length(boxes); i!=j])
    @variable(m, static_lambda_box_cover1[i=1:length(boxes), j=1:length(boxes), k=1:10; i!=j] >= 0)
    @expression(m, static_lambda_box_cover1_row[i=1:length(boxes), j=1:length(boxes); i!=j],
        [
        m[:static_lambda_box_cover1][i, j, 1],
        m[:static_lambda_box_cover1][i, j, 2],
        m[:static_lambda_box_cover1][i, j, 3],
        m[:static_lambda_box_cover1][i, j, 4],
        m[:static_lambda_box_cover1][i, j, 5],
        m[:static_lambda_box_cover1][i, j, 6],
        m[:static_lambda_box_cover1][i, j, 7],
        m[:static_lambda_box_cover1][i, j, 8],
        m[:static_lambda_box_cover1][i, j, 9],
        m[:static_lambda_box_cover1][i, j, 10]
        ]
    )
    @constraint(m, static_lambda_box_cover_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j],
        sum(m[:static_lambda_box_cover1_row][i, j]) == 1
    )
    @constraint(m, static_lambda_box_cover1_sos[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:static_lambda_box_cover1_row][i, j] in SOS2()
    )

    @constraint(m, static_product_lambda_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j], (
        -m[:static_new_variable_1][i, j] +
        lin_plus[1]*m[:static_lambda_box_cover1][i, j, 2] +
        lin_plus[2]*m[:static_lambda_box_cover1][i, j, 3] +
        lin_plus[3]*m[:static_lambda_box_cover1][i, j, 4] +
        lin_plus[4]*m[:static_lambda_box_cover1][i, j, 5] +
        lin_plus[5]*m[:static_lambda_box_cover1][i, j, 6] +
        lin_plus[6]*m[:static_lambda_box_cover1][i, j, 7] +
        lin_plus[7]*m[:static_lambda_box_cover1][i, j, 8] +
        lin_plus[8]*m[:static_lambda_box_cover1][i, j, 9] +
        lin_plus[9]*m[:static_lambda_box_cover1][i, j, 10] == 0)
    )
    @constraint(m, static_piecewise_lambda_constraint1[i=1:length(boxes), j=1:length(boxes); i!=j], (
        -m[:piecewise1][i, j] +
        lin_plus_2[1]*m[:static_lambda_box_cover1][i, j, 2] +
        lin_plus_2[2]*m[:static_lambda_box_cover1][i, j, 3] +
        lin_plus_2[3]*m[:static_lambda_box_cover1][i, j, 4] +
        lin_plus_2[4]*m[:static_lambda_box_cover1][i, j, 5] +
        lin_plus_2[5]*m[:static_lambda_box_cover1][i, j, 6] +
        lin_plus_2[6]*m[:static_lambda_box_cover1][i, j, 7] +
        lin_plus_2[7]*m[:static_lambda_box_cover1][i, j, 8] +
        lin_plus_2[8]*m[:static_lambda_box_cover1][i, j, 9] +
        lin_plus_2[9]*m[:static_lambda_box_cover1][i, j, 10] == 0)
    )

    @variable(m, piecewise2[i=1:length(boxes), j=1:length(boxes); i!=j], )
    @variable(m, static_lambda_box_cover2[i=1:length(boxes), j=1:length(boxes), k=1:10; i!=j] >= 0)
    @expression(m, static_lambda_box_cover2_row[i=1:length(boxes), j=1:length(boxes); i!=j],
        [
        m[:static_lambda_box_cover2][i, j, 1],
        m[:static_lambda_box_cover2][i, j, 2],
        m[:static_lambda_box_cover2][i, j, 3],
        m[:static_lambda_box_cover2][i, j, 4],
        m[:static_lambda_box_cover2][i, j, 5],
        m[:static_lambda_box_cover2][i, j, 6],
        m[:static_lambda_box_cover2][i, j, 7],
        m[:static_lambda_box_cover2][i, j, 8],
        m[:static_lambda_box_cover2][i, j, 9],
        m[:static_lambda_box_cover2][i, j, 10]
        ]
    )
    @constraint(m, static_lambda_box_cover_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        sum(m[:static_lambda_box_cover2_row][i, j]) == 1
    )
    @constraint(m, static_lambda_box_cover2_sos[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:static_lambda_box_cover2_row][i, j] in SOS2()
    )

    @constraint(m, static_product_lambda_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        (
        -m[:static_new_variable_2][i, j] +
        lin_min[1]*m[:static_lambda_box_cover2][i, j, 2] +
        lin_min[2]*m[:static_lambda_box_cover2][i, j, 3] +
        lin_min[3]*m[:static_lambda_box_cover2][i, j, 4] +
        lin_min[4]*m[:static_lambda_box_cover2][i, j, 5] +
        lin_min[5]*m[:static_lambda_box_cover2][i, j, 6] +
        lin_min[6]*m[:static_lambda_box_cover2][i, j, 7] +
        lin_min[7]*m[:static_lambda_box_cover2][i, j, 8] +
        lin_min[8]*m[:static_lambda_box_cover2][i, j, 9] +
        lin_min[9]*m[:static_lambda_box_cover2][i, j, 10] == 0
        )
    )
    @constraint(m, static_piecewise_lambda_constraint2[i=1:length(boxes), j=1:length(boxes); i!=j],
        (
        -m[:piecewise2][i, j] +
        lin_min_2[1]*m[:static_lambda_box_cover2][i, j, 2] +
        lin_min_2[2]*m[:static_lambda_box_cover2][i, j, 3] +
        lin_min_2[3]*m[:static_lambda_box_cover2][i, j, 4] +
        lin_min_2[4]*m[:static_lambda_box_cover2][i, j, 5] +
        lin_min_2[5]*m[:static_lambda_box_cover2][i, j, 6] +
        lin_min_2[6]*m[:static_lambda_box_cover2][i, j, 7] +
        lin_min_2[7]*m[:static_lambda_box_cover2][i, j, 8] +
        lin_min_2[8]*m[:static_lambda_box_cover2][i, j, 9] +
        lin_min_2[9]*m[:static_lambda_box_cover2][i, j, 10] == 0
        )
    )

    # Finding total partial stability of one box start. (length(boxes)-1 as we do not want to get support from sentinel).
    # Multiply by if, as if tells whether i,j is to be considered.
    @expression(m, pbs_box_ij[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if_on_top][i, j] * (m[:piecewise1][i, j] - m[:piecewise2][i, j])
    )
    @expression(m, pbs_box_i[i=1:length(boxes)],
        sum([m[:pbs_box_ij][i, j] for j in 1:length(boxes) if i!=j])
    )

    return m
end

end


text:
type PersistenceOperationResult{T}
	# the result of a persistence operation
	issuccessful::Bool
	errorencountered::Bool

	successfulset::Set{T}
	failedset::Set{T}

	lastitemattempted::T

	errormessage::String

	function PersistenceOperationResult()
		r = new()

		r.issuccessful = false
		r.errorencountered = false

		r.successfulset = Set{T}()
		r.failedset = Set{T}()

		return r
	end
end


text:
@testset "content" begin
	doc = read("../src/content.html",String)
	@test occursin("""<script type="text/javascript">""", doc)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using BaseTestNext

# test file to test testing

# Test @test
@test true
@test 1 == 1
@test 1 != 2
@test strip("\t  hi   \n") == "hi"
@test strip("\t  this should fail   \n") != "hi"

a = Array(Float64, 2, 2, 2, 2, 2)
a[1,1,1,1,1] = 10
@test a[1,1,1,1,1] == 10
@test a[1,1,1,1,1] != 2

@test rand() != rand()

# Test printing of result types
# Pass - constant
@test contains(sprint(show, @test true), "Expression: true")
# Pass - expression
@test contains(sprint(show, @test 10 == 2*5), "Evaluated: 10 == 10")
@test contains(sprint(show, @test !false), "Expression: !false")
# Pass - exception
@test contains(sprint(show, @test_throws ErrorException error()),
                "Thrown: ErrorException")

# Test printing of a TestSetException
tse_str = sprint(show, BaseTestNext.TestSetException(1,2,3))
@test contains(tse_str, "1 passed")
@test contains(tse_str, "2 failed")
@test contains(tse_str, "3 errored")

@test BaseTestNext.finish(BaseTestNext.FallbackTestSet()) != nothing

OLD_STDOUT = STDOUT
catch_out = IOStream("")
rd, wr = redirect_stdout()

# Check that the fallback test set throws immediately
@test_throws ErrorException (@test 1 == 2)

@testset "no errors" begin
    @test true
    @test 1 == 1
end

# Test entirely empty test set
@testset "outer" begin
    @testset "inner" begin
    end
end

try

@testset "outer" begin
    @testset "inner1" begin
        @test true
        @test false
        @test 1 == 1
        @test 2 == :foo
        @test 3 == 3
        @testset "d" begin
            @test 4 == 4
        end
        @testset begin
            @test :blank != :notblank
        end
    end
    @testset "inner1" begin
        @test 1 == 1
        @test 2 == 2
        @test 3 == :bar
        @test 4 == 4
        @test_throws ErrorException 1+1
        @test_throws ErrorException error()
        @test_throws RemoteException error()
        @testset "errrrr" begin
            @test "not bool"
            @test error()
        end

        error("exceptions in testsets should be caught")
        @test 1 == 1 # this test will not be run
    end

    @testset "loop with desc" begin
        @testset "loop1 $T" for T in (Float32, Float64)
            @test 1 == T(1)
        end
    end
    @testset "loops without desc" begin
        @testset for T in (Float32, Float64)
            @test 1 == T(1)
        end
        @testset for T in (Float32, Float64), S in (Int32,Int64)
            @test S(1) == T(1)
        end
    end
    srand(123)
    @testset "some loops fail" begin
        @testset for i in 1:5
            @test i <= rand(1:10)
        end
        # should add 3 errors and 3 passing tests
        @testset for i in 1:6
            iseven(i) || error("error outside of test")
            @test true # only gets run if the above passed
        end
    end
end
    # These lines shouldn't be called
    redirect_stdout(OLD_STDOUT)
    error("No exception was thrown!")
catch ex

    @test isa(ex, BaseTestNext.TestSetException)
    @test ex.pass  == 24
    @test ex.fail  == 6
    @test ex.error == 6
end

# Test @test_approx_eq
# TODO
@test isapprox(.1+.1+.1, .3)
@test !isapprox(.1+.1+.1, .4)

@test_throws ErrorException BaseTestNext.test_approx_eq(ones(10),ones(11),1e-8,"a","b")
@test_throws ErrorException BaseTestNext.test_approx_eq(ones(10),zeros(10),1e-8,"a","b")

# Test @test_approx_eq_eps
# TODO

ts = @testset "@testset should return the testset" begin
    @test true
end
@test typeof(ts) == BaseTestNext.DefaultTestSet
@test typeof(ts.results[1]) == BaseTestNext.Pass

tss = @testset "@testset/for should return an array of testsets: $i" for i in 1:3
    @test true
end
@test length(tss) == 3
@test typeof(tss[1]) == BaseTestNext.DefaultTestSet
@test typeof(tss[1].results[1]) == BaseTestNext.Pass

# now we're done running tests with DefaultTestSet so we can go back to STDOUT
redirect_stdout(OLD_STDOUT)

# import the methods needed for defining our own testset type
import BaseTestNext: record, finish
using BaseTestNext: get_testset_depth, get_testset
using BaseTestNext: AbstractTestSet, Result, Pass, Fail, Error
immutable CustomTestSet <: BaseTestNext.AbstractTestSet
    description::AbstractString
    foo::Int
    results::Vector
    # constructor takes a description string and options keyword arguments
    CustomTestSet(desc; foo=1) = new(desc, foo, [])
end

record(ts::CustomTestSet, child::AbstractTestSet) = push!(ts.results, child)
record(ts::CustomTestSet, res::Result) = push!(ts.results, res)
function finish(ts::CustomTestSet)
    # just record if we're not the top-level parent
    if get_testset_depth() > 0
        record(get_testset(), ts)
    end
    ts
end

ts = @testset CustomTestSet "Testing custom testsets" begin
    # this testset should inherit the parent testset type
    @testset "custom testset inner 1" begin
        @test true
        @test false
        @test error("this error will be reported as an error")
        @test_throws ErrorException nothing
        @test_throws ErrorException error("this error is a success")
    end
    # this testset has its own testset type
    @testset CustomTestSet foo=4 "custom testset inner 2" begin
        # this testset should inherit the type, but not the argument. If a particular
        # testset type wants inheritance behavior they should implement it themselves
        # using get_testset() in the constructor
        @testset "custom testset inner 2 inner 1" begin
            @test true
        end
        # make sure the RHS can use computed values, also tests options without
        # specifying the testset type
        @testset foo=(1+2) "custom testset inner 2 inner 2" begin
            @test true
        end
    end
end

@test typeof(ts) == CustomTestSet
@test ts.foo == 1
@test ts.description == "Testing custom testsets"
@test typeof(ts.results[1]) == CustomTestSet
@test ts.results[1].description == "custom testset inner 1"
@test ts.results[1].foo == 1
@test typeof(ts.results[1].results[1]) == Pass
@test typeof(ts.results[1].results[2]) == Fail
@test typeof(ts.results[1].results[3]) == Error
@test typeof(ts.results[1].results[4]) == Fail
@test typeof(ts.results[1].results[5]) == Pass

@test typeof(ts.results[2]) == CustomTestSet
@test ts.results[2].description == "custom testset inner 2"
@test ts.results[2].foo == 4
@test typeof(ts.results[2].results[1]) == CustomTestSet
@test ts.results[2].results[1].foo == 1
@test typeof(ts.results[2].results[1].results[1]) == Pass
@test typeof(ts.results[2].results[2]) == CustomTestSet
@test ts.results[2].results[2].foo == 3

# test custom testset types on testset/for
tss = @testset CustomTestSet foo=3 "custom testset $i" for i in 1:6
    @testset "inner testset $i-$j" for j in 1:3
        @test iseven(i + j)
    end
    # make sure a testset within a testset/for works
    @testset "inner testset $i" begin
        @test iseven(i)
    end
end

for i in 1:6
    @test typeof(tss[i]) == CustomTestSet
    @test tss[i].foo == 3
    for j in 1:3
        @test typeof(tss[i].results[j]) == CustomTestSet
        @test tss[i].results[j].foo == 1
        @test typeof(tss[i].results[j].results[1]) == (iseven(i+j) ? Pass : Fail)
    end
    @test typeof(tss[i].results[4]) == CustomTestSet
    @test typeof(tss[i].results[4].results[1]) == (iseven(i) ? Pass : Fail)
end