text:
"""
Forecast the state of the process at the next time step.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x : AbstractMvNormal.  Current state estimate.
- u : Vector, optional.  External input.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
	model is not time-varying, it will be ignored.)
#### Returns
- MvNormal distribution representing the forecast and its associated uncertainty.
"""
function predict(m::AbstractGaussianSSM, x::AbstractMvNormal;
		u::Vector=zeros(m.nu), t::Real=0.0)
    F = process_matrix(m, x, t)
    CI = control_input(m, u, t)
    return MvNormal(F * mean(x) + CI, F * cov(x) * F' + m.V(t))
end

"""
Observe the state process with uncertainty.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x : AbstractMvNormal.  Current state estimate.
- u : Vector, optional.  External input.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
model is not time-varying, it will be ignored.)
#### Returns
- MvNormal distribution, representing the probability of recording any
particular observation data.
"""
function observe(m::AbstractGaussianSSM, x::AbstractMvNormal, t::Real=0.0)
	G = observation_matrix(m, x, t)
	return MvNormal(G * mean(x), G * cov(x) * G' + m.W(t))
end

function _filter{T}(m::AbstractGaussianSSM, y::Array{T}, x0::AbstractMvNormal,
		u::Array{T}, times::Vector{T}, filter::AbstractKalmanFilter)
	x_filtered = Array(AbstractMvNormal, size(y, 2))
	loglik = 0.0
	x_pred = predict(m, x0, u=u[:, 1], t=times[1])
	x_filtered[1] = update(m, x_pred, y[:, 1], filter, 1)
	loglik = logpdf(x_filtered[1], mean(x_pred))
	if ! any(isnan(y[:, 1]))
		loglik += logpdf(observe(m, x_filtered[1]), y[:,1])
	end
	for i in 2:size(y, 2)
		x_pred = predict(m, x_filtered[i-1], u=u[:, i], t=times[i])
		# Check for missing values in observation
		if any(isnan(y[:, i]))
			x_filtered[i] = x_pred
		else
			x_filtered[i] = update(m, x_pred, y[:, i], filter, i)
			loglik += logpdf(observe(m, x_filtered[i]), y[:, i])
		end
		loglik += logpdf(x_pred, mean(x_filtered[i]))
	end
	return FilteredState(y, x_filtered, u, times, loglik, false)
end

"""
Given a process/observation model and a set of data, estimate the states of the
hidden process at each time step.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- y : Array of data.  Each column is a set of observations at one time, while
each row is one observed variable.  Can contain NaNs to represent missing data.
- x0 : AbstractMvNorm. Multivariate normal distribution representing our estimate
of the hidden state and our uncertainty about that estimate before our first observation.
(i.e., if the first data arrives at t=1, this is where we think the process is at t=0.)
- u : Matrix, optional.  Array of external inputs.
- t : Real, optional. Time supplied to the time-varying model matrices (if the
model is not time-varying, it will be ignored.)
#### Returns
- A FilteredState object holding the estimates of the hidden state.

#### Notes
This function only does forward-pass filtering--that is, the state estimate at time
t incorporates data from 1:t, but not from t+1:T.  For full forward-and-backward
filtering, run `smooth` on the FilteredState produced by this function.
"""
function filter{T}(m::LinearGaussianSSM, y::Array{T}, x0::AbstractMvNormal;
		filter::LinearKalmanFilter=KF(), u::Matrix{T}=zeros(m.nu, size(y, 2)),
		times::Vector{T}=zeros(size(y, 2)))
	return _filter(m, y, x0, u, times, filter)
end

function filter{T}(m::NonlinearGaussianSSM, y::Array{T}, x0::AbstractMvNormal;
		filter::NonlinearKalmanFilter=EKF(), u::Matrix{T}=zeros(m.nu, size(y, 2)),
		times::Vector{T}=zeros(size(y, 2)))
	return _filter(m, y, x0, u, times, filter)
end

"""
Given a process/observation model and a set of data, estimate the states of the
hidden process at each time step using forward and backward passes.

This function has two methods.  It can take either an AbstractGaussianSSM and a the
FilteredState produced by a previous call to `filter`, or it can take a model, data,
set, and initial state estimate, in which case it does the forward and backward
passes at once.

"""
function _smooth{T}(m::AbstractGaussianSSM, fs::FilteredState{T})
	n = size(fs.observations, 2)
	smooth_dist = Array(AbstractMvNormal, n)
	smooth_dist[end] = fs.state[end]
	if ! any(isnan(fs.observations[:, 1]))
		loglik = logpdf(observe(m, smooth_dist[end], fs.times[end]),
			fs.observations[:, end])
	end
	for i in (n - 1):-1:1
		state_pred = predict(m, fs.state[i], u=fs.input[:, i], t=fs.times[i])
		P = cov(fs.state[i])
		F = process_matrix(m, fs.state[i], fs.times[i])
		J = P * F' * inv(cov(state_pred))
		x_smooth = mean(fs.state[i]) + J *
			(mean(smooth_dist[i+1]) - mean(state_pred))
		P_smooth = P + J * (cov(smooth_dist[i+1]) - cov(state_pred)) * J'
		smooth_dist[i] = MvNormal(x_smooth, P_smooth)
		loglik += logpdf(predict(m, smooth_dist[i], u=fs.input[:, i], t=fs.times[i]),
			mean(smooth_dist[i+1]))
		if ! any(isnan(fs.observations[:, i]))
			loglik += logpdf(observe(m, smooth_dist[i], fs.times[i]), fs.observations[:, i])
		end
	end
	return FilteredState(fs.observations, smooth_dist, fs.input, fs.times, loglik, true)
end

function smooth{T}(m::AbstractGaussianSSM, fs::FilteredState{T})
	return _smooth(m, fs)
end

## Linear Gaussian method
function smooth{T}(m::LinearGaussianSSM, y::Array{T}, x0::AbstractMvNormal;
		u::Matrix{T}=zeros(m.nu, size(y, 2)), filter::LinearKalmanFilter=KF())
	fs = filter(m, y, x0, u=u, filter=filter)
	return _smooth(m, fs)
end

## Nonlinear Gaussian method
function smooth{T}(m::NonlinearGaussianSSM, y::Array{T}, x0::AbstractMvNormal;
		u::Matrix{T}=zeros(m.nu, size(y, 2)), filter::NonlinearKalmanFilter=EKF())
	fs = filter(m, y, x0, u=u, filter=filter)
	return _smooth(m, fs)
end

"""
Simulate a realization of a state-space process and observations of it.

#### Parameters
- m : AbstractGaussianSSM.  Model of the state evolution and observation processes.
- x0 : AbstractMvNorm. Distribution of the state process's starting value.

#### Returns
- (x, y) : Tuple of the process and observation arrays.
"""
function simulate(m::AbstractGaussianSSM, n::Int64, x0::AbstractMvNormal)
	F = process_matrix(m, x0, 1)
	G = observation_matrix(m, x0, 1)
	x = zeros(size(F, 1), n)
	y = zeros(size(G, 1), n)
	x[:, 1] = rand(MvNormal(F * mean(x0), m.V(1)))
	y[:, 1] = rand(MvNormal(G * x[:, 1], m.W(1)))
	for i in 2:n
		F = process_matrix(m, x[:, i-1], i)
		G = observation_matrix(m, x[:, i-1], i)
		x[:, i] = F * x[:, i-1] + rand(MvNormal(m.V(i)))
		y[:, i] = G * x[:, i] + rand(MvNormal(m.W(i)))
	end
	return (x, y)
end

"""
Fit a state-space model to data using maximum likelihood.

"""
function fit(build_func, y, x0, params0, args...)

	function objective(params)
		ssm = build_func(params)
		return -loglikelihood(smooth(m, y, x0))
	end
	fit = optimize(objective, params0, args...)
end


text:
"""
	RandomVariable(dist::Sampleable{Univariate}, name::Symbol)

Defines a random variable, with a univariate distribution from Distributions.jl and a name.

# Examples

```jldoctest
julia> RandomVariable(Normal(), :x)
RandomVariable(Normal{Float64}(μ=0.0, σ=1.0), :x)

julia> RandomVariable(Exponential(1), :x)
RandomVariable(Exponential{Float64}(θ=1.0), :x)
```
"""
struct RandomVariable <: RandomUQInput
    dist::Sampleable{Univariate}
    name::Symbol
end

"""
	sample(rv::RandomVariable, n::Integer)

Generates n samples from a random variable. Returns a DataFrame.

# Examples

```jldoctest
x = RandomVariable(Normal(), :x)
sample(x, 3)

# output

3×1 DataFrame
│ Row │ x          │
│     │ Float64    │
├─────┼────────────┤
│ 1   │ 1.00763    │
│ 2   │ -0.643618  │
│ 3   │ -0.0826132 │
```

See also: [`RandomVariable`](@ref)
"""
function sample(rv::RandomVariable, n::Integer=1)
    return DataFrame(rv.name => rand(rv.dist, n))
end

function to_physical_space!(rv::RandomVariable, x::DataFrame)
    x[!, rv.name] = quantile.(rv.dist, cdf.(Normal(), x[:, rv.name]))
    return nothing
end

function to_standard_normal_space!(rv::RandomVariable, x::DataFrame)
    x[!, rv.name] = quantile.(Normal(), cdf.(rv.dist, x[:, rv.name]))
    return nothing
end

mean(rv::RandomVariable) = DataFrame(rv.name => Distributions.mean(rv.dist))
mean(rvs::Array{RandomVariable}) = mapreduce(mean, hcat, rvs)

dimensions(rv::RandomVariable) = 1

ext: jl

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
function ∇ₛx∇ₛx end
function ∇ₛxμⁱ∇ₛx end
function ∇x∇x end
function ∇xμⁱ∇x end

function (sim::Simulation)(::typeof(∇x∇x))
    sim(∇x, d̂) * sim(∇x, p̂)
end

function (sim::Simulation)(::typeof(∇xμⁱ∇x))
    μi =  sim(convert_to_diagonal_matrix,:μᵣ)
    sim(∇x, d̂) *  μi  * sim(∇x, p̂)
end

function (sim::Simulation)(::typeof(∇ₛx∇ₛx))
    sim(∇ₛx, d̂) * sim(∇ₛx, p̂)
end

function (sim::Simulation)(::typeof(∇ₛxμⁱ∇ₛx))
    μi =  sim(convert_to_diagonal_matrix,:μᵣ)
    sim(∇ₛx, d̂) *  μi  * sim(∇ₛx, p̂)
end


text:
using Test
using Neo4j
using Graphs
using ProtoBuf
using JSON
using CloudGraphs
# using LibBSON
using Mongoc

import Base: convert

# Have we loaded the library?
@test isdefined(Main, :CloudGraphs) == true
@test typeof(CloudGraphs) == Module

# Testing type registration
mutable struct DataTest
  matrix::Array{Float64, 2}
  string::AbstractString #ASCIIString
  boolmatrix::Array{Int32,2}
  DataTest() = new()
  DataTest(m,s,b) = new(m,s,b)
end
mutable struct PackedDataTest
  vecmat::Vector{Float64}
  matrows::Int64
  string::AbstractString #ASCIIString
  boolvecmat::Array{Int32,1}
  boolmatrows::Int64
  PackedDataTest() = new()
  PackedDataTest(m,i1,s,b,i2) = new(m[:],i1,s,b[:],i2)
  PackedDataTest(d::DataTest) = new(d.matrix[:],
                                  size(d.matrix,1),
                                  d.string,
                                  d.boolmatrix[:],
                                  size(d.boolmatrix,1))
end

function convert(::Type{PackedDataTest}, d::DataTest) # encoder
  return PackedDataTest(d)
end
function convert(T::Type{DataTest}, d::PackedDataTest) # decoder
  r1 = d.matrows
  c1 = floor(Int,length(d.vecmat)/r1)
  M1 = reshape(d.vecmat,r1,c1)
  r2 = d.matrows
  c2 = floor(Int,length(d.boolvecmat)/r2)
  M2 = reshape(d.boolvecmat,r2,c2)
  return DataTest(M1,d.string,M2)
end

# Highly simplified packers and unpackers for testing.
function testEncodePackedType(a::DataTest)
  return PackedDataTest(a)
end
function testGetpackedtype(a::Any)
  return PackedDataTest()
end
function testDecodePackedType(a::Any, b::Any)
  return convert(DataTest, a)
end

# Defaults
if !haskey(ENV, "NEO4JUN")
    ENV["NEO4JUN"] = "neo4j"
end
if !haskey(ENV, "NEO4JPW")
    ENV["NEO4JPW"] = "marine"
end
if !haskey(ENV, "MONGOUN")
    ENV["MONGOUN"] = ""
end
if !haskey(ENV, "MONGOPW")
    ENV["MONGOPW"] = ""
end
configuration = CloudGraphs.CloudGraphConfiguration("localhost", 7474, ENV["NEO4JUN"], ENV["NEO4JPW"], "localhost", 27017, false, ENV["MONGOUN"], ENV["MONGOPW"]);
cloudGraph = connect(configuration, testEncodePackedType, testGetpackedtype, testDecodePackedType);
println("Success!");


text:
module Julog

include("structs.jl")
include("parse.jl")
include("utils.jl")
include("main.jl")

export Const, Var, Compound, Term, Clause, Subst, ClauseTable
export get_args, is_ground, substitute, has_subterm, find_subterms
export to_nnf, to_cnf, to_dnf, flatten_conjs, flatten_disjs, deuniversalize
export eval_term, unify, resolve, derivations, derive, fwd_chain, bwd_chain
export regularize_clauses, index_clauses, deindex_clauses, retrieve_clauses
export insert_clauses!, insert_clauses, subtract_clauses!, subtract_clauses
export parse_prolog, write_prolog
export @julog, @prolog, @varsub

end # module


text:
module SimplePlots

  using UUIDs
  using JSON

  using Printf
  using DataStructures

  using IJulia.CommManager
  using Observables

  using IJulia

  abstract type AbstractPlot end

  include("simple_plot.jl")
  include("composite_plot.jl")

  include("show.jl")
  include("render.jl")

  include("text.jl")
  include("parse_layout.jl")

  include("plots/index.jl")
  include("interact/index.jl")

  include("init.jl")
  include("strip_escape.jl")

end


text:
##############################################################################
#
# The hash trick: use a hash function instead of a lexicon to determine the
# columns of a DocumentTermMatrix-like encoding of the data
#
# To produce a DTM for a Corpus for which we do not have an existing lexicon,
# we need someway to map the terms from each document into column indices.
#
# We use the now standard "Hash Trick" (CITATION NEEDED), in which we
# hash strings and then reduce the resulting integers modulo N, which
# defines the numbers of columns we want our DTM to have. This amounts to
# doing a non-linear dimensionality reduction with low probability that similar
# terms hash to the same dimension.
#
# To make things easier, we wrap Julia's hash functions in a new type,
# TextHashFunction, which maintains information about the desired cardinality
# of the hashes.
#
##############################################################################

type TextHashFunction
    hash_function::Function
    cardinality::Int
end

TextHashFunction(cardinality::Int) = TextHashFunction(hash, cardinality)

TextHashFunction() = TextHashFunction(hash, 100)

cardinality(h::TextHashFunction) = h.cardinality

function index_hash(s::AbstractString, h::TextHashFunction)
    return @compat(Int(rem(h.hash_function(s), h.cardinality))) + 1
end


text:
################################################################################
function init(::Type{Metropolis}, config::ConfigMetropolis, objective::Objective, difftune)
    return Metropolis(ModelWrappers.ℓDensityResult(objective))
end

function infer(
    _rng::Random.AbstractRNG,
    diagnostics::Type{AbstractDiagnostics},
    kernel::Metropolis,
    model::ModelWrapper,
    data::D,
) where {D}
    return DiagnosticsMetropolis{model.info.reconstruct.default.output}
end

################################################################################
# Export
export infer, init


text:
import Base.intersect

complex_empty = Set{Complex{Float64}}()
identical_warn = "Intersection of identical clines returns the empty set"

"""
    intersect(C::Cline,D::Cline)
    C ∩ D

Find the points common to the clines `C` and `D` as a set of complex numbers.

Note that if `C==D` a warning is issued and the emptyset is returned.
"""
function intersect(C::Circle, D::Circle)::Set{Complex{Float64}}
    if C == D
        @warn identical_warn
        return complex_empty
    end

    z = center(C)
    r = radius(C)

    w = center(D)
    s = radius(D)

    d = abs(z - w)

    if d > r + s || d + s < r || d + r < s
        return Set{Complex{Float64}}()
    end

    θ = acos((r^2 + d^2 - s^2) / (2 * r * d))
    u = (w - z) / abs(w - z)
    uu = u * exp(θ * im)
    a = z + uu * r
    uu = u * exp(-θ * im)
    b = z + uu * r
    return Set([a, b])
end

function intersect(L::Line, M::Line)::Set{Complex{Float64}}
    if L == M
        @warn identical_warn
        return complex_empty
    end

    # check if parallel
    m1 = slope(L)
    m2 = slope(M)
    if (isinf(m1) && isinf(m2)) || (abs(m1 - m2) < get_tolerance())
        return complex_empty
    end

    a = L.a
    b = L.b

    ax, ay = reim(a)
    bx, by = reim(b)

    aa = M.a
    bb = M.b

    aax, aay = reim(aa)
    bbx, bby = reim(bb)

    MAT = zeros(2, 2)
    MAT[1, 1] = ax - bx
    MAT[1, 2] = bbx - aax
    MAT[2, 1] = ay - by
    MAT[2, 2] = bby - aay

    RHS = [bbx - bx, bby - by]
    st = MAT \ RHS

    s = st[1]
    z = s * L.a + (1 - s) * L.b
    return Set(z)
end

function intersect(C::Cline, D::Cline)::Set{Complex{Float64}}
    if C == D
        @warn identical_warn
        return complex_empty
    end

    F = LFT(C)  # map C to x-axis
    G = inv(F)
    DD = F(D)  # apply that to D and find where DD crosses the x-axis

    # if the images of D is a circle
    if isa(DD, Circle)
        x, y = reim(center(DD))
        r = radius(DD) + get_tolerance()

        if r < abs(y)
            return complex_empty = Set{Complex{Float64}}()
        end

        δ = sqrt(r^2 - y^2)
        a = G(x + δ)
        b = G(x - δ)
        if abs(a - b) < get_tolerance()  # if the points are too close, we say they're the same
            ab = (a + b) / 2
            return Set{Complex{Float64}}(ab)
        else
            return Set{Complex{Float64}}([a, b])
        end
    end

    # otherwise, DD, the image of D, is a line

    X = Line(0, 1)  # x-axis
    S = X ∩ DD
    return G.(S)

    # x, y = reim(DD.a)
    # xx, yy = reim(DD.b)

    # if abs(y - yy) < get_tolerance()  # DD is a horizontal line
    #     return complex_empty()
    # end

    # if abs(x - xx) < get_tolerance()  # DD is a vertical line
    #     z = G((x + xx) / 2)
    #     return Set{Complex{Float64}}(z)
    # end

    # slope = (y - yy) / (x - xx)
    # δ = y / slope
    # z = G(x - δ)
    # return Set{Complex{Float64}}(z)
end


text:
""" ExprModels

provides functionality of extracting and manipulating models from julia source code

"""
module ExprModels

include("parse.jl")
using SemanticModels.ExprModels.Parsers

using SemanticModels
import SemanticModels: model

export model, isexpr, invoke, callsites, structured, bodyblock,
       pusharg!, setarg!, funcarg, issome, head, isblock, isfunc,
       or, and, isexpr, iscall, isusing, isimport, funcarg, funcname

""" model(::AbstractModel, expr::Expr)
dig into the expression that describes a model and break it down into components. This allows you to construct a structured representation of the modeling problem at the expression level. Just like how julia modeling frameworks build structured representations of the problems in data structures. This version builds them at the expression level.
The first argument is the type you want to construct, the second argument is the expression that you want to analyze. For example
```
model(ExpODEModel, expr)::ExpODEModel
```
"""
function model(::Type{T}, expr::Expr) where T<:SemanticModels.AbstractModel
    error("NotImplemented: model(::$T,::Expr")
end

""" isexpr(x)
predicate for isa(x, Expr).
"""
isexpr(x) = isa(x, Expr)

function invoke(m::AbstractModel, args...)
    Mod = eval(m.expr)
    Base.invokelatest(Mod.main, args...)
end

"""    callsites(expr::Expr, name::Symbol)

extract the location where the function `name` is called in `expr`.
"""
function callsites(expr::Expr, name::Symbol)
    matches = Expr[]
    f(x::Any) = x
    f(x::Expr) = begin
        if x.head == :call && x.args[2] == name
            push!(matches, x)
            return x
        else
            return walk(x, f, g)
        end
    end
    g(x) = x
    walk(expr, f, g)
    return matches
end

"""    structured(func, var::Symbol, assign=true)

extract the expressions that use structuring/destructuring assignment
to name the components of `var`
"""
function structured(func::Expr, var::Symbol, assign=true)
    body = func.args[2]
    paramvec = Expr[]
    for line in body.args
        if isa(line, LineNumberNode)
            continue
        end
        if  assign && line.head == :(=) && inexpr(line.args[2], var)
            push!(paramvec, line)
        end
        if !assign && line.head == :(=) && inexpr(line.args[1], var)
            push!(paramvec, line)
        end
    end
    return paramvec
end

"""    bodyblock(expr::Expr)

get the array of args representing the body of a function definition.
"""
function bodyblock(expr::Expr)
    expr.head == :function || error("$expr is not a function definition")
    return expr.args[2].args
end

"""    pusharg!(expr::Expr, s::Symbol)

push a new argument onto the definition of a function.

See also [`argslist`](@ref), [`setarg!`](@ref)
"""
function pusharg!(ex::Expr, s::Symbol)
    ex.head == :function || error("ex is not a function definition")
    push!(argslist(ex), s)
    return ex
end

"""    setarg!(expr::Expr, s::Symbol)

replace the argument in a function call.

See also [`argslist`](@ref), [`pusharg!`](@ref)
"""
function setarg!(ex::Expr, old, new)
    ex.head == :call || error("ex is not a function call")
    for (i, x) in enumerate(ex.args)
        if x == old
            ex.args[i] = new
        end
    end
    return ex
end

function funcarg(ex::Expr)
    return ex.args[1].args[2]
end

include("exprs.jl")
include("Transformations.jl")
include("SimpleModels.jl")
include("ExpODEModels.jl")
include("ExpStateModels.jl")
include("MonomialRegressionModels.jl")

end


text:
# Type variable deconstruction patterns
module TypeVarDecons
using MLStyle
using MLStyle.Infras
using MLStyle.MatchCore
using MLStyle.Qualification
using MLStyle.TypeVarExtraction

export type_matching

macro type_matching(t, forall)
    quote
        NAME = mangle(mod)
        __T__ = $t
        __FORALL__ = $forall
        if !($any_constraint(__T__, __FORALL__))
            function (body)
                @format [body, tag, NAME, TARGET, __T__] quote
                    @inline __L__ function NAME(TARGET :: __T__) where {$(__FORALL__...)}
                        __T__ # if not put this here, an error would be raised : "local variable XXX cannot be used in closure declaration"
                        body
                    end
                    NAME(tag)
                end
            end
        else
            function (body)
                @format [body, tag, NAME, TARGET, __T__] quote
                    @inline __L__ function NAME(TARGET :: __T__) where {$(__FORALL__...)}
                        __T__
                        body
                    end
                    @inline __L__ function NAME(_)
                        failed
                    end
                    NAME(tag)
                end
            end
        end
    end |> esc
end

def_pattern(TypeVarDecons,
        predicate = x -> x isa Expr && x.head == :(::),
        rewrite = (tag, case, mod) ->
            let args   = (case.args..., ),
                TARGET = mangle(mod)
                function for_type(t)
                    @match t begin
                        :($typ where {$(tvars...)}) => (@type_matching typ tvars)
                        _ => @typed_as t

                    end
                end

                function f(args :: NTuple{2, Any})
                    pat, t = args
                    for_type(t) ∘ mk_pattern(TARGET, pat, mod)
                end

                function f(args :: NTuple{1, Any})
                    t = args[1]
                    for_type(t)
                end
                f(args)
            end)

end


text:
using Test
@test 1==1


text:
using Test, Random
using LinearAlgebra
using MATLAB
using Arpack

rng = MersenneTwister(123)
for n in [2, 5, 30, 100, 1000]
    P = randn(rng, n, n); P = (P + P')/2;
    q = randn(rng, n)/10
    C = randn(n, n); C = (C + C')/2;
    # Make C positive definite
    C -= 1.01*min(minimum(eigvals(C)), 0)*I
    r = [1e-4 1e-2 1 1000]
    eye = Matrix{Float64}(I, n, n)
    for i = 1:length(r)
        if rand(rng) < .1
            # Make P positive definite
            P = P + 1.1*abs(eigs(P, nev=1, which=:SR)[1][1])*I;
        elseif rand(rng) < .1
            # Make P negative definite
            P = P - 1.1*abs(eigs(P, nev=1, which=:LR)[1][1])*I;
        end
        if n < 30
            X, info = TRS.trs_small(P, q, r[i], C, compute_local=true)
        else
            X, info = TRS.trs(P, q, r[i], C, compute_local=true)
        end
        x_g = X[:, 1]
        x_matlab, λ_matlab = mxcall(:TRSgep, 2, P, q, C, r[i])
        @show norm(P*x_g + q + info.λ[1]*C*x_g), norm(P*x_matlab[:, 1] + q + λ_matlab*C*x_matlab[:, 1])
        # @show info.λ[1] + (dot(x_g, P*x_g) + dot(q,x_g))/r[i]^2
        # @show info.λ[1], λ_matlab
        str = "Ellipsoidal Trs - r:"*string(r[i])
        @testset "$str" begin
            @test info.λ[1] - λ_matlab <= 1e-6*λ_matlab
            if size(x_matlab, 2) > 1
                diff = min(norm(x_g - x_matlab[:, 1]), norm(x_g - x_matlab[:, 2]))
            else
                diff = norm(x_g - x_matlab)
            end
            @test diff <= 1e-3*r[i]
        end
    end
    # hard case
    λ_min, v, _ = eigs(-P, nev=1, which=:LR)
    v = v/norm(v)
    q = (I - v*v')*q
    if n < 0
        X, info = TRS.trs_small(P, q, r[end], C, compute_local=true)
    else
        X, info = TRS.trs(P, q, r[end], C, compute_local=true)
    end
    x_g = X[:, 1]
    x_matlab, λ_matlab = mxcall(:TRSgep, 2, P, q, C, r[end])
    @show norm(P*x_g + q + info.λ[1]*C*x_g), norm(P*x_matlab[:, 1] + q + λ_matlab*C*x_matlab[:, 1])
    # @show info.λ[1] + (dot(x_g, P*x_g) + dot(q,x_g))/r[end]^2
    # @show info.λ[1], λ_matlab
    @testset "Ellipsoidal Trs - hard case" begin
        @test info.λ[1] - λ_matlab <= 1e-6*λ_matlab
        if size(x_matlab, 2) > 1
            diff = min(norm(x_g - x_matlab[:, 1]), norm(x_g - x_matlab[:, 2]))
        else
            diff = norm(x_g - x_matlab)
        end
        @test diff <= 1e-3*r[end]
    end
end

nothing


text:
using BSplineKit
using LinearAlgebra
using QuadGK
using SpecialFunctions: airyai

# Solution of Airy differential equation using Galerkin method.
# Problem taken from Olver & Townsend SIAM Rev. 2013.
# The original equation is for a function `u` that satisfies inhomogeneous
# Dirichlet BCs at x = ±1.
# Here we perform the change of variables v(x) = u(x) - u₀(x), such that `v`
# satisfies *homogeneous* Dirichlet BCs.
function test_airy_equation()
    ε = 1e-4
    N = 100
    breaks = collect(range(-1, 1; length = N + 1))
    B = @inferred BSplineBasis(BSplineOrder(6), breaks)
    R = @inferred RecombinedBSplineBasis(Derivative(0), B)  # for Dirichlet BCs

    # Exact solution
    u_exact(x) = airyai(x / cbrt(ε))

    # Affine part of the solution, for change of variables.
    # This is the (unique) affine curve that satisfies the two inhomogeneous BCs.
    u_0(x) = ((1 - x) * u_exact(-1) + (1 + x) * u_exact(1)) / 2

    # This is to obtain coefficients from data at collocation points
    xs = @inferred collocation_points(B)
    Cfact = @inferred lu!(collocation_matrix(B, xs))

    # Variable coefficient operator: c(x) * v(x), with c(x) = x.
    # Note that we could expand c(x) = x using linear splines (k = 2), but we
    # use the same basis as for `v` for simplicity.
    # Here Mc is the linear Galerkin operator associated to multiplication by
    # c(x).
    Mc = let
        cs = @inferred Cfact \ xs  # coefficients for c(x) = x
        Tc = @inferred galerkin_tensor((R, R, B), Derivative.((0, 0, 0)))
        Mc = @inferred Tc * cs
        @test Mc isa BandedMatrix
        @test Mc ≈ Mc'  # the matrix is practically symmetric
        Hermitian(Mc)
    end

    # Second-order term: ε v''(x)
    # Taking advantage of integration by parts, we construct the matrix
    # M′[i, j] = ⟨ ϕ′[i] ϕ′[j] ⟩.
    # The boundary term is zero due to Dirichlet BCs.
    L = let
        L = @inferred galerkin_matrix(R, Derivative.((1, 1)))
        @test L isa Hermitian
        @test L.uplo == 'U'
        parent(L) .*= ε
        L
    end

    # We add the two LHS operators
    # We only need to add the upper triangular part, because of symmetry.
    parent(L) .+= UpperTriangular(Mc)

    # Factorisation
    # Apparently we can't do Cholesky because `Mc` is not SPD.
    # For LU, we can't pass a Hermitian matrix wrapper...
    Lfact = lu(BandedMatrix(L))

    # RHS resulting from change of variables
    f(x) = x * u_0(x)
    rhs = let
        fs = ldiv!(Cfact, f.(xs))     # coefficients of RHS term
        Mf = galerkin_matrix((R, B))  # RHS operator
        -Mf * fs
    end

    # Numerical solution
    vR = Lfact \ rhs  # coefficients of the solution in recombined basis R
    vB = @inferred recombination_matrix(R) * vR  # same in full basis B
    vsol = @inferred Spline(B, vB)  # solution as a spline
    usol(x) = vsol(x) + u_0(x)  # undo change of variables

    # Check error (this will depend on the choices of ε, N, k, …)
    err, _ = quadgk(x -> (usol(x) - u_exact(x))^2, -1, 1; rtol = 1e-8)
    @test err < 1e-6

    nothing
end

@testset "Airy equation" begin
    test_airy_equation()
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
Compiled templating language for Genie.
"""
module Flax

import Revise
import SHA, Reexport, OrderedCollections, Logging, FilePaths
import Genie, Genie.Configuration
Reexport.@reexport using HttpCommon

export @foreach, @vars, @yield
export partial, template

import Base.string
import Base.show
import Base.==
import Base.hash

include("HTMLRenderer.jl")
using .HTMLRenderer

include("JSONRenderer.jl")
using .JSONRenderer

const Html = @__MODULE__

const BUILD_NAME    = "FlaxViews"

init_task_local_storage() = (haskey(task_local_storage(), :__vars) || task_local_storage(:__vars, Dict{Symbol,Any}()))
init_task_local_storage()

task_local_storage(:__yield, "")

"""
    partial(path::String; context::Module = @__MODULE__, vars...) :: String

Renders (includes) a view partial within a larger view or layout file.
"""
function partial(path::String; context::Module = @__MODULE__, vars...) :: String
  for (k,v) in vars
    try
      task_local_storage(:__vars)[k] = v
    catch
      init_task_local_storage()
      task_local_storage(:__vars)[k] = v
    end
  end

  template(path, partial = true, context = context)
end

"""
    parseview(data::String; partial = false, context::Module = @__MODULE__) :: Function

Parses a view file, returning a rendering function. If necessary, the function is JIT-compiled, persisted and loaded into memory.
"""
function parseview(data::String; partial = false, context::Module = @__MODULE__) :: Function
  data_hash = hash(data)
  path = "Flax_" * string(data_hash)

  func_name = function_name(string(data_hash, partial)) |> Symbol
  mod_name = m_name(string(path, partial)) * ".jl"
  f_path = joinpath(Genie.config.path_build, BUILD_NAME, mod_name)
  f_stale = build_is_stale(f_path, f_path)

  if f_stale || ! isdefined(context, func_name)
    f_stale && build_module(string_to_flax(data, partial = partial), path, mod_name)

    return Base.include(context, joinpath(Genie.config.path_build, BUILD_NAME, mod_name))
  end

  getfield(context, func_name)
end

"""
    template(path::String; partial::Bool = true, context::Module = @__MODULE__) :: String

Renders a template file.
"""
@inline function template(path::String; partial::Bool = true, context::Module = @__MODULE__) :: String
  HTMLRenderer.get_template(path, partial = partial, context = context) |> Base.invokelatest
end

"""
    build_is_stale(file_path::String, build_path::String) :: Bool

Checks if the view template has been changed since the last time the template was compiled.
"""
function build_is_stale(file_path::String, build_path::String) :: Bool
  isfile(file_path) || return true

  file_mtime = stat(file_path).mtime
  build_mtime = stat(build_path).mtime
  status = file_mtime > build_mtime

  Genie.config.log_views && status && @debug("🚨  Flax view $file_path build $build_path is stale")

  status
end

"""
    registervars(vars...) :: Nothing

Loads the rendering vars into the task's scope
"""
@inline function registervars(vars...) :: Nothing
  init_task_local_storage()
  task_local_storage(:__vars, merge(Dict{Symbol,Any}(vars), task_local_storage(:__vars)))

  nothing
end

@inline function vars_signature() :: String
  task_local_storage(:__vars) |> keys |> collect |> sort |> string
end

"""
    function_name(file_path::String)

Generates function name for generated Flax views.
"""
@inline function function_name(file_path::String) :: String
  "func_$(SHA.sha1( relpath(isempty(file_path) ? " " : file_path) * vars_signature() ) |> bytes2hex)"
end

"""
    m_name(file_path::String)

Generates module name for generated Flax views.
"""
@inline function m_name(file_path::String) :: String
  string(SHA.sha1( relpath(isempty(file_path) ? " " : file_path) * vars_signature()) |> bytes2hex)
end

"""
    html_to_flax(file_path::String; partial = true) :: String

Converts a HTML document to a Flax document.
"""
@inline function html_to_flax(file_path::String; partial = true) :: String
  to_flax(file_path, parse_template, partial = partial)
end

"""
    string_to_flax(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts string view data to Flax code
"""
@inline function string_to_flax(content::String; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String
  to_flax(content, parse_string, partial = partial, f_name = f_name, prepend = prepend)
end

"""
    to_flax(input::String, f::Function; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "") :: String

Converts an input file to Flax code
"""
@inline function to_flax(input::String, f::Function; partial = true, f_name::Union{Symbol,Nothing} = nothing, prepend = "\n") :: String
  f_name = (f_name === nothing) ? function_name(string(input, partial)) : f_name

  string("function $(f_name)() \n",
          injectvars(),
          prepend,
          f(input, partial = partial),
          "\nend \n")
end

function injectvars() :: String
  output = ""
  for kv in task_local_storage(:__vars)
    output *= "$(kv[1]) = @vars($(repr(kv[1]))) \n"
  end

  output
end

function injectvars(context::Module) :: Nothing
  for kv in task_local_storage(:__vars)
    isdefined(context, Symbol(kv[1])) || Core.eval(context, Meta.parse("$(kv[1]) = @vars($(repr(kv[1])))"))
  end

  nothing
end

"""
    build_module(content::String, path::String, mod_name::String) :: String

Persists compiled Flax view data to file and returns the path
"""
function build_module(content::String, path::String, mod_name::String) :: String
  module_path = joinpath(Genie.config.path_build, BUILD_NAME, mod_name)

  isdir(dirname(module_path)) || mkpath(dirname(module_path))

  open(module_path, "w") do io
    write(io, "# $path \n\n")
    write(io, content)
  end

  module_path
end

"""
    read_template_file(file_path::String) :: String

Reads `file_path` template from disk.
"""
function read_template_file(file_path::String) :: String
  io = IOBuffer()
  open(file_path) do f
    for line in enumerate(eachline(f))
      print(io, parsetags(line), "\n")
    end
  end

  String(take!(io))
end

"""
    parse_template(file_path::String; partial = true) :: String

Parses a HTML file into Flax code.
"""
@inline function parse_template(file_path::String; partial::Bool = true) :: String
  parse(read_template_file(file_path), partial = partial)
end

"""
    parse_string(data::String; partial = true) :: String

Parses a HTML string into Flax code.
"""
@inline function parse_string(data::String; partial::Bool = true) :: String
  parse(parsetags(data), partial = partial)
end

@inline function parse(input::String; partial::Bool = true) :: String
  HTMLRenderer.parsehtml(input, partial = partial)
end

"""
    parsetags(line::Tuple{Int64,String}, strip_close_tag = false) :: String

Parses special Flax tags.
"""
@inline function parsetags(line::Tuple{Int64,String}) :: String
  parsetags(line[2])
end
@inline function parsetags(code::String) :: String
  code = replace(code, "<%"=>"""<script type="julia/eval">""")
  replace(code, "%>"=>"""</script>""")
end

"""
    register_elements() :: Nothing

Generated functions that represent Flax functions definitions corresponding to HTML elements.
"""
@inline function register_elements() :: Nothing
  for elem in HTMLRenderer.NORMAL_ELEMENTS
    register_normal_element(elem)
  end

  for elem in HTMLRenderer.VOID_ELEMENTS
    register_void_element(elem)
  end

  nothing
end

@inline function register_element(elem::Union{Symbol,String}, elem_type::Union{Symbol,String} = :normal; context = Flax) :: Nothing
  elem = string(elem)
  occursin('-', elem) && (elem = HTMLRenderer.denormalize_element(elem))

  elem_type == :normal ? register_normal_element(elem) : register_void_element(elem)
end

function register_normal_element(elem::Union{Symbol,String}; context = Flax) :: Nothing
  Core.eval(context, """
    function $elem(f::Function, args...; attrs...) :: HTMLString
      \"\"\"\$(HTMLRenderer.normal_element(f, "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  Core.eval(context, """
    function $elem(children::Union{String,Vector{String}} = "", args...; attrs...) :: HTMLString
      \"\"\"\$(HTMLRenderer.normal_element(children, "$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  nothing
end

function register_void_element(elem::Union{Symbol,String}; context = Flax) :: Nothing
  Core.eval(context, """
    function $elem(args...; attrs...) :: HTMLString
      \"\"\"\$(HTMLRenderer.void_element("$(string(elem))", [args...], Pair{Symbol,Any}[attrs...]))\"\"\"
    end
  """ |> Meta.parse)

  nothing
end

"""
    @foreach(f, arr)

Iterates over the `arr` Array and applies function `f` for each element.
The results of each iteration are concatenated and the final string is returned.

## Examples

@foreach(@vars(:translations)) do t
  t
end
"""
macro foreach(f, arr)
  e = quote
    isempty($(esc(arr))) && return ""

    mapreduce(*, $(esc(arr))) do _s
      $(esc(f))(_s)
    end
  end

  quote
    Core.eval($__module__, $e)
  end
end

register_elements()

"""
    @vars

Utility macro for accessing view vars
"""
macro vars()
  :(task_local_storage(:__vars))
end

"""
    @vars(key)

Utility macro for accessing view vars stored under `key`
"""
macro vars(key)
  :(task_local_storage(:__vars)[$key])
end

"""
    @vars(key, value)

Utility macro for setting a new view var, as `key` => `value`
"""
macro vars(key, value)
  quote
    try
      task_local_storage(:__vars)[$key] = $(esc(value))
    catch
      init_task_local_storage()
      task_local_storage(:__vars)[$key] = $(esc(value))
    end
  end
end

"""
    @yield

Outputs the rendering of the view within the template.
"""
macro yield()
  quote
    try
      task_local_storage(:__yield)
    catch
      task_local_storage(:__yield, "")
    end
  end
end
macro yield(value)
  :(task_local_storage(:__yield, $value))
end

function el(; vars...)
  OrderedCollections.OrderedDict(vars)
end

"""
    preparebuilds() :: Bool

Sets up the build folder and the build module file for generating the compiled views.
"""
function preparebuilds(subfolder = BUILD_NAME) :: Bool
  build_path = joinpath(Genie.config.path_build, subfolder)
  isdir(build_path) || mkpath(build_path)

  true
end

function purgebuilds(subfolder = BUILD_NAME) :: Bool
  rm(joinpath(Genie.config.path_build, subfolder), force = true, recursive = true)

  true
end

function changebuilds(subfolder = BUILD_NAME) :: Bool
  Genie.config.path_build = Genie.Configuration.buildpath()
  preparebuilds()
end

end


text:
# DirectedGraph

"Representing a fixed graph with edges stored in an array"
abstract type FixedGraph end

"Representing a fixed directed graph with directed edges stored in an array"
abstract type FixedDirectedGraph <: FixedGraph end

"""
    size(graph::FixedGraph)

Return the size of the graph `graph`
"""
size(graph::FixedGraph) = graph.size

"Return the array of edges of the graph"
edges(graph::FixedGraph) = graph.edges

"Return the nodes of the graph as an iterator"
nodes(graph::FixedGraph) = 1:graph.size

"Representing a fixed tree with directed edges stored in an array"
mutable struct FixedTreeGraph <: FixedDirectedGraph
    size::Int
    edges::Vector{Tuple{Int, Int}}
    positions::Union{Array{Int, 2}, Nothing}
end

"""
    FixedTreeGraph(size::Int)

Construct a `FixedTreeGraph` with no edges of `size`.
"""
FixedTreeGraph(size::Int) = FixedTreeGraph(size, [], nothing)

"""
    FixedTreeGraph(size::Int, edges::Array{Tuple{Int, Int}})

Construct a `FixedTreeGraph` with edges `edges` of size `size`.
"""
FixedTreeGraph(size::Int, edges::Array{Tuple{Int, Int}}) = FixedTreeGraph(size, edges, nothing)

"""
    addedge!(tree::FixedGraph, from, to)

Add an edge  `(from, to)` in `tree`.
"""
function addedge!(tree::FixedGraph, from, to)
    push!(edges(tree), (from, to))
end

"""
    positions(tree::FixedTreeGraph)::Array{Int, 2}

First sort the edges in the tree and then return an array containing the position and length of the
segment in the edge array corresponding to each node.
"""
function positions(tree::FixedTreeGraph)::Array{Int, 2}
    # Do not repeat computing this
    if tree.positions !== nothing
        return tree.positions
    end

    # Sort the edges by starting node
    sortedge!(tree)

    # The first number indicates where the segment belonging to a node starts.
    # The second number is how many edges this node has.
    position_arr = zeros(Int, (size(tree),2))
    current_node = 0
    for i in 1:(size(tree)-1)
        u, v = tree.edges[i]
        if u!= current_node
            current_node = u
            position_arr[current_node,1]=i
        end
        position_arr[current_node, 2]+=1
    end
    tree.positions = position_arr
end

"""
    sortedge!(tree::FixedTreeGraph)

Sort the edges in `tree` according to starting nodes.
"""
function sortedge!(tree::FixedTreeGraph)
    vertex_count = zeros(Int, size(tree))

    edge_array = edges(tree)

    # Count the number of edges
    for edge in edge_array
        @inbounds vertex_count[edge[1]] += 1
    end

    # Transform to cumulative sum
    for node in 2:size(tree)
        @inbounds vertex_count[node] += vertex_count[node-1]
    end

    # copy edges to the right position
    edge_sorted = Vector{Tuple{Int, Int}}(undef, length(edge_array))
    for edge in edge_array
        @inbounds from = edge[1]
        @inbounds position = vertex_count[from]
        @inbounds edge_sorted[position] = edge
        @inbounds vertex_count[from] -= 1
    end

    tree.edges = edge_sorted
end

function iterate(tree::FixedTreeGraph)
    root_degree = positions(tree)[1, 2]

    stack = Stack{Tuple{Int, Int}}()
    if root_degree > 0
        push!(stack, (1, 0))
    end

    return (root_degree, stack)
end

function iterate(tree::FixedTreeGraph, stack)
    stack
    position_array = positions(tree)

    if isempty(stack)
        return nothing
    end

    parent_node, child_visited = pop!(stack)

    @inbounds next_edge = edges(tree)[position_array[parent_node]+child_visited]
    @inbounds next_node = next_edge[2]

    # Do we need to put parent back to the stack?
    @inbounds parent_degree = position_array[parent_node, 2]
    child_visited += 1
    if child_visited < parent_degree
        push!(stack, (parent_node, child_visited))
    end

    # Put the new node in the stack
    @inbounds current_degree = position_array[next_node, 2]
    if current_degree >= 1
        push!(stack, (next_node, 0))
    end

    return (current_degree, stack)
end

degrees(tree::FixedTreeGraph) = tree

length(tree::FixedTreeGraph) = size(tree)
eltype(tree::FixedTreeGraph) = Int


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

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

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

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

include("../../../src/Sparso.jl")
using Sparso

include("utils.jl")

# TEMPORARY: this should be removed in future once structure analysis can
# recognize that L and U are part of A
# Copied from context-test2-ilu0.jl to avoid affecting other tests if it were changed.
function pcg_symgs_ilu0(x, A, b, tol, maxiter)
    # These are the only 3 properties we need to discover automatically.
    # We also require SA_CONST_STRUCTURED for A, but that has already been discovered.
    set_matrix_property(Dict(
        :A => SA_SYMM_VALUED
        )
    )

    Ap = zeros(size(A, 1))

    total_time = time()

    trsv_time = 0.
    spmv_time = 0.
    blas1_time = 0.

    L, U = Sparso.ilu(A)

    spmv_time -= time()
    r = b - A * x
    spmv_time += time()

    blas1_time -= time()
    normr0 = norm(r)
    z = copy(r)
    blas1_time += time()

    rel_err = 1

    trsv_time -= time()
    Base.SparseMatrix.fwdTriSolve!(L, z)
    Base.SparseMatrix.bwdTriSolve!(U, z)
    trsv_time += time()

    blas1_time -= time()
    p = copy(z) #NOTE: do not write "p=z"! That would make p and z aliased (the same variable)
    rz = dot(r, z)
    blas1_time += time()

    # For statistics only. We get these values so that L, U and A won't be live out of
    # the loop -- to save reverse-reordering of them.

    bytesTrsvPerIter = 12.*(nnz(L)+nnz(U)) + 2.*4*size(L,1) + 2.*8*(size(L,1) + size(L,2))
    bytesSpMVPerIter = 12.*nnz(A) + 4*size(A,1) + 8*(size(A,1) + size(A,2))

    k = 1

    # Pre-allocation of space. Should get rid of once Linxiang's object removal works
    # TODO: remove it.
    Ap = Array(Cdouble, length(p))
    while k <= maxiter
        set_matrix_property(Dict(
            :r  => SA_HAS_DEDICATED_MEMORY,
            :x  => SA_HAS_DEDICATED_MEMORY,
            :p  => SA_HAS_DEDICATED_MEMORY,
          )
        )

        old_rz = rz

        spmv_time -= time()
        #Ap = A*p # Ap = Sparso.SpMV(A, p) # This takes most time. Compiler can reorder A to make faster
        #Sparso.SpMV!(Ap, A, p)
        # TODO: remove A_mul_B!. Use A*p once object removal works
        A_mul_B!(Ap, A, p)
        spmv_time += time()

        blas1_time -= time()
        alpha = old_rz / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        rel_err = norm(r)/normr0
        blas1_time += time()

        if rel_err < tol
            break
        end

        blas1_time -= time()
        z = copy(r)
        blas1_time += time()

        trsv_time -= time()
        Base.SparseMatrix.fwdTriSolve!(L, z)
          # Could have written as z=L\z if \ is specialized for triangular
        Base.SparseMatrix.bwdTriSolve!(U, z)
          # Could have wrriten as z=U\z if \ is specialized for triangular
        trsv_time += time()

        blas1_time -= time()
        rz = dot(r, z)
        beta = rz/old_rz
        p = z + beta * p
        blas1_time += time()

        k += 1
    end

    total_time = time() - total_time

    bytesTrsv = bytesTrsvPerIter*k
    bytesSpMV = bytesSpMVPerIter*(k + 1)

    println("total = $(total_time)s ($((bytesTrsv + bytesSpMV)/total_time/1e9) gbps) trsv_time = $(trsv_time)s ($(bytesTrsv/trsv_time/1e9) gbps) spmv_time = $(spmv_time)s ($(bytesSpMV/spmv_time/1e9) gbps) blas1_time = $blas1_time")

    return x, k, rel_err
end

if length(ARGS) == 0
    m = 10
    A = generate_symmetric_nonzerodiagonal_sparse_matrix(m)
else
    A = matrix_market_read(ARGS[1], true, true)
    m = size(A, 1)
end
x       = zeros(m)
b       = ones(m)
tol     = 1e-7
maxiter = 20000

if length(ARGS) == 2
  test = ARGS[2]
else
  test = "julia"
end

if test == "call-repl"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_REPLACE_CALLS)
elseif test == "context"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_CONTEXT, SA_REPLACE_CALLS)
elseif test == "reorder"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)
elseif test == "verbose"
  set_options(SA_ENABLE, SA_USE_SPMP, SA_VERBOSE, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)
end

A2 = copy(A)
println("compiler warmup (ignored): ")
if test == "julia"
  x, k, rel_err = pcg_symgs_ilu0(x, A2, b, tol, maxiter)
else
  @acc x, k, rel_err = pcg_symgs_ilu0(x, A2, b, tol, maxiter)
end

x = zeros(m)
println("\nRUN: ")
if test == "julia"
  x, k, rel_err = pcg_symgs_ilu0(x, A, b, tol, maxiter)
else
  @acc x, k, rel_err = pcg_symgs_ilu0(x, A, b, tol, maxiter)
end

println("k=", k, " rel_err=", rel_err)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
import DFControl.Utils: searchdir, strip_split, getfirst
using DFControl.Structures: ismagnetic
#does this really need the checking for corruption stuff?
function read_xsf_header(::Type{T}, filename::String) where T
    open(filename) do f
        nx, ny, nz = 0, 0, 0
        origin     = zero(Point3{T})
        supercell  = zero(Mat3{T})
        read_split_parse(typ) = parse.(typ, split(readline(f)))
        lines = 0
        while true
            lines += 1
            line = readline(f)
            if line == "BEGIN_DATAGRID_3D_UNKNOWN"
                nx, ny, nz = read_split_parse(Int)
                origin     = Point3{T}(read_split_parse(T))
                t = collect(hcat([read_split_parse(T) for i=1:3]...)')::Matrix{T}
                supercell  = Mat3{T}(t)
                lines += 1 + 1 + 3
                break
            end
        end
        return nx, ny, nz, origin, supercell, lines
    end
end

function read_points_from_xsf(::Type{T}, filename::String) where {T <: AbstractFloat}
    open(filename) do f
        while !eof(f)
            line = readline(f)
            if line == " DATAGRID_3D_DENSITY" ||
               occursin("DATAGRID_3D_UNKNOWN", line)

                nx, ny, nz = parse.(Int, split(readline(f)))
                origin     = Point3{T}(parse.(T, split(readline(f))))
                a_vec      = Vec3{T}(parse.(T, split(readline(f))))
                b_vec      = Vec3{T}(parse.(T, split(readline(f))))
                c_vec      = Vec3{T}(parse.(T, split(readline(f))))
                return [origin +
                        ia * a_vec +
                        ib * b_vec +
                        ic * c_vec for ia in range(0, 1, length=nx),
                                       ib in range(0, 1, length=ny),
                                       ic in range(0, 1, length=nz)]
            end
        end
    end
end
read_points_from_xsf(filename::String) = read_points_from_xsf(Float64, filename)

"""
read_values_from_xsf(filename::String, atom::Atom, T=Float64)

Returns an Array from reading a Wannier wavefunction file.
"""
function read_values_from_xsf(::Type{T}, filename::String) where {T <: AbstractFloat}
    open(filename) do f
        while !eof(f)
            line = readline(f)
            if line == "PRIMVEC"
                cell  = [Point3{T}.(parse.(T, split(readline(f)))) for i=1:3]
            end

            if line == " DATAGRID_3D_DENSITY" || occursin("DATAGRID_3D_UNKNOWN", line)
                nx, ny, nz = parse.(Int, split(readline(f)))
	            for i = 1:4
		            readline(f)
	            end
                out     = Array{T}(undef, nx, ny, nz)
                line    = readline(f)
				counter = 1
                while line != "END_DATAGRID_3D"
                    for t in parse.(T, split(line))
	                    out[counter] = t
	                    counter += 1
                    end
                    line = readline(f)
                end
                return out
            end
        end
    end
end
read_values_from_xsf(filename::String) = read_values_from_xsf(Float64, filename)

"""
write_xsf(filename::String, wfc::Wfc3D{T}) where T<:AbstractFloat

Writes the real part of the Wfc3D to a .xsf file that is readable by XCrysden or VESTA.
"""
function write_xsf(filename::String, wfc, structure; value_func=x -> norm(x))
    open(filename,"w") do f
        origin = wfc.points[1,1,1]
        write(f,"# Generated from PhD calculations\n")
        write(f, "CRYSTAL\n")
        c = ustrip.(structure.cell')
        write(f, "PRIMVEC\n")
        write(f, "$(c[1,1]) $(c[1,2]) $(c[1,3])\n")
        write(f, "$(c[2,1]) $(c[2,2]) $(c[2,3])\n")
        write(f, "$(c[3,1]) $(c[3,2]) $(c[3,3])\n")
        write(f, "CONVVEC\n")
        write(f, "$(c[1,1]) $(c[1,2]) $(c[1,3])\n")
        write(f, "$(c[2,1]) $(c[2,2]) $(c[2,3])\n")
        write(f, "$(c[3,1]) $(c[3,2]) $(c[3,3])\n")
        write(f, "PRIMCOORD\n")
        write(f, "$(length(structure.atoms)) 1\n")
        for at in structure.atoms
            n = at.element.symbol
            p = ustrip.(at.position_cart)
            write(f, "$n $(p[1]) $(p[2]) $(p[3])\n")
        end
        write.((f,),["", "BEGIN_BLOCK_DATAGRID_3D\n", "3D_FIELD\n", "BEGIN_DATAGRID_3D_UNKNOWN\n"])
        write(f,"$(size(wfc)[1])    $(size(wfc)[2])     $(size(wfc)[3])\n")
        write(f,"$(origin[1])   $(origin[2])   $(origin[3])\n")
        write(f,"$(wfc.points[end,1,1][1]-origin[1])   $(wfc.points[end,1,1][2]-origin[2])   $(wfc.points[end,1,1][3]-origin[3])\n")
        write(f,"$(wfc.points[1,end,1][1]-origin[1])   $(wfc.points[1,end,1][2]-origin[2])   $(wfc.points[1,end,1][3]-origin[3])\n")
        write(f,"$(wfc.points[1,1,end][1]-origin[1])   $(wfc.points[1,1,end][2]-origin[2])   $(wfc.points[1,1,end][3]-origin[3])\n")
        for wfp in wfc.values
            write(f,"$(value_func(wfp)) ")
        end
        write(f,"\n")
        write.((f,), ["END_DATAGRID_3D\n", "END_BLOCK_DATAGRID_3D\n"])
    end
end

@doc raw"""
	readhami(chk, eig_file::AbstractString)

Reads `eig_file` and uses the w90 checkpoint info in `chk` to return a vector of TbBlocks with the hopping parameters of the Wannier Tight Binding Hamiltonian.
"""
function readhami(chk, eig_file::AbstractString)
    v_mat = chk.V_matrix
    eigvals = read_eig(eig_file)
    R_cryst, degens = chk.ws_R_cryst, chk.ws_degens
    ws_shifts, ws_nshifts = chk.ws_shifts_cryst, chk.ws_nshifts
    Hq = map(1:length(chk.kpoints)) do ik
        v = v_mat[1:num_states(chk, ik), 1:chk.n_wann, ik]

        return v' * diagm(eigvals[disentanglement_range(chk, ik), ik]) * v
    end
    c = chk.cell'
    LT = eltype(c)
    T = eltype(eigvals)
    Hr_t = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for R in R_cryst]
    fourier_q_to_R(chk.kpoints, R_cryst) do iR, ik, phase
        @inbounds Hr_t[iR] .+= phase .* Hq[ik]
    end
    for o in Hr_t
        o ./= length(chk.kpoints)
    end

    out = [TbBlock(R, c*R, zeros(ComplexF64, chk.n_wann, chk.n_wann), zeros(ComplexF64, chk.n_wann, chk.n_wann)) for R in R_cryst]

    for (iR, (h, R, shifts, nshifts, d)) in enumerate(zip(Hr_t, R_cryst, ws_shifts, ws_nshifts, degens))
        for i in eachindex(h)
            ns = nshifts[i]
            frac = 1/(ns * d)
            for is = 1:ns
                rcryst = R + shifts[i][is]
                h1 = out[rcryst]
                if h1 === nothing
                    h1 = TbBlock{T, LT, Matrix{Complex{T}}}(rcryst, c*rcryst, zeros(ComplexF64, chk.n_wann, chk.n_wann), zeros(ComplexF64, chk.n_wann, chk.n_wann))
                    push!(out, h1)
                end
                h1.block[i] += h[i] * frac
                h1.tb_block[i] = h[i]
            end
        end
    end
    return out
end

#super not optimized
#TODO Test: new wigner seitz shift stuff
@doc raw"""
	read_colin_hami(up_chk, down_chk, up_eig_file::AbstractString, down_eig_file::AbstractString)

Returns the colinear TbHami representing the up-down blocks of the Wannier Tight Binding Hamiltonian.
"""
function read_colin_hami(up_chk, down_chk, up_eig_file::AbstractString, down_eig_file::AbstractString)
	uphami   = readhami(up_chk, up_eig_file)
	downhami = readhami(down_chk, down_eig_file)
	dim = blocksize(uphami)
	@assert dim == blocksize(downhami) "Specified files contain Hamiltonians with different dimensions of the Wannier basis."

	u1 = uphami[1]
	d1 = downhami[u1.R_cryst]

	first = TbBlock(u1.R_cryst,
	                u1.R_cart,
	                ColinMatrix(block(u1), block(d1)),
	                ColinMatrix(u1.tb_block, d1.tb_block))

	outhami  = [first]
	for u in uphami[2:end]
    	d = downhami[u.R_cryst]
    	if d !== nothing
    		push!(outhami, TbBlock(u.R_cryst, u.R_cart, ColinMatrix(block(u), block(d)), ColinMatrix(u.tb_block, d.tb_block)))
		end
	end
	return outhami
end

"""
    readhami(job::Job)

Goes through the job and will attempt to read the hamiltonian files.
If it finds a colinear calculation in the job it will read the up and down hamiltonians,
if the job was either nonmagnetic or noncolinear it will read only one hamiltonian file (there should be only one).
"""
function readhami(job::Job)

    @assert any(x -> x isa Calculation{Wannier90}, job.calculations) "No wannier90 calculations found in the job."

    seedname = getfirst(x -> x isa Calculation{Wannier90}, job.calculations).name
    jld_file = joinpath(job, "hami.jld2")
    if ispath(jld_file)
        return DFC.load(jld_file)["hami"]
    end
	eig_files = reverse(searchdir(job, ".eig"))
	chk_files = reverse(searchdir(job, ".chk"))
	@assert !isempty(eig_files) "No eig files ($(seedname).eig) found."
	@assert !isempty(chk_files) "No chk files ($(seedname).chk) found."
	if DFC.Structures.ismagnetic(job.structure)
    	if !DFC.Structures.iscolin(job.structure) || any(x -> DFC.Calculations.hasflag(x, :lspinorb) && x[:lspinorb], job.calculations)
    		return make_noncolin.(readhami(read_chk(chk_files[1]), joinpath(job, eig_files[1])))
    	else
    		return read_colin_hami(read_chk.(chk_files)..., eig_files...)
    	end
	else
		return readhami(read_chk(chk_files[1]), joinpath(job, eig_files[1]))
	end
end

"""
read_rmn_file(filename::String, structure::AbstractStructure)

Returns and array of tuples that define the dipoles between the Wannier functions in different unit cells.
"""
function read_rmn_file(filename::String, structure::DFC.Structures.Structure)
    open(filename) do  f
        out = RmnBlock{Float64}[]
        readline(f)
        n_wanfun = parse(Int, readline(f))
        readline(f)
        while !eof(f)
            l= split(readline(f))
            R_cryst = Vec3(parse.(Int, l[1:3]))
            block = getfirst(x -> x.R_cryst == R_cryst, out)

            if block == nothing
                block = RmnBlock{Float64}(cell(structure)' * R_cryst, R_cryst, Matrix{Point3{Float64}}(I, n_wanfun, n_wanfun))
                push!(out, block)
            end
            dipole = Point3(parse.(Float64, l[6:2:10]))
            block.block[parse.(Int, l[4:5])...] = dipole
        end
        return out
    end
end

#stuff that is not used currently
# """
# read_potential_file(filename::String, T=Float64)
#
# Reads a Quantum Espresso potential output file.
# """
# function read_potential_file(filename::String, T=Float64)
#     tmp_data=nothing
#     open(filename) do f
#         primCell = zeros(I, 3, 3)
#         while !eof(f)
#             line = readline(f)
#             if line == "PRIMVEC"
#                 for i=1:3
#                     primCell[i,:] = split(readline(f)) |> parse(T)
#                 end
#             end
#             if line == "DATAGRID_3D_UNKNOWN"
#                 nx, ny, nz = parse.(Int, split(readline(f)))
#                 for i = 1:4
#                     split(readline(f))
#                 end
#                 tmp_data=Array(T, nx, ny, nz)
#                 line = readline(f)
#                 k, k1, k2=1, 1, 1
#                 while line!= "END_DATAGRID_3D"
#                     tmp = Array{T}(map(x->(v = tryparse(T,x); isnull(v) ? 0.0 : get(v)),split(line)))
#                     for t in tmp
#                         tmp_data[k,k1,k2] = t
#                         if k<nx
#                             k+=1
#                         else
#                             k=1
#                             k1+=1
#                             if k1>ny
#                                 k1=1
#                                 k2+=1
#                                 if k2>nz
#                                     k2=1
#                                 end
#                             end
#                         end
#                     end
#                     line = readline(f)
#                 end
#                 return interpolate(tmp_data,BSpline(Cubic(Line())),OnGrid())
#             end
#         end
#     end
# end

#-------------------------Not currently used beyond here!------------------#
function write_dip_file(filename::String,points,cell,atoms,names,direction)
    open(filename,"w") do f
        origin = points[1,1,1][1]
        write(f,["# Generated from PhD calculations\n"])
        write(f,["CRYSTAL\n","PRIMVEC\n"])
        write(f,["$(cell[1][1])  $(cell[1][2]) $(cell[1][3])\n"])
        write(f,["$(cell[2][1])  $(cell[2][2]) $(cell[2][3])\n"])
        write(f,["$(cell[3][1])  $(cell[3][2]) $(cell[3][3])\n"])
        write(f,["CONVVEC\n"])
        write(f,["$(cell[1][1])  $(cell[1][2]) $(cell[1][3])\n"])
        write(f,["$(cell[2][1])  $(cell[2][2]) $(cell[2][3])\n"])
        write(f,["$(cell[3][1])  $(cell[3][2]) $(cell[3][3])\n"])
        write(f,["PRIMCOORD\n"])
        write(f,["  $(length(atoms)) 1\n"])
        for (atom,name) in zip(atoms,names)
            write(f,["$name  $(atom.center[1])  $(atom.center[2])  $(atom.center[3])\n"])
        end

        write(f,["", "BEGIN_BLOCK_DATAGRID_3D\n", "3D_FIELD\n",
        "BEGIN_DATAGRID_3D_UNKNOWN\n"])
        write(f,"$(size(points)[1])    $(size(points)[2])     $(size(points)[3])\n")
        write(f,"$(origin[1])   $(origin[2])   $(origin[3])\n")
        write(f,"$(points[end,1,1][1][1]-origin[1])   $(points[end,1,1][1][2]-origin[2])   $(points[end,1,1][1][3]-origin[3])\n")
        write(f,"$(points[1,end,1][1][1]-origin[1])   $(points[1,end,1][1][2]-origin[2])   $(points[1,end,1][1][3]-origin[3])\n")
        write(f,"$(points[1,1,end][1][1]-origin[1])   $(points[1,1,end][1][2]-origin[2])   $(points[1,1,end][1][3]-origin[3])\n")
        for point in points
            write(f,"$(point[2][direction]) ")
        end
        write(f,"\n")
        write(f,["END_DATAGRID_3D\n", "END_BLOCK_DATAGRID_3D\n"])
    end
end

function read_xsf_file_GPU(filename::String, T=Float64)
    open(filename) do f
        while !eof(f)
            line = readline(f)
            if line == " DATAGRID_3D_DENSITY"
                line0 = split(readline(f))
                line1 = split(readline(f))
                line2 = split(readline(f))
                line3 = split(readline(f))
                line4 = split(readline(f))
                nx = parse(Int,line0[1])
                ny = parse(Int,line0[2])
                nz = parse(Int,line0[3])
                out = Complex(T,nx,ny,nz)
                line = readline(f)
                k=1
                k1=1
                k2=1
                while line!= "END_DATAGRID_3D"
                    tmp = Array{Complex{T}}(map(x->(v = tryparse(T,x); isnull(v) ? Complex(0.0,0.0) : Complex{T}(get(v),0.0)),split(line)))
                    for t in tmp
                        out[k,k1,k2] = t
                        if k<nx
                            k+=1
                        else
                            k=1
                            k1+=1
                            if k1>ny
                                k1=1
                                k2+=1
                                if k2>nz
                                    k2=1
                                end
                            end
                        end
                    end
                    line = readline(f)
                end
                return out
            end
            if line == "BEGIN_DATAGRID_3D_UNKNOWN"
                line0 = split(readline(f))
                line1 = split(readline(f))
                line2 = split(readline(f))
                line3 = split(readline(f))
                line4 = split(readline(f))
                nx = parse(Int,line0[1])
                ny = parse(Int,line0[2])
                nz = parse(Int,line0[3])
                out = Array(Complex{T},nx,ny,nz)
                line = readline(f)
                k=1
                k1=1
                k2=1
                while line!= "END_DATAGRID_3D"
                    tmp = Array{Complex{T}}(map(x->(v = tryparse(T,x); isnull(v) ? Complex(0.0,0.0) : Complex{T}(get(v),0.0)),split(line)))
                    for t in tmp
                        out[k,k1,k2] = t
                        if k<nx
                            k+=1
                        else
                            k=1
                            k1+=1
                            if k1>ny
                                k1=1
                                k2+=1
                                if k2>nz
                                    k2=1
                                end
                            end
                        end
                    end
                    line = readline(f)
                end
                return out
            end
        end
    end
end

# function write_dipole_mesh(filename,mesh::Array{Tuple{Point3{T},Point3{T}},3},direction) where T
#     tmp_points = similar(mesh,WfPoint3{T})
#     for (ip,p) in enumerate(mesh)
#         tmp_points[ip] = WfcPoint3{T}(getfield(p[2],direction),p[1])
#     end
# end

# function write_exchanges(filename::String, structure::Structure)
#     open(filename, "w") do f
#         exchanges = data(structure)[:exchanges]
#         for (i, atom1) in enumerate(atoms(structure))
#             for (i2, atom2) in enumerate(atoms(structure)[i + 1:end])
#                 J = exchange_between(atom1, atom2, exchanges)
#                 if J != 0
#                     write(f, "$i:$(name(atom1)) -> $(i2+i):$(name(atom2)) $J\n")
#                 end
#             end
#         end
#     end
# end

function read_exchanges(filename::String, T=Float64)
    open(filename, "r") do f
        out = T[]
        indices = Tuple{Int,Int}[]
        max = 0
        while !eof(f)
            spl = split(readline(f))
            i, j = parse.(Int, spl[1:2])
            if j> max
                max = j
            end
            push!(indices, (i,j))
            exch = parse(T, spl[3])
            push!(out, exch)
        end
        real_out = zeros(T, (max, max))
        for (o, (i, j)) in zip(out, indices)
            real_out[i, j] = o
        end
        return real_out
    end
end

# function read_chk(filename)
#     # try
#         f = FortranFile(filename)
#         header = String(read(f, FString{33}))
#         n_bands = Int(read(f, Int32))
#         n_excluded_bands = Int(read(f, Int32))
#         # exclude_bands_t = zeros(Int32, n_excluded_bands)
#         read(f, (Int32, n_excluded_bands))
#         # exclude_bands= convert.(Int, exclude_bands_t)
#         read(f)
#         real_lattice = Mat3(read(f, (Float64, 3, 3))...)
#         recip_lattice = K_CART_TYPE{Float64}.(Mat3(read(f, (Float64, 3, 3))...)')
#         n_kpoints = read(f, Int32)
#         mp_grid = Vec3(Int.(read(f, (Int32, 3)))...)
#         kpt_lattice_t = read(f, (Float64, 3, n_kpoints))
#         kpt_lattice = [Vec3(kpt_lattice_t[:, i]...) for i = 1:size(kpt_lattice_t)[2]]
#         k_nearest_neighbors = read(f, Int32)

#         n_wann = read(f, Int32)
#         chkpt = strip(String(read(f, FString{20})))
#         have_disentangled = read(f, Int32) == 1 ? true : false
#         if have_disentangled
#             omega_invariant = read(f, Float64)
#             lwindow = map(x-> x==1 ? true : false, read(f, (Int32, n_bands, n_kpoints)))
#             ndimwin = read(f, (Int32, n_kpoints))
#             U_matrix_opt = read(f, (Complex{Float64}, n_bands, n_wann, n_kpoints))
#         else
#             omega_invariant = 0.0
#             lwindow = fill(true, 1, n_kpoints)
#             ndimwin = fill(n_wann, n_kpoints)
#             U_matrix_opt= Array{Complex{Float64}, 3}()
#         end
#         U_matrix = permutedims(read(f, (Complex{Float64}, n_wann, n_wann, n_kpoints)),(2,1,3))

#         # Combined effect of disentanglement and localization
#         V_matrix = Array{Complex{Float64},3}(undef, n_bands, n_wann, n_kpoints)
#         if have_disentangled
#             for ik in 1:n_kpoints
#                 V_matrix[:, :, ik] = U_matrix_opt[:, :, ik] * U_matrix[:, :, ik]
#             end
#         else
#             V_matrix = U_matrix
#         end

#         m_matrix = read(f, (Complex{Float64}, n_wann, n_wann, k_nearest_neighbors, n_kpoints))
#         wannier_centers_t = read(f, (Float64, 3, n_wann))
#         wannier_centers = [Point3(wannier_centers_t[:, i]...) for i = 1:size(wannier_centers_t)[2]]
#         wannier_spreads = read(f, (Float64, n_wann))

#         wb = read(f, (Float64, k_nearest_neighbors)) #this requires patched w90

#         return (
#             n_bands=n_bands,
#             n_excluded_bands=n_excluded_bands,
#             cell=real_lattice',
#             recip_cell=recip_lattice,
#             n_kpoints=n_kpoints,
#             mp_grid=mp_grid,
#             kpoints=kpt_lattice,
#             n_nearest_neighbors=k_nearest_neighbors,
#             neighbor_weights=wb,
#             n_wann=n_wann,
#             have_disentangled=have_disentangled,
#             Ω_invariant = omega_invariant,
#             lwindow=lwindow,
#             ndimwin =ndimwin,
#             U_matrix_opt=U_matrix_opt,
#             U_matrix=U_matrix,
#             V_matrix=V_matrix,
#             m_matrix=m_matrix,
#             wannier_centers=wannier_centers,
#             wannier_spreads=wannier_spreads
#         )
    # catch
    #     open(filename, "r") do f
    #         header = readline(f)
    #         n_bands = parse(Int, readline(f))
    #         n_excluded_bands = parse(Int,readline(f))
    #         # exclude_bands_t = zeros(Int32, n_excluded_bands)
    #         for i=1:n_excluded_bands
    #             readline(f)
    #         end
    #         # exclude_bands= convert.(Int, exclude_bands_t)
    #         # read(f)
    #         ssl = () -> strip_split(readline(f))
    #         real_lattice = Mat3(parse.(Float64, ssl()))
    #         recip_lattice = K_CART_TYPE{Float64}.(parse.(Float64, ssl())')
    #         n_kpoints = parse(Int, readline(f))
    #         mp_grid = parse(Vec3{Int}, ssl())
    #         kpt_lattice = [zero(Vec3{Float64}) for i=1:n_kpoints]
    #         for i = 1:n_kpoints
    #             kpt_lattice[i] = parse(Vec3{Float64}, ssl())
    #         end
    #         k_nearest_neighbors = parse(Int, readline(f))

    #         n_wann = parse(Int, readline(f))
    #         chkpt = readline(f)
    #         have_disentangled = parse(Bool, readline(f))
    #         if have_disentangled
    #             omega_invariant = parse(Float64, readline(f))
    #             lwindow = [parse(Bool, readline(f)) for i=1:n_bands, j=1:n_kpoints]
    #             # lwindow = map(x-> x == 1 ? true : false, )
    #             ndimwin = [parse(Int, readline(f)) for i=1:n_kpoints]
    #             U_matrix_opt = [complex(parse.(Float64, ssl())...) for i=1:n_bands, j=1:n_wann, k=1:n_kpoints]
    #         else
    #             omega_invariant = 0.0
    #             lwindow = fill(true, 1, n_kpoints)
    #             ndimwin = fill(n_wann, n_kpoints)
    #             U_matrix_opt= Array{Complex{Float64}, 3}()
    #         end
    #         U_matrix = permutedims([complex(parse.(Float64, ssl())...) for i=1:n_wann, j=1:n_wann, k=1:n_kpoints],(2,1,3))

    #         # Combined effect of disentanglement and localization
    #         V_matrix = Array{Complex{Float64},3}(undef, n_bands, n_wann, n_kpoints)
    #         if have_disentangled
    #             for ik in 1:n_kpoints
    #                 V_matrix[:, :, ik] = U_matrix_opt[:, :, ik] * U_matrix[:, :, ik]
    #             end
    #         else
    #             V_matrix = U_matrix
    #         end

    #         m_matrix = [complex(parse.(Float64, ssl())...) for i=1:n_wann, j=1:n_wann, k=1:k_nearest_neighbors, l=1:n_kpoints]
    #         wannier_centers = [parse(Point3{Float64}, ssl()) for i=1:n_wann]
    #         wannier_spreads = [parse(Float64, readline(f)) for i=1:n_wann]

    #         # wb = [parse(Float64, readline(f)) for i=1:k_nearest_neighbors] #this requires patched w90

    #         return (
    #             n_bands=n_bands,
    #             n_excluded_bands=n_excluded_bands,
    #             cell=real_lattice',
    #             recip_cell=recip_lattice,
    #             n_kpoints=n_kpoints,
    #             mp_grid=mp_grid,
    #             kpoints=kpt_lattice,
    #             n_nearest_neighbors=k_nearest_neighbors,
    #             neighbor_weights=wb,
    #             n_wann=n_wann,
    #             have_disentangled=have_disentangled,
    #             Ω_invariant = omega_invariant,
    #             lwindow=lwindow,
    #             ndimwin =ndimwin,
    #             U_matrix_opt=U_matrix_opt,
    #             U_matrix=U_matrix,
    #             V_matrix=V_matrix,
    #             m_matrix=m_matrix,
    #             wannier_centers=wannier_centers,
    #             wannier_spreads=wannier_spreads
    #         )
    #     end
    #     # return read_chk_formatted(filename)
    # end
# end
#TODO: cleanup
function read_chk(filename)
    f = FortranFile(filename)
    header = String(read(f, FString{33}))
    n_bands = Int(read(f, Int32))
    n_excluded_bands = Int(read(f, Int32))
    if n_excluded_bands != 0
        read(f, (Int32, n_excluded_bands))
    else
        read(f)
        read(f)
    end
    real_lattice = 1DFControl.angstrom.*(Mat3(read(f, (Float64, 3, 3))...))
    recip_lattice = K_CART_TYPE{Float64}.(Mat3(read(f, (Float64, 3, 3))...)')
    n_kpoints = read(f, Int32)
    mp_grid = Vec3(Int.(read(f, (Int32, 3)))...)
    kpt_lattice_t = read(f, (Float64, 3, n_kpoints))
    kpt_lattice = [Vec3(kpt_lattice_t[:, i]...) for i = 1:size(kpt_lattice_t)[2]]
    k_nearest_neighbors = read(f, Int32)

    n_wann = read(f, Int32)
    chkpt = strip(String(read(f, FString{20})))
    have_disentangled = read(f, Int32) != 0 ? true : false
    if have_disentangled
        omega_invariant = read(f, Float64)
        lwindow = map(x-> x!=0 ? true : false, read(f, (Int32, n_bands, n_kpoints)))
        ndimwin = read(f, (Int32, n_kpoints))
        U_matrix_opt = read(f, (Complex{Float64}, n_bands, n_wann, n_kpoints))
    else
        omega_invariant = 0.0
        lwindow = fill(true, n_bands, n_kpoints)
        ndimwin = fill(n_wann, n_kpoints)
        U_matrix_opt= [i == j ? 1.0im : 0.0im for i=1:n_bands, j=1:n_wann, ik = 1:n_kpoints]
    end
    U_matrix = read(f, (Complex{Float64}, n_wann, n_wann, n_kpoints))

    # Combined effect of disentanglement and localization
    V_matrix = Array{Complex{Float64},3}(undef, n_bands, n_wann, n_kpoints)
    if have_disentangled
        for ik in 1:n_kpoints
            V_matrix[:, :, ik] = U_matrix_opt[:, :, ik] * U_matrix[:, :, ik]
        end
    else
        V_matrix = U_matrix
    end

    m_matrix = read(f, (Complex{Float64}, n_wann, n_wann, k_nearest_neighbors, n_kpoints))
    wannier_centers_t = read(f, (Float64, 3, n_wann))
    wannier_centers = [Point3(wannier_centers_t[:, i]...) for i = 1:size(wannier_centers_t)[2]]
    wannier_spreads = read(f, (Float64, n_wann))
    wb = nothing
    try
        wb = read(f, (Float64, k_nearest_neighbors)) #this requires patched w90
    catch
        @warn "neighbor weights not found, berry calculations won't work. Patch your w90 if this functionality is wanted"
        wb = nothing
    end
    ws_R_cryst, ws_degens = wigner_seitz_points(mp_grid, metrics(ustrip.(real_lattice), ustrip.(recip_lattice)).real)
    ws_shifts_cryst, ws_nshifts = find_wigner_seitz_shifts(ws_R_cryst, wannier_centers, real_lattice, mp_grid)
    return (
        n_bands=n_bands,
        n_excluded_bands=n_excluded_bands,
        cell=real_lattice',
        recip_cell=recip_lattice,
        n_kpoints=n_kpoints,
        mp_grid=mp_grid,
        kpoints=kpt_lattice,
        n_nearest_neighbors=k_nearest_neighbors,
        neighbor_weights=wb,
        n_wann=n_wann,
        have_disentangled=have_disentangled,
        Ω_invariant = omega_invariant,
        lwindow=lwindow,
        ndimwin =ndimwin,
        U_matrix_opt=U_matrix_opt,
        U_matrix=U_matrix,
        V_matrix=V_matrix,
        m_matrix=m_matrix,
        wannier_centers=wannier_centers,
        wannier_spreads=wannier_spreads,
        ws_R_cryst = ws_R_cryst,
        ws_degens = ws_degens,
        ws_shifts_cryst = ws_shifts_cryst,
        ws_nshifts = ws_nshifts
    )
end
function read_chk(job::Job)
    if DFC.iscolin(job.structure)
        return map(s -> read_chk(joinpath(job, "$(name(getfirst(x -> eltype(x)==Wannier90&& x[:spin] == s, job.calculations))).chk")), ["up", "down"])
    else
        return read_chk(joinpath(job, "$(name(getfirst(x -> eltype(x)==Wannier90, job.calculations))).chk"))
    end
end

function read_eig(filename)
    t = readdlm(filename)
    n_bands = maximum(t[:,1])
    n_kpoints = maximum(t[:,2])

    #we follow order of w90 eigval matrix
    Hk = Matrix{Float64}(undef, Int(n_bands), Int(n_kpoints))
    for x in 1:size(t)[1]
        tv = t[x, :]
        Hk[Int(tv[1]), Int(tv[2])] = tv[3]
    end
    return Hk
end

metrics(chk) = metrics(ustrip.(chk.cell), ustrip.(chk.recip_cell))
function metrics(cell, recip_cell)
    real  = zeros(3, 3)
    recip = zeros(3, 3)
    for j in 1:3, i in 1:j
        for l in 1:3
            real[i, j]  += cell[i, l]  * cell[j, l]
            recip[i, j] += recip_cell[i, l] * recip_cell[j, l]
        end
        if i < j
            real[j, i]  = real[i, j]
            recip[j, i] = recip[j, i]
        end
    end
    return (real = real, recip = recip)
end

#This is a straight translation from the function in W90, this give the wigner_seitz R points
# The point of this is to determine the R_cryst but also the degeneracies i.e. the periodic images that have
# the exact same distance and will thus have exactly the same Tb hamiltonian block.
# This means that if one would be interpolating kpoings without dividing by the degeneracies, the periodic images
# would be "Double counted", which is why we divide by degen. In the actual tb hamiltonian this is fine though, no division needed.
wigner_seitz_points(chk) = wigner_seitz_points(chk.mp_grid, metrics(chk).real)
function wigner_seitz_points(mp_grid, real_metric)
    nrpts   = 0
    r_degens  = Int[]
    r = Vec3{Int}[]
    for n1 in -mp_grid[1] : mp_grid[1], n2 in -mp_grid[2] : mp_grid[2], n3 in -mp_grid[3] : mp_grid[3]
        R = Vec3(n1, n2, n3)
        dist_R0  = 0.0
        min_dist = typemax(Float64)
        ndegen   = 1
        best_R = copy(R)
        for i1 in -2:2, i2 in -2:2, i3 in -2:2
            ndiff = R .- Vec3(i1, i2, i3) .* mp_grid
            dist = ndiff' * real_metric * ndiff
            if abs(dist - min_dist) < 1e-7
                ndegen += 1
            elseif dist < min_dist
                min_dist = dist
                ndegen   = 1
            end
            if i1 == i2 == i3 == 0
                dist_R0 = dist
            end
        end
        # Only if R is actually the smallest distance it gets added to the R_cryst.
        if abs(min_dist - dist_R0) < 1e-7
            push!(r, R)
            push!(r_degens, ndegen)
        end
    end
    return r, r_degens
end

const WS_DISTANCE_TOL = 1e-5
find_wigner_seitz_shifts(chk) = find_wigner_seitz_shifts(chk.ws_R_cryst, chk.wannier_centers, chk.cell, chk.mp_grid)
function find_wigner_seitz_shifts(R_cryst, wannier_centers, cell, mp_grid)
    nwann = length(wannier_centers)
    ws_shifts_cryst = [[Vec3{Int}[zero(Vec3{Int})] for i=1:nwann, j=1:nwann] for iR = 1:length(R_cryst)]
    ws_nshifts      = [zeros(Int, nwann,nwann) for iR = 1:length(R_cryst)]
    c = ustrip.(cell')
    ic = inv(c)
    for (iR, R) in enumerate(R_cryst)
        r_cart = c * R
        for i in 1:nwann, j=1:nwann
            best_r_cart = -wannier_centers[i] + r_cart + wannier_centers[j]
            nr = norm(best_r_cart)

            r_cryst = ic * best_r_cart

            for l in -3:3, m in -3:3, n in -3:3
                lmn = Vec3(l, m, n)
                test_r_cryst = r_cryst + lmn .* mp_grid
                test_r_cart  = c * test_r_cryst
                if norm(test_r_cart) < nr
                    best_r_cart = test_r_cart
                    nr = norm(test_r_cart)
                    ws_shifts_cryst[iR][i, j][1] = lmn .* mp_grid
                end
            end

            if nr < WS_DISTANCE_TOL
                ws_nshifts[iR][i, j] = 1
                ws_shifts_cryst[iR][i, j][1] = Vec3(0, 0, 0)
            else
                best_r_cryst = ic * best_r_cart
                orig_shift = ws_shifts_cryst[iR][i,j][1]
                for l in -3:3, m in -3:3, n in -3:3
                    lmn = Vec3(l, m, n)
                    test_r_cryst = best_r_cryst + lmn .* mp_grid
                    test_r_cart  = c * test_r_cryst
                    if abs(norm(test_r_cart) - nr) < WS_DISTANCE_TOL
                        ws_nshifts[iR][i, j] += 1
                        if ws_nshifts[iR][i, j] == 1
                            ws_shifts_cryst[iR][i, j][ws_nshifts[iR][i, j]] = orig_shift + lmn .* mp_grid
                        else
                            push!(ws_shifts_cryst[iR][i, j], orig_shift + lmn .* mp_grid)
                        end
                    end
                end
            end
        end
    end
    return ws_shifts_cryst, ws_nshifts
end

#TODO: speedup, we don't need all of them if things were disentangled

mutable struct ReciprocalOverlaps{T}
    k_id           ::Int
    neighbor_ids   ::Vector{Int}
    b_vectors_cryst::Vector{Vec3{Int}} #connecting vector in crystalline coordinates
    overlaps       ::Vector{Matrix{Complex{T}}}
end

function ReciprocalOverlaps{T}(k_id::Int, n_nearest_neighbors::Int, nbands::Int) where {T}
    overlap_matrices = [Matrix{Complex{T}}(undef, nbands, nbands) for i=1:n_nearest_neighbors]
    ReciprocalOverlaps{T}(k_id, zeros(Int, n_nearest_neighbors), zeros(Vec3{Int}, n_nearest_neighbors), overlap_matrices)
end

function read_uHu(file)
    try
        uHu_file = FortranFile(file)
        read(uHu_file, FString{20})
        nbands, nkpoints, n_nearest_neighbors = read(uHu_file, (Int32, 3))
        out = [Matrix{ComplexF64}(undef, nbands, nbands) for i=1:nkpoints*n_nearest_neighbors^2]
        for i=1:nkpoints*n_nearest_neighbors^2
            out[i] = transpose(read(uHu_file, (ComplexF64, nbands, nbands)))
        end
        return out
    catch
        open(file, "r") do f
            readline(f)
            nbands, nkpoints, n_nearest_neighbors = parse.(Int, strip_split(readline(f)))
            out = [Matrix{ComplexF64}(undef, nbands, nbands) for i=1:nkpoints*n_nearest_neighbors^2]
            for i=1:nkpoints*n_nearest_neighbors^2
                for j = 1:nbands
                    for k = 1:nbands
                        out[i][j,k] = complex(parse.(Float64, strip_split(readline(f)))...)
                    end
                end
            end
            return out
        end
    end
end

function fill_overlaps!(grid::Vector{AbInitioKPoint{T}}, mmn_filename::AbstractString, uHu_filename::AbstractString, wannier_chk_params) where {T}
    num_wann = wannier_chk_params.n_wann
    uHu_file = FortranFile(uHu_filename)
    read(uHu_file, FString{20})
    read(uHu_file, (Int32, 3))

    open(mmn_filename, "r") do f
        readline(f) #header
        nbands, nkpoints, n_nearest_neighbors = parse.(Int, strip_split(readline(f)))
        #pre setup uHu
        for k in grid
            k.uHu = [Matrix{Complex{T}}(undef, num_wann, num_wann) for m=1:n_nearest_neighbors, n=1:n_nearest_neighbors]
        end
        neighbor_counter = 1
        for i in 1:nkpoints*n_nearest_neighbors
            sline = strip_split(readline(f))
            cur_neighbor = mod1(neighbor_counter, n_nearest_neighbors)

            ik, ik2 = parse.(Int, sline[1:2])

            overlap_ab_initio_gauge = Matrix{Complex{T}}(undef, nbands, nbands)
            for n in eachindex(overlap_ab_initio_gauge)
                overlap_ab_initio_gauge[n] = complex(parse.(T, strip_split(readline(f)))...)
            end

            vmat_ik = wannier_chk_params.V_matrix[:, :, ik]
            vmat_ik2 = wannier_chk_params.V_matrix[:, :, ik2]
            first_band_id_ik = findfirst(wannier_chk_params.lwindow[:, ik])
            first_band_id_ik2 = findfirst(wannier_chk_params.lwindow[:, ik2])

            num_states_ik = wannier_chk_params.ndimwin[ik]
            num_states_ik2 = wannier_chk_params.ndimwin[ik2]

            V1 = vmat_ik[1:num_states_ik, 1:num_wann]
            V2 = vmat_ik2[1:num_states_ik2, 1:num_wann]

            disentanglement_range_k1 = first_band_id_ik:first_band_id_ik+num_states_ik-1
            disentanglement_range_k2 = first_band_id_ik2:first_band_id_ik2+num_states_ik2-1
            S12 = overlap_ab_initio_gauge[disentanglement_range_k1, disentanglement_range_k2]

            kpoint = grid[ik]

            vr = (wannier_chk_params.recip_cell * parse(Vec3{Int}, sline[3:5]) + grid[ik2].k_cart) - kpoint.k_cart

            V1_T = V1'
            S12_V2 = S12 * V2
            kpoint.overlaps[cur_neighbor] =  V1_T * S12_V2
            k_eigvals_mat = diagm(kpoint.eigvals[disentanglement_range_k1])
            kpoint.hamis[cur_neighbor] = V1_T * k_eigvals_mat * S12_V2
            neighbor_counter += 1
            for nearest_neighbor2 in 1:n_nearest_neighbors
                ik3 = kpoint.neighbors[nearest_neighbor2].k_id2
                first_band_id_ik3 = findfirst(wannier_chk_params.lwindow[:, ik3])
                num_states_ik3    = wannier_chk_params.ndimwin[ik3]

                V3 = wannier_chk_params.V_matrix[1:num_states_ik3, 1:num_wann, ik3]

                uHu_k2_k3 = transpose(read(uHu_file, (ComplexF64, nbands, nbands)))
                disentanglement_range_k3 = first_band_id_ik3:first_band_id_ik3+num_states_ik3-1
                kpoint.uHu[nearest_neighbor2, cur_neighbor] =  V3' * uHu_k2_k3[disentanglement_range_k3, disentanglement_range_k2] * V2
            end
        end
        return grid
    end
end

function read_KBonds(nnkpfile::AbstractString)
    open(nnkpfile, "r") do f
        line = readline(f)
        while line != "begin recip_lattice"
            line = readline(f)
        end
        recip_cell = Matrix{K_CART_TYPE{Float64}}(undef, 3, 3)
        recip_cell[:, 1] = parse.(Float64, split(readline(f))) .* 1/1DFControl.angstrom
        recip_cell[:, 2] = parse.(Float64, split(readline(f))) .* 1/1DFControl.angstrom
        recip_cell[:, 3] = parse.(Float64, split(readline(f))) .* 1/1DFControl.angstrom

        while line != "begin kpoints"
            line = readline(f)
        end
        nkpoints = parse(Int, strip(readline(f)))
        kpoints = Vector{Vec3{K_CART_TYPE{Float64}}}(undef, nkpoints)
        counter = 1
        line = readline(f)
        while line != "end kpoints"
            kpoints[counter] = recip_cell * Vec3(parse.(Float64, split(line)))
            line = readline(f)
            counter += 1
        end
        out = [KBond[] for i = 1:nkpoints]
        while line != "begin nnkpts"
            line = readline(f)
        end
        n_nearest_neighbors = parse(Int, strip(readline(f)))
        counter = 1
        line = readline(f)
        while line != "end nnkpts"
            counter = mod1(counter, n_nearest_neighbors)
            sline = strip_split(line)
            ik, ik2 = parse.(Int, sline[1:2])

            vr = (recip_cell * parse(Vec3{Int}, sline[3:5]) + kpoints[ik2]) - kpoints[ik]
            push!(out[ik], KBond(ik, ik2, vr))
            counter += 1
            line = readline(f)
        end
        return out
    end
end

function fill_k_neighbors!(kpoints::Vector{AbInitioKPoint{T}}, file::AbstractString, recip_cell::Mat3) where {T}
    kbonds = read_KBonds(file)
    @assert length(kbonds) == length(kpoints) "Number kpoints in seedname.nnkp doesn't match with the number of kpoints in seedname.chk."
    for ik = 1:length(kbonds)
        nntot = length(kbonds[ik])
        kpoints[ik].overlaps = [Matrix{Complex{T}}(undef, 0, 0) for ib = 1:nntot]
        kpoints[ik].hamis = [Matrix{Complex{T}}(undef, 0, 0) for ib = 1:nntot]
        kpoints[ik].neighbors = kbonds[ik]
        kpoints[ik].uHu = Matrix{Matrix{Complex{T}}}(undef, nntot, nntot)
    end
    return kpoints
end

function read_wannier_functions(job)
    wancalc_ids = findall(x -> eltype(x) == Wannier90, job.calculations)
    wanfuncs = Vector{WannierFunction}[]
    for wancalc in job.calculations[wancalc_ids]
        t_funcs = WannierFunction[]
        # here we assume that if it's a spinor calculation, and all the upre, upim, downre, downim are there
        if DFC.hasflag(wancalc, :spinors) && wancalc[:spinors]
            upre_files   = filter(x -> occursin(name(wancalc), x), searchdir(job, "upre"))
            upim_files   = filter(x -> occursin(name(wancalc), x), searchdir(job, "upim"))
            downre_files = filter(x -> occursin(name(wancalc), x), searchdir(job, "downre"))
            downim_files = filter(x -> occursin(name(wancalc), x), searchdir(job, "downim"))
            any(isempty.((upre_files, upim_files, downre_files, downim_files))) && continue
            points = read_points_from_xsf(upre_files[1])

            for (ur, ui, dr, di) in zip(upre_files, upim_files, downre_files, downim_files)
                push!(t_funcs, WannierFunction(ur, ui, dr, di, points))
            end

        else
            xsf_files = filter(x -> occursin(".xsf", x), DFC.find_files(job, name(wancalc)))
            isempty(xsf_files) && continue

            points = read_points_from_xsf(xsf_files[1])
            for f in xsf_files
                push!(t_funcs, WannierFunction(f, points))
            end
        end
        push!(wanfuncs, t_funcs)
    end
    return wanfuncs
end

read_unk(file) = occursin("NC", file) ? read_unk_noncollinear(file) : read_unk_collinear(file)

function read_unk_collinear(file)
    f = FortranFile(file)
    ngx, ngy, ngz, nk, nbnd = read(f, (Int32, 5))
    Uk = zeros(ComplexF64, ngx, ngy, ngz, nbnd, 1)
    for i in 1:nbnd
        record = FortranFiles.Record(f)
        read!(record, view(Uk,:,:,:,i,1))
        close(record)
    end
    return Uk
end

function read_unk_noncollinear(file)
    f = FortranFile(file)
    ngx, ngy, ngz, nk, nbnd = read(f, (Int32, 5))
    Uk = zeros(ComplexF64, ngx, ngy, ngz, nbnd, 2)
    for i in 1:nbnd
        record = FortranFiles.Record(f)
        read!(record, view(Uk,:,:,:,i,1))
        close(record)
        record = FortranFiles.Record(f)
        read!(record, view(Uk,:,:,:,i,2))
        close(record)
    end
    return Uk
end

function read_spn(filename)
    f = FortranFile(filename)
    read(f)
    nbnd, nk = read(f, (Int32, 2))
    Sx, Sy, Sz = [zeros(ComplexF64, nbnd, nbnd) for k = 1:nk],
                 [zeros(ComplexF64, nbnd, nbnd) for k = 1:nk],
                 [zeros(ComplexF64, nbnd, nbnd) for k = 1:nk]
    for ik in 1:nk
        t = read(f, (ComplexF64, 3, div(nbnd * (nbnd + 1), 2)))
        counter = 1
        for ib1 in 1:nbnd, ib2 in 1:ib1
            Sx[ik][ib2, ib1] = t[1, counter]
            Sx[ik][ib1, ib2] = conj(t[1, counter])
            Sy[ik][ib2, ib1] = t[2, counter]
            Sy[ik][ib1, ib2] = conj(t[2, counter])
            Sz[ik][ib2, ib1] = t[3, counter]
            Sz[ik][ib1, ib2] = conj(t[3, counter])
            counter += 1
        end
    end
    return (Sx, Sy, Sz)
end

function S_R(chk, Sx, Sy, Sz)
    #q is in wannier gauge
    nk = length(chk.kpoints)
    Sx_q = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nk]
    Sy_q = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nk]
    Sz_q = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nk]
    vmat = chk.V_matrix
    nwann = chk.n_wann
    for ik in 1:nk
        v = vmat[1:num_states(chk, ik), 1:nwann, ik]
        disr = disentanglement_range(chk, ik)
        Sx_q[ik] = v' * Sx[ik][disr, disr] * v
        Sy_q[ik] = v' * Sy[ik][disr, disr] * v
        Sz_q[ik] = v' * Sz[ik][disr, disr] * v
    end
    R_cryst, degens = chk.ws_R_cryst, chk.ws_degens
    nR = length(R_cryst)
    Sx_R = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nR]
    Sy_R = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nR]
    Sz_R = [zeros(ComplexF64, chk.n_wann, chk.n_wann) for i=1:nR]

    fourier_q_to_R(chk.kpoints, R_cryst) do iR, ik, phase
        Sx_R[iR] .+= Sx_q[ik] .* phase
        Sy_R[iR] .+= Sy_q[ik] .* phase
        Sz_R[iR] .+= Sz_q[ik] .* phase
    end
    for iR in 1:nR
        Sx_R[iR] ./= nk
        Sy_R[iR] ./= nk
        Sz_R[iR] ./= nk
    end
    return Sx_R, Sy_R, Sz_R
end

function readspin(spn_file, chk_file)
    Sx_dft, Sy_dft, Sz_dft = read_spn(spn_file)
    return S_R(read_chk(chk_file), Sx_dft, Sy_dft, Sz_dft)
end

function readspin(job::Job)
	chk_files = reverse(searchdir(job, ".chk"))
	spn_files = reverse(searchdir(job, ".spn"))
	isempty(chk_files) && error("No .chk files found in job dir: $(job.local_dir)")
	isempty(spn_files) && error("No .spn files found in job dir: $(job.local_dir)")
	if length(chk_files) > 1
    	error("Not implemented for collinear spin-polarized calculations")
    end
    return readspin(spn_files[1], chk_files[1])
end

wan_hash(job::Job) = hash(read_chk(job))

function read_nnkp(nnkp_file) #not everything, just what I need for now
    open(nnkp_file, "r") do f
        blocks = read_wannier_blocks(f)
        nkpoints = parse(Int, blocks[:kpoints][1])
        nntot = parse(Int, blocks[:nnkpts][1])
        nnlist = Matrix{Int32}(undef, nkpoints, nntot)
        nncell = Array{Int32, 3}(undef, 3, nkpoints, nntot)
        counter = 2
        for ik = 1:nkpoints
            for nn = 1:nntot
                sl = strip_split(blocks[:nnkpts][counter])
                nnlist[ik, nn] = parse(Int, sl[2])
                for i = 1:3
                    nncell[i, ik, nn] = parse(Int, sl[i+2])
                end
                counter += 1

            end
        end
        return nnlist, nncell
    end
end

function Rmn(chk, nnkp_file)
    kbonds = read_KBonds(nnkp_file)
    m_matrix = chk.m_matrix
    R_cryst, degens = chk.ws_R_cryst, chk.ws_degens
    nR = length(R_cryst)
    nwann = chk.n_wann
    nntot = chk.n_nearest_neighbors
    wb = chk.neighbor_weights
    r_R = [zeros(Vec3{ComplexF64}, nwann, nwann) for i=1:nR]
    fourier_q_to_R(chk.kpoints, R_cryst) do iR, ik, phase
        for m = 1:nwann
            for n = 1:nwann
                for nn = 1:nntot
                    r_R[iR][n, m] += m == n ? - wb[nn] .* ustrip.(kbonds[ik][nn].vr) .* imag(log(m_matrix[n, m, nn, ik])) .* phase : 1im*wb[nn] .* ustrip.(kbonds[ik][nn].vr) .* m_matrix[n, m, nn, ik] .* phase
                end
            end
        end
    end
    nk = length(chk.kpoints)
    for iR in 1:nR
        r_R[iR] ./= nk
    end
    return r_R
end

function read_wannier_blocks(f)
    out = Dict{Symbol, Any}()
    while !eof(f)
        l = readline(f)
        if occursin("begin", l)
            s = Symbol(split(l)[2])
            lines = AbstractString[]
            l = readline(f)
            while !occursin("end", l)
                push!(lines, l)
                l = readline(f)
            end
            out[s] = lines
        end
    end
    return out
end


text:
using VLMinorityGameKit
using BSON
using DataFrames
using Plots

path_to_game_world = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/GW-T-M10-S5-NA1001.bson"
path_to_results_dir = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/results/SIM-T-M10-S5-NA1001-L2.bson"
d = BSON.load(path_to_game_world)
gameWorld = d[:gameworld]

println("Building the game world - finished. Starting the simulation ...")

# setup the number of time steps -
number_of_samples = 10
numberOfTimeSteps = 14
numberOfTraders = 1001
price_array = Array{Float64,2}(undef, (numberOfTimeSteps), number_of_samples)
simulation_result_array = Array{DataFrame,1}(undef, number_of_samples)
for sample_index = 1:number_of_samples
    result = execute(gameWorld, numberOfTimeSteps; liquidity=2.0 * numberOfTraders);
    simulation_result_array[sample_index] = result.data
    println("Completed sample = $(sample_index) ...")
end

# grab the first table, and get the number of rows -
(number_of_table_rows,_) = size(simulation_result_array[1])

# package -
simulation_dictionary = Dict()
simulation_dictionary["number_of_samples"] = number_of_samples
simulation_dictionary["number_of_timesteps"] = number_of_table_rows
simulation_dictionary["number_of_traders"] = numberOfTraders
simulation_dictionary["memory"] = 10
simulation_dictionary["strategies_per_agent"] = 5
simulation_dictionary["liquidity"] = 2.0 * numberOfTraders
simulation_dictionary["simulation_result_array"] = simulation_result_array

# dump -
bson(path_to_results_dir, data=simulation_dictionary)

# result = simulate(gameWorld, numberOfTimeSteps; liquidity=10.0 * numberOfTraders);
println("Simulation is finished .... enjoy.")