text:
@testset "inference with abstract global variable" begin
    let
        vmod = gen_virtual_module()
        res = @analyze_toplevel context = vmod virtualize = false begin
            s = "julia"
            sum(s)
        end

        @test is_concrete(vmod, :s)
        test_sum_over_string(res)
    end

    @testset "union assignment" begin
        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                global globalvar
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
        end

        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported
            end

            @test is_abstract(vmod, :globalvar)
            @test vmod.globalvar.t === Union{String,Symbol}
            @test length(res.inference_error_reports) === 1
            er = first(res.inference_error_reports)
            @test er isa NoMethodErrorReport
            @test isa(er.t, Vector) # should be true
        end

        # sequential
        let
            vmod = gen_virtual_module()
            res = @analyze_toplevel context = vmod virtualize = false begin
                if rand(Bool)
                    globalvar = "String"
                else
                    globalvar = :Symbol
                end

                foo(s::AbstractString) = length(s)
                foo(globalvar) # union-split no method matching error should be reported

                globalvar = 10
                foo(globalvar) # no method matching error should be reported
            end

            @test is_concrete(vmod, :globalvar)
            @test is_analyzed(vmod, :globalvar)
            @test isa_analyzed(vmod.globalvar, Int)
            @test length(res.inference_error_reports) === 2
            let er = first(res.inference_error_reports)
                @test er isa NoMethodErrorReport
                @test isa(er.t, Vector)
            end
            let er = last(res.inference_error_reports)
                @test er isa NoMethodErrorReport
                @test !isa(er.t, Vector)
            end
        end
    end
end

@testset "getfield with abstract global variable" begin
    # nested module access will be resolved as a direct call of `getfield`
    let
        res = @analyze_toplevel begin
            module foo

            const bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    # this should work even if the accessed variable is not constant
    let
        res = @analyze_toplevel begin
            module foo

            bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end
end

@testset "cache separation from native execution" begin
    # we shouldn't use global code cache for native execution,
    # since it has really not been analyzed by JET

    let
        # # check if the sysimg contains the cache for `getproperty(::Nothing, ::Symbol)`
        # local found = false
        # for mi in only(methods(getproperty, (Nothing,Symbol))).specializations
        #     isnothing(mi) && continue
        #     if mi.specTypes === Tuple{typeof(getproperty),Nothing,Symbol}
        #         found |= true
        #     end
        # end
        # @assert found

        # we can still get error report from that frame
        result = report_call((Nothing,)) do a
            a.field
        end
        @test length(get_reports(result)) === 1
    end
end

@testset "simple invalidation" begin
    let # simple invalidation
        m = Module()

        # analyze the first definition
        @eval m foo(a, b) = (sum(a), b)
        result = @report_call m.foo([1,2,3], "julia")
        @test isempty(get_reports(result))

        # renew the definition and invalidate it
        @eval m foo(a, b) = (a, sum(b))
        result = @report_call m.foo([1,2,3], "julia")
        test_sum_over_string(result)
    end

    let # backedge invalidation
        m = Module()
        @eval m callf(f, args...) = f(args...)

        # analyze the first definition
        @eval m foo(a, b) = (sum(a), b)
        result = @report_call m.callf(m.foo, [1,2,3], "julia")
        @test isempty(get_reports(result))

        # renew the definition and invalidate it
        @eval m foo(a, b) = (a, sum(b))
        result = @report_call m.callf(m.foo, [1,2,3], "julia")
        test_sum_over_string(result)
    end
end

# COMBAK this test is very fragile, think about the alternate tests
# @testset "end to end invalidation" begin
#     # invalidation from deeper call site should still refresh JET analysis
#     let
#         # NOTE: branching on https://github.com/JuliaLang/julia/pull/38830
#         symarg = last(first(methods(Base.show_sym)).sig.parameters) === Symbol ?
#                  :(sym::Symbol) :
#                  :(sym)
#
#         l1, l2, l3 = @freshexec begin
#             # ensure we start with this "errorneous" `show_sym`
#             @eval Base begin
#                 function show_sym(io::IO, $(symarg); allow_macroname=false)
#                     if is_valid_identifier(sym)
#                         print(io, sym)
#                     elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
#                         print(io, '@')
#                         show_sym(io, sym_str[2:end]) # NOTE: `sym_str[2:end]` here is errorneous
#                     else
#                         print(io, "var", repr(string(sym)))
#                     end
#                 end
#             end
#
#             # should have error reported
#             result1 = @report_call println(QuoteNode(nothing))
#
#             # should invoke invalidation in the deeper call site of `println(::QuoteNode)`
#             @eval Base begin
#                 function show_sym(io::IO, $(symarg); allow_macroname=false)
#                     if is_valid_identifier(sym)
#                         print(io, sym)
#                     elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
#                         print(io, '@')
#                         show_sym(io, Symbol(sym_str[2:end]))
#                     else
#                         print(io, "var", repr(string(sym)))
#                     end
#                 end
#             end
#
#             # now we shouldn't have reports
#             result2 = @report_call println(QuoteNode(nothing))
#
#             # again, invoke invalidation
#             @eval Base begin
#                 function show_sym(io::IO, $(symarg); allow_macroname=false)
#                     if is_valid_identifier(sym)
#                         print(io, sym)
#                     elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
#                         print(io, '@')
#                         show_sym(io, sym_str[2:end])
#                     else
#                         print(io, "var", repr(string(sym)))
#                     end
#                 end
#             end
#
#             # now we should have reports, again
#             result3 = @report_call println(QuoteNode(nothing))
#
#             (length ∘ JET.get_reports).((result1, result2, result3)) # return
#         end
#
#         @test l1 > 0
#         @test l2 == 0
#         @test l3 == l1
#     end
# end

@testset "integration with global code cache" begin
    # analysis for `sum(::String)` is already cached, `sum′` and `sum′′` should use it
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            sum′(s) = sum(s)
            sum′′(s) = sum′(s)
            $report_call() do
                sum′′("julia")
            end
        end)
        test_sum_over_string(get_reports(result))
    end

    # incremental setup
    let
        m = gen_virtual_module()

        result = Core.eval(m, quote
            $report_call() do
                sum("julia")
            end
        end)
        test_sum_over_string(get_reports(result))

        result = Core.eval(m, quote
            sum′(s) = sum(s)
            $report_call() do
                sum′("julia")
            end
        end)
        test_sum_over_string(get_reports(result))

        result = Core.eval(m, quote
            sum′′(s) = sum′(s)
            $report_call() do
                sum′′("julia")
            end
        end)
        test_sum_over_string(get_reports(result))
    end

    # should not error for virtual stacktrace traversing with a frame for inner constructor
    # https://github.com/aviatesk/JET.jl/pull/69
    let # FIXME https://github.com/JuliaLang/julia/pull/41885
        res = @analyze_toplevel begin
            struct Foo end
            println(Foo())
        end
        @test isempty(res.inference_error_reports)
    end
end

@testset "integration with local code cache" begin
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            struct Foo{T}
                bar::T
            end
            $report_call((Foo{Int},)) do foo
                foo.baz # typo
            end
        end)

        @test !isempty(get_reports(result))
        @test !isempty(get_cache(result.analyzer))
        @test any(get_cache(result.analyzer)) do analysis_result
            analysis_result.argtypes==Any[CC.Const(getproperty),m.Foo{Int},CC.Const(:baz)]
        end
    end

    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            struct Foo{T}
                bar::T
            end
            getter(foo, prop) = getproperty(foo, prop)
            $report_call((Foo{Int}, Bool)) do foo, cond
                getter(foo, :bar)
                cond ? getter(foo, :baz) : getter(foo, :qux) # non-deterministic typos
            end
        end)

        # there should be local cache for each errorneous constant analysis
        @test !isempty(get_reports(result))
        @test !isempty(get_cache(result.analyzer))
        @test any(get_cache(result.analyzer)) do analysis_result
            analysis_result.argtypes==Any[CC.Const(m.getter),m.Foo{Int},CC.Const(:baz)]
        end
        @test any(get_cache(result.analyzer)) do analysis_result
            analysis_result.argtypes==Any[CC.Const(m.getter),m.Foo{Int},CC.Const(:qux)]
        end
    end
end

@testset "constant analysis" begin
    # constant prop should limit false positive union-split no method reports
    let
        m = @fixturedef begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        result = Core.eval(m, :($report_call(foo, (P, Int))))
        @test isempty(get_reports(result))

        # works for cache
        result = Core.eval(m, :($report_call(foo, (P, Int))))
        @test isempty(get_reports(result))
    end

    # more cache test, constant prop should re-run in deeper level
    let
        m = @fixturedef begin
            mutable struct P
                i::Int
                s::String
            end
            foo(p, i) = p.i = i
            bar(args...) = foo(args...)
        end

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away
        result = Core.eval(m, :($report_call(bar, (P, Int))))
        @test isempty(get_reports(result))

        # works for cache
        result = Core.eval(m, :($report_call(bar, (P, Int))))
        @test isempty(get_reports(result))
    end

    # constant prop should not exclude those are not related
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $report_call(foo, (P, Int, #= invalid =# Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::Int64)" should be threw away, while
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int64)" should be kept
        @test length(get_reports(result)) === 1
        er = first(get_reports(result))
        @test er isa NoMethodErrorReport
        @test er.t === Tuple{typeof(convert), Type{String}, Int}
    end

    # constant prop should narrow down union-split no method error to single no method matching error
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            mutable struct P
                i::Int
                s::String
            end
            function foo(p, i, s)
                p.i = i
                p.s = s
            end

            $report_call(foo, (P, String, Int))
        end)

        # "for one of the union split cases, no matching method found for signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Union{Type{Int64}, Type{String}}, v::String)" should be narrowed down to "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{Int}, v::String)"
        @test !isempty(get_reports(result))
        @test any(get_reports(result)) do report
            report isa NoMethodErrorReport &&
            report.t === Tuple{typeof(convert), Type{Int}, String}
        end
        # "no matching method found for call signature: Base.convert(Base.fieldtype(Base.typeof(x::P)::Type{P}, f::Symbol)::Type{String}, v::Int)"
        # won't be reported since `typeinf` early escapes on `Bottom`-annotated statement
    end

    # report-throw away with constant analysis shouldn't throw away reports from the same
    # frame but with the other constants
    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            bar() = foo(-1), foo(1) # constant analysis on `foo(1)` shouldn't throw away reports from `foo(-1)`
            $report_call(bar)
        end)
        @test !isempty(get_reports(result))
        @test any(r->isa(r,NoMethodErrorReport), get_reports(result))
    end

    let
        m = gen_virtual_module()
        result = Core.eval(m, quote
            foo(a) = a<0 ? a+string(a) : a
            function bar(b)
                a = b ? foo(-1) : foo(1)
                b = foo(-1)
                return a, b
            end
            $report_call(bar, (Bool,))
        end)
        @test !isempty(get_reports(result))
        # FIXME our report uniquify logic might be wrong and it wrongly singlifies the different reports here
        @test_broken count(isa(report, NoMethodErrorReport) for report in get_reports(result)) == 2
    end

    @testset "constant analysis throws away false positive reports" begin
        let
            m = @fixturedef begin
                foo(a) = a > 0 ? a : "minus"
                bar(a) = foo(a) + 1
            end

            # constant propagation can reveal the error pass can't happen
            result = Core.eval(m, :($report_call(()->bar(10))))
            @test isempty(get_reports(result))

            # for this case, no constant prop' doesn't happen, we can't throw away error pass
            result = Core.eval(m, :($report_call(bar, (Int,))))
            @test length(get_reports(result)) === 1
            er = first(get_reports(result))
            @test er isa NoMethodErrorReport
            @test er.t == [Tuple{typeof(+),String,Int}]

            # if we run constant prop' that leads to the error pass, we should get the reports
            result = Core.eval(m, :($report_call(()->bar(0))))
            @test length(get_reports(result)) === 1
            er = first(get_reports(result))
            @test er isa NoMethodErrorReport
            @test er.t === Tuple{typeof(+),String,Int}
        end

        # we should throw-away reports collected from frames that are revealed as "unreachable"
        # by constant prop'
        let
            m = @fixturedef begin
                foo(a) = bar(a)
                function bar(a)
                    return if a < 1
                        baz1(a, "0")
                    else
                        baz2(a, a)
                    end
                end
                baz1(a, b) = a ? b : b
                baz2(a, b) = a + b
            end

            # no constant prop, just report everything
            result = Core.eval(m, :($report_call(foo, (Int,))))
            @test length(get_reports(result)) === 1
            er = first(get_reports(result))
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # constant prop should throw away the non-boolean condition report from `baz1`
            result = Core.eval(m, quote
                $report_call() do
                    foo(1)
                end
            end)
            @test isempty(get_reports(result))

            # constant prop'ed, still we want to have the non-boolean condition report from `baz1`
            result = Core.eval(m, quote
                $report_call() do
                    foo(0)
                end
            end)
            @test length(get_reports(result)) === 1
            er = first(get_reports(result))
            @test er isa NonBooleanCondErrorReport &&
                er.t === Int

            # so `Bool` is good for `foo` after all
            result = Core.eval(m, :($report_call(foo, (Bool,))))
            @test isempty(get_reports(result))
        end

        # end to end
        let
            res = @analyze_toplevel begin
                function foo(n)
                    if n < 10
                        return n
                    else
                        return "over 10"
                    end
                end

                function bar(n)
                    if n < 10
                        return foo(n) + 1
                    else
                        return foo(n) * "+1"
                    end
                end

                bar(1)
                bar(10)
            end

            @test isempty(res.inference_error_reports)
        end
    end
end

@testset "additional analysis pass for task parallelism code" begin
    # general case with `schedule(::Task)` pattern
    result = report_call() do
        t = Task() do
            sum("julia")
        end
        schedule(t)
        fetch(t)
    end
    test_sum_over_string(get_reports(result))

    # handle `Threads.@spawn` (https://github.com/aviatesk/JET.jl/issues/114)
    result = report_call() do
        fetch(Threads.@spawn 1 + "foo")
    end
    @test length(get_reports(result)) == 1
    let r = first(get_reports(result))
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), Int, String}
    end

    # handle `Threads.@threads`
    result = report_call((Int,)) do n
        a = String[]
        Threads.@threads for i in 1:n
            push!(a, i)
        end
        return a
    end
    @test !isempty(get_reports(result))
    @test any(get_reports(result)) do r
        isa(r, NoMethodErrorReport) &&
        r.t === Tuple{typeof(convert), Type{String}, Int}
    end

    # multiple tasks in the same frame
    result = report_call() do
        t1 = Threads.@spawn 1 + "foo"
        t2 = Threads.@spawn "foo" + 1
        fetch(t1), fetch(t2)
    end
    @test length(get_reports(result)) == 2
    let r = get_reports(result)[1]
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), Int, String}
    end
    let r = get_reports(result)[2]
        @test isa(r, NoMethodErrorReport)
        @test r.t === Tuple{typeof(+), String, Int}
    end

    # nested tasks
    result = report_call() do
        t0 = Task() do
            t = Threads.@spawn sum("julia")
            fetch(t)
        end
        schedule(t0)
        fetch(t0)
    end
    test_sum_over_string(get_reports(result))

    # when `schedule` call is separated from `Task` definition
    make_task(s) = Task() do
        sum(s)
    end
    function run_task(t)
        schedule(t)
        fetch(t)
    end
    result = report_call() do
        t = make_task("julia")

        run_task(t)
    end
    test_sum_over_string(get_reports(result))
    let r = first(get_reports(result))
        # we want report to come from `run_task`, but currently we invoke JET analysis on `Task` construction
        @test_broken any(r.vst) do vf
            vf.linfo.def.name === :run_task
        end
    end

    # report uncaught exception happened in a task
    # TODO currently uncaught exceptions are erased by return type check at caller `Task(::Function)`
    result = report_call() do
        fetch(Threads.@spawn throw("foo"))
    end
    @test_broken length(get_reports(result)) == 1
    @test_broken isa(first(get_reports(result)), UncaughtExceptionReport)

    # don't fail into infinite loop (rather, don't spoil inference termination)
    m = @fixturedef begin
        # adapated from https://julialang.org/blog/2019/07/multithreading/
        import Base.Threads.@spawn

        # sort the elements of `v` in place, from indices `lo` to `hi` inclusive
        function psort!(v, lo::Int=1, hi::Int=length(v))
            if lo >= hi                       # 1 or 0 elements; nothing to do
                return v
            end
            if hi - lo < 100000               # below some cutoff, run in serial
                sort!(view(v, lo:hi), alg = MergeSort)
                return v
            end

            mid = (lo+hi)>>>1                 # find the midpoint

            half = @spawn psort!(v, lo, mid)  # task to sort the lower half; will run
            psort!(v, mid+1, hi)              # in parallel with the current call sorting
                                              # the upper half
            wait(half)                        # wait for the lower half to finish

            temp = v[lo:mid]                  # workspace for merging

            i, k, j = 1, lo, mid+1            # merge the two sorted sub-arrays
            @inbounds while k < j <= hi
                if v[j] < temp[i]
                    v[k] = v[j]
                    j += 1
                else
                    v[k] = temp[i]
                    i += 1
                end
                k += 1
            end
            @inbounds while k < j
                v[k] = temp[i]
                k += 1
                i += 1
            end

            return v
        end
    end
    result = report_call(m.psort!, (Vector{Int},))
    @test true
end

@testset "opaque closure" begin
    # can cache const prop' result with varargs
    function oc_varargs_constprop()
        oc = Base.Experimental.@opaque (args...)->args[1]+args[2]+arg[3] # typo on `arg[3]`
        return Val{oc(1,2,3)}()
    end
    result = @report_call oc_varargs_constprop()
    @test !isempty(get_cache(result.analyzer))
end

@testset "https://github.com/aviatesk/JET.jl/issues/133" begin
    res = @analyze_toplevel begin
        @ccall strlen("foo"::Cstring)::Csize_t
    end
    @test isempty(res.inference_error_reports)
end


text:
#!/usr/bin/env julia

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

using Pkg
Pkg.activate(normpath(@__DIR__))

using Colors, Images, ImageIO, ArgParse, Polynomials, Documenter
using ColorTypes:RGB
import FileIO: @format_str, query
using Raytracing

FILE_NAME = split(PROGRAM_FILE, "/")[end]

main(x::Union{String, Float64, Int64}...) = main([string(var) for var in [x...]])
function main(args)
	parsed_arguments = ArgParse_command_line(args) # the result is a Dict{String,Any}
	(isnothing(parsed_arguments)) && (return nothing)
	#print_ArgParseSettings(parsed_arguments)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

	if parsed_command=="demo"
		#println(parse_demo_settings(parsed_settings))
		demo(parse_demo_settings(parsed_settings)...)
	elseif parsed_command=="tonemapping"
		#println(parse_tonemapping_settings(parsed_settings))
		tone_mapping(parse_tonemapping_settings(parsed_settings)...)
	elseif parsed_command=="demo_animation"
		#println(parse_tonemapping_settings(parsed_settings))
		demo_animation(parse_demoanimation_settings(parsed_settings)...)
	elseif parsed_command=="render"
		#println(parse_render_settings(parsed_settings))
		render(parse_render_settings(parsed_settings)...)
	elseif parsed_command=="animation"
		#println(parse_render_settings(parsed_settings))
		render_animation(parse_render_animation_settings(parsed_settings)...)
	else
		throw(ArgumentError("unknown command $(parsed_command)"))
	end

	return nothing
end

if (ARGS==String[])
	println("\nwithout input arguments/commands, show this help message and exit\n")
	main(["--help"])
else
	main(ARGS)
end

ext: jl


text:
using Parameters
using Random
using Logging

import LinearAlgebra: mul!, ldiv!, transpose!
import Random: randn!

using ..UnitCells: UnitCell
using ..Lattices: Lattice, sorted_neighbor_table_perm!, loc_to_site, calc_neighbor_table
using ..Checkerboard: checkerboard_order, checkerboard_groups, checkerboard_mul!, checkerboard_transpose_mul!
using ..IterativeSolvers: GMRES, ConjugateGradient, BiCGStab
using ..Utilities: get_index, get_τ, δ, reshaped

export write_K_matrix!

struct SSHBond{T1<:AbstractFloat,T2<:Continuous} <: AbstractBond

    "Average hopping energy."
    t::T2

    "Standard deviation hopping energy."
    σt::T2

    "Average phonon frequency"
    ω::T1

    "Standard deviation phonon frequency."
    σω::T1

    "Anharmonic term coefficient for ω₄⋅X⁴ term."
    ω₄::T1

    "Standard deviation in ω₄"
    σω₄::T1

    "Average linear electron-phonon coupling"
    α::T1

    "Standard deviation linear electron-phonon coupling"
    σα::T1

    "Non-linear electron phonon coupling of the form α₂⋅X²."
    α₂::T1

    "Standard deviation in α₂."
    σα₂::T1

    "orbital"
    o₁::Int

    "orbital"
    o₂::Int

    "displacement in unit cells"
    v::Vector{Int}

    "whether or not a phonon lives on the bond."
    has_phonon::Bool

    "name of the phonon."
    name::String

    function SSHBond(t::T2,σt::T2,ω::T1,σω::T1,ω₄::T1,σω₄::T1,α::T1,σα::T1,α₂::T1,σα₂::T1,
                     o₁::Int,o₂::Int,v::AbstractVector{Int},name::String="") where {T1<:AbstractFloat,T2<:Continuous}

        @assert length(v)==3

        if name==""
            name = randstring(5)
        end

        v′ = zeros(Int,3)
        copyto!(v′,v)
        has_phonon = (!iszero(ω))||(!iszero(σω))
        return new{T1,T2}(t,σt,ω,σω,ω₄,σω₄,α,σα,α₂,σα₂,o₁,o₂,v′,has_phonon,name)
    end
end

mutable struct SSHModel{T1,T2,T3,T4} <: AbstractModel{T1,T2,T3,T4}

    "Inverse temperature"
    β::T1

    "Discretization in imaginary time direction"
    Δτ::T1

    "Length of imaginary time axis"
    Lτ::Int

    "Number of sites in physical lattice"
    Nsites::Int

    "Number of bonds/hoppings in lattice"
    Nbonds::Int

    "Number of phonons in lattice"
    Nph::Int

    "Dimension of M matrix"
    Ndim::Int

    "Number of degrees of freedom"
    Ndof::Int

    "Number of types of bonds in lattice."
    nbonds::Int

    "Number of types of phonons in lattice."
    nph::Int

    "Represents lattice"
    lattice::Lattice{T1}

    "Bond defintions, including description of phonon living on the bond if there is one."
    bond_definitions::Vector{SSHBond{T1,T2}}

    "Ndof=Nph*Lτ phonon fields. Ordered in vector by phonon defintion."
    x::Vector{T1}

    "Nbond bare electron hoppings."
    t::Vector{T2}

    "Nph phonon frequencies. Ordered in vector by phonon defintion."
    ω::Vector{T1}

    "Anharmonic phonon term ω₄⋅X⁴."
    ω₄::Vector{T1}

    "Nph electron-phonon couplings. Ordered in vector by phonon defintion."
    α::Vector{T2}

    "Non-linear electron-phonon coupling."
    α₂::Vector{T2}

    "Chemical potential"
    μ::Vector{T1}

    "Diagonal matrix exp{Δτ⋅μ}"
    expΔτμ::Vector{T1}

    """
    Modulated hoppings t′ = t-α⋅X.
    """
    t′::Matrix{T2}

    "Map field to phonon"
    field_to_phonon::Vector{Int}

    "Map Ndof phonon fields onto a time slice τ"
    field_to_τ::Vector{Int}

    "Maps every field to its equivalent primary field.
    If a field maps onto itself then it is primary."
    primary_field::Vector{Int}

    "Map phonon to bond"
    phonon_to_bond::Vector{Int}

    "Map bond to phonon. If 0 then no phonon to map to."
    bond_to_phonon::Vector{Int}

    "Maps a given phonon onto an SSHPhonon defintion."
    bond_to_definition::Vector{Int}

    """
    Checkerboard bond order. Maps bond onto neighbor table entry.
    """
    checkerboard_perm::Vector{Int}

    """
    Inverse checkerboard order. Map neighbor table entry onto bond.
    """
    inv_checkerboard_perm::Vector{Int}

    "Neighbor table telling which sites in lattice are connect by bonds.
    Ordered according to checkerboard decomposition."
    neighbor_table::Matrix{Int}

    "Matrix of dimensions (Lτ,Nbonds) containing cosh(t′) where t′=Δτ(t-α⋅x).
    Ordered in vector according to checkerboard decomposition."
    cosht::Matrix{T2}

    "Matrix of dimensions (Lτ,Nbonds) containing sinh(t′) where t′=Δτ(t-α⋅x).
    Ordered in vector according to checkerboard decomposition."
    sinht::Matrix{T2}

    "If true the default mul! routine multiplies by the M (or Mᵀ) matrix.
    If false the default mul! routine multiplies by the symmetric matrix MᵀM (or MMᵀ) instead."
    mul_by_M::Bool

    "If true multiply by Mᵀ (or MMᵀ) instead of M (or MᵀM)."
    transposed::Bool

    "A vector of length `Ndim` to temporarily store data.
    Used in multiplication by MᵀM or MMᵀ."
    v′::Vector{T2}

    "A vector of length `Ndim` to temporarily store data.
    Used to store Mᵀ⋅b result when solving M⋅x=b via MᵀM⋅x=Mᵀ⋅b."
    v″::Vector{T2}

    "A vector of length `Ndim` to temporarily store data.
    Used for calculating true residual error after linear solve."
    v‴::Vector{T2}

    """
    Iterative Solver
    """
    solver::T3

    """
    Random Number Generator
    """
    rng::T4

    """
    data folder
    """
    datafolder::String

    function SSHModel(lattice::Lattice{T1}, β::T1, Δτ::T1;
                     is_complex::Bool=false, iterativesolver::String="cg",rng::AbstractRNG,
                     tol::T1=1e-4, maxiter::Int=10000, restart::Int=-1) where {T1<:AbstractFloat}

        # determines data type of matrix elements of M
        if is_complex
            T2 = Complex{T1}
        else
            T2 = T1
        end

        # length of imaginary time axis
        Lτ = round(Int,β/Δτ)

        # number of sites in lattice
        Nsites = lattice.nsites

        # dimension of M matrix
        Ndim = Nsites*Lτ

        # zero phonons defined initially
        Nph = 0

        # zero degrees of freedom initially
        Ndof = 0

        # zero bonds in lattice initially
        Nbonds = 0

        # zero types of bonds initially
        nbonds = 0

        # zero types of phonons initially
        nph = 0

        # intialize vector to contain phonon defintions
        bond_definitions = Vector{SSHBond{T1,T2}}(undef,0)

        # initialize arrays
        x                     = Vector{T1}(undef,0)
        t                     = Vector{T2}(undef,0)
        ω                     = Vector{T1}(undef,0)
        ω₄                    = Vector{T1}(undef,0)
        α                     = Vector{T2}(undef,0)
        α₂                    = Vector{T2}(undef,0)
        μ                     = zeros(T1,Nsites)
        expΔτμ                = ones(T1,Nsites)
        t′                    = Matrix{T2}(undef,Lτ,0)
        checkerboard_perm     = Vector{Int}(undef,0)
        inv_checkerboard_perm = Vector{Int}(undef,0)
        neighbor_table        = Matrix{Int}(undef,2,0)
        cosht                 = Matrix{T2}(undef,Lτ,0)
        sinht                 = Matrix{T2}(undef,Lτ,0)
        field_to_phonon       = Vector{Int}(undef,0)
        field_to_τ            = Vector{Int}(undef,0)
        primary_field         = Vector{Int}(undef,0)
        phonon_to_bond        = Vector{Int}(undef,0)
        bond_to_phonon        = Vector{Int}(undef,0)
        bond_to_definition    = Vector{Int}(undef,0)

        # declaring temporary storage vectors
        v′ = zeros(T2,Ndim)
        v″ = zeros(T2,Ndim)
        v‴ = zeros(T2,Ndim)

        # construct solver
        if lowercase(iterativesolver)=="cg"
            mul_by_M   = false
            transposed = false
            solver = ConjugateGradient(v″,tol=tol,maxiter=maxiter)
        elseif lowercase(iterativesolver)=="gmres"
            mul_by_M   = true
            transposed = false
            solver = GMRES(v″,tol=tol,restart=restart,maxiter=maxiter)
        elseif lowercase(iterativesolver)=="bicgstab"
            mul_by_M   = true
            transposed = false
            solver = BiCGStab(v″,tol=tol,maxiter=maxiter)
        end
        T3 = typeof(solver)

        # type of random number generator
        T4  = typeof(rng)

        return new{T1,T2,T3,T4}(β,Δτ,Lτ,Nsites,Nbonds,Nph,Ndim,Ndof,nbonds,nph,
                                lattice, bond_definitions,
                                x,t,ω,ω₄,α,α₂,μ,expΔτμ,t′,
                                field_to_phonon,field_to_τ,primary_field,
                                phonon_to_bond, bond_to_phonon,bond_to_definition,
                                checkerboard_perm,inv_checkerboard_perm,neighbor_table,cosht,sinht,
                                mul_by_M, transposed, v′, v″, v‴, solver,rng,"")
    end
end

"""
Assign hopping (and phonon) in SSH model.
"""
function assign_hopping!(ssh::SSHModel{T1,T2},t::T2,σt::T1,ω::T1,σω::T1,ω₄::T1,σω₄::T1,α::T1,σα::T1,α₂::T1,σα₂::T1,
                         o₁::Int,o₂::Int,v::Vector{Int},name::String="") where {T1<:AbstractFloat,T2<:Continuous}

    # define ssh bond
    sshbond = SSHBond(t,σt,ω,σω,ω₄,σω₄,α,σα,α₂,σα₂,o₁,o₂,v,name)
    push!(ssh.bond_definitions,sshbond)

    return nothing
end

"""
Assign chemical potential in SSH model.
"""
function assign_μ!(ssh::SSHModel,μ::T1,σμ::T1,orbit::Int) where {T1<:AbstractFloat}

    site_to_orbit = ssh.lattice.site_to_orbit::Vector{Int}
    for i in 1:ssh.Nsites
        if orbit==0 || orbit==site_to_orbit[i]
            ssh.μ[i]      = μ + σμ * randn(ssh.rng)
            ssh.expΔτμ[i] = exp( ssh.Δτ * ssh.μ[i] )
        end
    end

    return nothing
end

"""
Initialize SSH Model.
"""
function initialize_model!(ssh::SSHModel{T1,T2}) where {T1,T2}

    # number of types of bonds in lattice
    ssh.nbonds = length(ssh.bond_definitions)

    # initially zero types of phonons
    ssh.nph = 0

    # names of phonons
    names = Vector{String}(undef,0)

    # iterate over types of bonds in lattices
    for i in 1:ssh.nbonds

        # getting parameters that define bond
        @unpack t, σt, α, σα, α₂, σα₂, ω, σω, ω₄, σω₄, o₁, o₂, v, has_phonon, name = ssh.bond_definitions[i]

        # calculate new neighbors
        new_neighbors = calc_neighbor_table(ssh.lattice,o₁,o₂,v)

        # number of new bonds
        Nnewbonds = size(new_neighbors,2)

        # adding new neighbors to neighbor tables for current bond defintion
        ssh.neighbor_table = hcat(ssh.neighbor_table, new_neighbors)

        # phase of hopping
        if iszero(t)
            phase = 1.0
        else
            phase = t/abs(t)
        end

        # adding new hopppings for current bond
        t_new = @. phase * ( fill(abs(t),Nnewbonds) + σt*randn(ssh.rng,Nnewbonds) )
        append!(ssh.t, t_new)

        # adding bond to definition mapping
        append!(ssh.bond_to_definition, fill(i,Nnewbonds))

        # if bond has phonon on it
        if has_phonon

            # incrementing count of number of types of phonons
            ssh.nph += 1

            # record phonon names
            push!(names,name)

            # adding new phonon frequencies
            ω_new = @. fill(ω,Nnewbonds) + σω*randn(ssh.rng,Nnewbonds)
            append!(ssh.ω, ω_new)

            # adding anharmonic phonon coefficient
            ω₄_new = @. fill(ω₄,Nnewbonds) + σω₄*randn(ssh.rng,Nnewbonds)
            append!(ssh.ω₄, ω₄_new)

            # adding new linear electron-phonon coupling
            if iszero(α)
                phase = 1.0
            else
                phase = α/abs(α)
            end
            α_new = @. phase * ( fill(abs(α),Nnewbonds) + σα*randn(ssh.rng,Nnewbonds) )
            append!(ssh.α, α_new)

            # adding new non-linear electron-phonon coupling
            if iszero(α₂)
                phase = 1.0
            else
                phase = α₂/abs(α₂)
            end
            α₂_new = @. phase * ( fill(abs(α₂),Nnewbonds) + σα₂*randn(ssh.rng,Nnewbonds) )
            append!(ssh.α₂, α₂_new)

            # adding phonon to bond mapping
            append!(ssh.phonon_to_bond, collect((i-1)*Nnewbonds+1:i*Nnewbonds))

            # adding bond to phonon mapping
            append!(ssh.bond_to_phonon, collect((ssh.nph-1)*Nnewbonds+1:ssh.nph*Nnewbonds))
        else

            # add null mapping from bond to no phonon
            append!(ssh.bond_to_phonon, zeros(Int,Nnewbonds))
        end
    end

    # sort neighbor_table
    perm                = sorted_neighbor_table_perm!(ssh.neighbor_table)
    ssh.neighbor_table .= ssh.neighbor_table[:,perm]

    # get checkerboard groups
    groups = checkerboard_groups(ssh.neighbor_table)

    # get checkerboard permutation and apply to neighbor table
    new_perm                  = checkerboard_order(groups)
    ssh.neighbor_table       .= ssh.neighbor_table[:,new_perm]
    ssh.inv_checkerboard_perm = perm[new_perm]
    ssh.checkerboard_perm     = sortperm(ssh.inv_checkerboard_perm)

    # number of bonds in lattice
    ssh.Nbonds = size(ssh.neighbor_table,2)

    # number of phonons in lattice
    ssh.Nph = length(ssh.ω)

    # number of phonon fields/degrees of freedom
    ssh.Ndof = ssh.Lτ * ssh.Nph

    # initialize phonon fields
    ssh.x = zeros(T1,ssh.Ndof)

    # mapping from field to phonon
    ssh.field_to_phonon = repeat(collect(1:ssh.Nph), inner=ssh.Lτ)
    ssh.field_to_τ      = repeat(collect(1:ssh.Lτ),  outer=ssh.Nph)

    # intialize checkerboard matrix elements
    ssh.t′    = zeros(T2,ssh.Lτ,ssh.Nbonds)
    ssh.cosht = zeros(T2,ssh.Lτ,ssh.Nbonds)
    ssh.sinht = zeros(T2,ssh.Lτ,ssh.Nbonds)
    for bond in 1:ssh.Nbonds
        t′    = ssh.t[bond]
        cosht = cosh(ssh.Δτ*t′)
        sinht = sinh(ssh.Δτ*t′)
        index = ssh.checkerboard_perm[bond]
        for τ in 1:ssh.Lτ
            ssh.t′[τ,bond]     = t′
            ssh.cosht[τ,index] = cosht
            ssh.sinht[τ,index] = sinht
        end
    end

    # initiailize primary fields mapping
    primary_field  = collect(1:ssh.Ndof)
    primary_field′ = reshaped(primary_field,(div(ssh.Ndof,ssh.nph),ssh.nph))

    # construct a tabulation of all the fields
    fields = reshaped(collect(1:ssh.Ndof),(div(ssh.Ndof,ssh.nph),ssh.nph))

    # iterate over types of phonons
    for ph_type in 1:ssh.nph
        # iterate over remaining types of phonons
        for ph_type′ in (ph_type+1):ssh.nph
            # if phonons share the same name
            if names[ph_type]==names[ph_type′]
                # record primary field
                if primary_field′[1,ph_type′] > fields[1,ph_type]
                    @views @. primary_field′[:,ph_type′] = fields[:,ph_type]
                end
            end
        end
    end

    # record primary fields
    ssh.primary_field = primary_field

    return nothing
end

"""
Update SSH model.
"""
function update_model!(ssh::SSHModel{T1,T2}) where {T1,T2}

    # account of updated chemical potential
    @. ssh.expΔτμ = exp( ssh.Δτ * ssh.μ )

    # warning flag for unphysically large phonon displacement
    flag_large_displacement = false

    # iterate of phonon fields
    for field in 1:ssh.Ndof
        # get phonon
        phonon = ssh.field_to_phonon[field]
        # get time slice
        τ = ssh.field_to_τ[field]
        # get bond
        bond = ssh.phonon_to_bond[phonon]
        # getting indexing for checkerboard order
        index = ssh.checkerboard_perm[bond]
        # get phonon field value
        xτ = ssh.x[field]
        # update matrix elements of exp{-Δτ⋅K} = exp{Δτ⋅(t-α⋅x)}
        v                  = ssh.α[phonon]*xτ + sign(xτ)*ssh.α₂[phonon]*xτ^2
        t′                 = ssh.t[bond] - v
        ssh.t′[τ,bond]     = t′
        ssh.cosht[τ,index] = cosh(ssh.Δτ*t′)
        ssh.sinht[τ,index] = sinh(ssh.Δτ*t′)
        # detect unphysically large phonon displacement
        if abs(ssh.t[bond]) < abs(v)
            flag_large_displacement = true
        end
    end

    # report unphysically large phonon displacement
    if flag_large_displacement
        @info("Unphysically Large Phonon Displacement\n")
        logger = global_logger()
        flush(logger.stream)
    end

    # make sure equivalent fields are equal
    for field in 1:ssh.Ndof
        # get primary field
        field′ = ssh.primary_field[field]
        # if current field is not primary
        if field != field′
            if !(ssh.x[field]≈ssh.x[field′])
                error("(x[$field]=$(ssh.x[field])) != (x[$field′]=$(ssh.x[field′]))\n")
            end
        end
    end

    return nothing
end

"""
Fill vector that will ultimate be used to update fields in model.
"""
function randn!(v::AbstractVector{T},ssh::SSHModel{T}) where {T}

    # fill with random numbers
    randn!(ssh.rng,v)

    # account for some fields being equivalent
    @views @. v = v[ssh.primary_field]

    return nothing
end

"""
Multiply by M matrix.
"""
function mulM!(Mv::AbstractVector{T2},ssh::SSHModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    ####################################
    ## PERFORM MULTIPLICATION y = M⋅v ##
    ####################################

    #           EXAMPLE OF M MATRIX CONVENTION
    #     |   1     0      0      0      0   +B(1) |
    #     | -B(2)   1      0      0      0     0   |
    # M = |   0   -B(3)    1      0      0     0   |
    #     |   0     0    -B(4)    1      0     0   |
    #     |   0     0      0    -B(5)    1     0   |
    #     |   0     0      0      0    -B(6)   1   |

    # Notes:
    # • [M⋅v](1) = v(1) + B(1)⋅v(Lτ)  for τ = 1
    # • [M⋅v](τ) = v(τ) - B(τ)⋅v(τ-1) for τ > 1
    # • B(τ) = exp{-Δτ⋅K[x(τ)]}⋅exp{Δτ⋅μ}
    # • exp{Δτ⋅μ} is represent an NxN diagonal matrix that is the same for all
    #   time slice τ and is stored as a vector.
    # • exp{-Δτ⋅K[x(τ)]} is applied to a vector/represented using the checkerboard approximation.

    @unpack cosht, sinht, neighbor_table, Nbonds, Nsites, Lτ, expΔτμ = ssh

    # Calculate [M⋅v](τ) = exp{Δτ⋅μ}⋅v(τ-1)

    @fastmath @inbounds for i in 1:Nsites
        expΔτμi = expΔτμ[i]
        for τ in 1:Lτ
            τm1    = mod1(τ-1,Lτ)
            iτ     = get_index(τ,  i,Lτ)
            iτm1   = get_index(τm1,i,Lτ)
            Mv[iτ] = expΔτμi * v[iτm1]
        end
    end

    # Calculate [M⋅v](τ) = B(τ)⋅v(τ-1) = exp{-Δτ⋅K[x(τ)]}⋅[exp{Δτ⋅μ}⋅v(τ-1)]

    checkerboard_mul!(Mv,neighbor_table,cosht,sinht)

    # Calculate [M⋅v](τ) = v(τ) - B(τ)⋅v(τ-1)

    # iterate over sites in lattice
    @fastmath @inbounds for i in 1:Nsites

        # [M⋅v](1) = v(1) + B(1)⋅v(Lτ)
        i1     = get_index(1,i,Lτ)
        Mv[i1] = v[i1] + Mv[i1]

        # iterate of time slices
        for τ in 2:Lτ

            # [M⋅v](τ) = v(τ) - B(τ)⋅v(τ-1)
            iτ     = get_index(τ,i,Lτ)
            Mv[iτ] = v[iτ] - Mv[iτ]
        end
    end

    return nothing
end

"""
Multiply by Mᵀ matrix.
"""
function mulMᵀ!(Mᵀv::AbstractVector{T2},ssh::SSHModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    #####################################
    ## PERFORM MULTIPLICATION y = Mᵀ⋅v ##
    #####################################

    #           EXAMPLE OF Mᵀ MATRIX CONVENTION
    #      |    1   -Bᵀ(2)   0      0      0      0    |
    #      |    0     1    -Bᵀ(3)   0      0      0    |
    # Mᵀ = |    0     0      1    -Bᵀ(4)   0      0    |
    #      |    0     0      0      1    -Bᵀ(5)   0    |
    #      |    0     0      0      0      1    -Bᵀ(6) |
    #      | +Bᵀ(1)   0      0      0      0      1    |

    # Notes:
    # • y(τ)  = [Mᵀ⋅v](τ)  = v(τ)  - Bᵀ(τ+1)⋅v(τ+1) for τ < Lτ
    # • y(Lτ) = [Mᵀ⋅v](Lτ) = v(Lτ) + Bᵀ(1)⋅v(1)     for τ = Lτ
    # • Bᵀ(τ) = exp{Δτ⋅μ}⋅exp{-Δτ⋅K[x(τ)]}ᵀ
    # • exp{Δτ⋅μ} is represent an NxN diagonal matrix that is the same for all
    #   time slice τ and is stored as a vector.
    # • exp{-Δτ⋅K[x(τ)]} is applied to a vector/represented using the checkerboard approximation.

    @unpack cosht, sinht, neighbor_table, Nbonds, Nsites, Lτ, expΔτμ = ssh

    # Calculate exp{-Δτ⋅K[x(τ+1)]}ᵀ⋅v(τ+1)

    copyto!(Mᵀv,v)
    checkerboard_transpose_mul!(Mᵀv,neighbor_table,cosht,sinht)

    # Calculate [Mᵀ⋅v](τ) = v(τ) - Bᵀ(τ+1)⋅v(τ+1)

    # iterate over sites in lattice
    for i in 1:Nsites

        expΔτμi = expΔτμ[i]

        # record [Mᵀ⋅v](L) = v(L) + Bᵀ(1)⋅v(1) = v(L) + exp{Δτ⋅μ}⋅exp{-Δτ⋅K[x(1)]}ᵀ⋅v(1)
        i1     = get_index(1 ,i,Lτ)
        iL     = get_index(Lτ,i,Lτ)
        Mᵀv_iL = v[iL] + expΔτμi * Mᵀv[i1]

        # iterate over time slices
        for τ in 1:Lτ-1

            # [Mᵀ⋅v](τ) = v(τ) - Bᵀ(τ+1)⋅v(τ+1) = v(τ) - exp{Δτ⋅μ}⋅exp{-Δτ⋅K[x(τ+1)]}ᵀ⋅v(τ+1)
            iτ      = get_index(τ,i,Lτ)
            iτp1    = get_index(τ+1,i,Lτ)
            Mᵀv[iτ] = v[iτ] - expΔτμi * Mᵀv[iτp1]
        end

        # [Mᵀ⋅v](L) = v(L) + Bᵀ(1)⋅v(1)
        Mᵀv[iL] = Mᵀv_iL
    end

    return nothing
end

"""
Calculates ⟨∂M/∂xᵢⱼ(τ)⟩ = uᵀ⋅[∂M/∂xᵢⱼ(τ)]⋅v for each xᵢⱼ(τ) degree of freedom,
writing each of the ⟨∂M/∂xᵢⱼ(τ)⟩ values to a vector dMdx.
"""
function muldMdx!(dMdx::AbstractVector{T2},u::AbstractVector{T2},ssh::SSHModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    #           EXAMPLE OF dM/dx MATRIX CONVENTION
    #          |     0        0        0        0        0       0    |
    #          |     0        0        0        0        0       0    |
    # dM/dx₃ = |     0    -dB₃/dx₃     0        0        0       0    |
    #          |     0        0        0        0        0       0    |
    #          |     0        0        0        0        0       0    |
    #          |     0        0        0        0        0       0    |

    # Notes:
    # • B(τ)          = exp{-Δτ⋅K(τ)}⋅exp{Δτ⋅μ}
    # • ∂B(τ)/∂xᵢⱼ(τ) = -Δτ⋅∂K(τ)/∂xᵢⱼ(τ)⋅exp{-Δτ⋅K(τ)}⋅exp{Δτ⋅μ} = -Δτ⋅∂K(τ)/∂xᵢⱼ(τ)⋅B(τ)
    # • Kᵢⱼ(τ)        = -(tᵢⱼ-αᵢⱼ⋅xᵢⱼ(τ)) = matrix elements of a symmetric NxN matrix K(τ)
    # • ∂K(τ)/∂xᵢⱼ(τ) = a symmetric NxN matrix with two non-zero matrix elements given by ∂Kᵢⱼ(τ)/∂xᵢⱼ(τ) = αᵢⱼ
    #
    # • ⟨∂M/∂xᵢⱼ(τ)⟩ = uᵀ⋅[∂M/∂xᵢⱼ(τ)]⋅v                       for τ > 1
    #                = uᵀ(τ)⋅[-∂B(τ)/∂xᵢⱼ(τ)]⋅v(τ-1)           for τ > 1
    #                = uᵀ(τ)⋅[Δτ⋅∂K(τ)/∂xᵢⱼ(τ)⋅B(τ)]⋅v(τ-1)    for τ > 1
    #                = +uᵀ(τ)⋅[Δτ⋅∂K(τ)/∂xᵢⱼ(τ)]⋅[B(τ)⋅v(τ-1)] for τ > 1
    #
    # • ⟨∂M/∂xᵢⱼ(1)⟩ = uᵀ⋅[∂M/∂xᵢⱼ(1)]⋅v                     for τ = 1
    #                = uᵀ(1)⋅[+∂B(1)/∂xᵢⱼ(1)]⋅v(L)           for τ = 1
    #                = uᵀ(1)⋅[-Δτ⋅∂K(1)/∂xᵢⱼ(τ)⋅B(1)]⋅v(L)   for τ = 1
    #                = -uᵀ(1)⋅[Δτ⋅∂K(1)/∂xᵢⱼ(τ)]⋅[B(1)⋅v(L)] for τ = 1

    # Checkerboard Decomposition:           exp{-Δτ⋅K(τ)}  ≈ ∏(n=1...N)[exp{-Δτ⋅Kₙ(τ)}]
    # Transpose Checkerboard Decomposition: exp{-Δτ⋅K(τ)}ᵀ ≈ ∏(n=N...1)[exp{-Δτ⋅Kₙ(τ)}]
    # Note: exp{-Δτ⋅Kₙ(τ)}ᵀ = exp{-Δτ⋅Kₙ(τ)} [Hermitian]

    @unpack cosht, sinht, neighbor_table, checkerboard_perm, inv_checkerboard_perm = ssh
    @unpack Nbonds, Nph, Ndof, Nsites, Lτ, Δτ, expΔτμ, α, α₂, x = ssh
    @unpack field_to_phonon, field_to_τ, phonon_to_bond, bond_to_phonon = ssh
    @unpack primary_field = ssh

    b = ssh.v′
    c = ssh.v″

    # |b₀(τ)⟩ = exp{Δτ⋅μ}|v(τ-1)⟩
    @fastmath @inbounds for i in 1:Nsites
        expΔτμi = expΔτμ[i]
        for τ in 1:Lτ
            τm1   = mod1(τ-1,Lτ)
            iτ    = get_index(τ,  i,Lτ)
            iτm1  = get_index(τm1,i,Lτ)
            b[iτ] = expΔτμi * v[iτm1]
        end
    end

    # ⟨c₀(τ)| = ⟨u(τ)|exp{-Δτ⋅K(τ)} ⟹ |c₀(τ)⟩ = exp{-Δτ⋅K(τ)}ᵀ|u(τ)⟩
    copyto!(c,u)
    checkerboard_transpose_mul!(c,neighbor_table,cosht,sinht)

    # initialize dMdx to zero
    fill!(dMdx,0.0)

    # iterate over bonds in checkerboard order
    for n in 1:Nbonds

        # get bond
        bond = inv_checkerboard_perm[n]

        # get phonon
        phonon = bond_to_phonon[bond]

        # get pair of sites
        i = neighbor_table[1,n]
        j = neighbor_table[2,n]

        # iterate over imaginary time
        for τ in 1:Lτ

            # get matrix elements of exp{-Δτ⋅Kₙ(τ)}
            cosht′ = cosht[τ,n]
            sinht′ = sinht[τ,n]

            # get indices into vectors
            iτ = get_index(τ,i,Lτ)
            jτ = get_index(τ,j,Lτ)

            # |bₙ(τ)⟩ = exp{-Δτ⋅Kₙ(τ)}|bₙ₋₁(τ)⟩
            biτ   = b[iτ]
            bjτ   = b[jτ]
            b[iτ] = cosht′ * biτ +      sinht′  * bjτ
            b[jτ] = cosht′ * bjτ + conj(sinht′) * biτ

            # |cₙ(τ)⟩ = exp{-Δτ⋅Kₙ(τ)}⁻¹|cₙ₋₁(τ)⟩
            ciτ   = c[iτ]
            cjτ   = c[jτ]
            c[iτ] = cosht′ * ciτ -      sinht′  * cjτ
            c[jτ] = cosht′ * cjτ - conj(sinht′) * ciτ

            # if there is a phonon on the bond
            if phonon != 0

                # get index to current field
                field = get_index(τ,phonon,Lτ)

                # get xₙ(τ)
                xₙ = x[field]

                # calculate ∂Kₙ(τ)/∂xₙ(τ)
                dKdx = α[phonon] + 2*α₂[phonon]*xₙ

                # ⟨∂M/∂xᵢⱼ(τ)⟩ = ⟨cₙ(τ)|Δτ⋅∂Kₙ(τ)/∂xₙ(τ)|bₙ(τ)⟩
                dmdx = conj(c[jτ])*Δτ*dKdx*b[iτ] + conj(c[iτ]*Δτ*dKdx)*b[jτ]

                # flip sign if τ=1
                if τ==1
                    dmdx = -dmdx
                end

                # record result
                dMdx[primary_field[field]] += dmdx
            end
        end
    end

    # account for equivalent fields
    @views @. dMdx = dMdx[primary_field]

    return nothing
end

##############################
## PHONON FIELD IO ROUTINES ##
##############################

"""
Writes the current phonon field configuration to file.
"""
function write_phonons!(ssh::SSHModel{T1,T2,T3},filename::String) where {T1,T2,T3}

    if ssh.Nph>0

        # get info about size of lattice and number of phonons
        L  = ssh.Lτ
        n  = ssh.nph
        N  = div(ssh.Nph,n)

        # get phonon fields
        x = reshaped(ssh.x,(L,N,n))

        # open file
        open(filename,"w") do file

            # write header to file
            write(file, "type loc tau x\n")

            # iterate over phonon fileds
            for phonon in 1:n
                for i in 1:N
                    for τ in 1:L
                        # get phonon field
                        x₀ = x[τ,i,phonon]
                        # write to file
                        write(file, @sprintf("%d %d %d %.6f\n",phonon,i,τ,x₀))
                    end
                end
            end
        end
    end

    return nothing
end

"""
Read phonon config from file.
"""
function read_phonons!(ssh::SSHModel{T1,T2,T3},filename::String) where {T1,T2,T3}

    # get info about size of lattice and number of phonons
    L  = ssh.Lτ
    n  = ssh.nph
    N  = div(ssh.Nph,n)

    # get phonon fields
    x = reshaped(ssh.x,(L,N,n))

    # open file
    open(filename,"r") do file

        # read in header
        header = readline(file)

        # iterate over lines in file
        for line in eachline(file)

            # split line at white space
            atoms = split(line," ")

            # extract info about phonon field and location
            phonon = parse(Int,atoms[1])
            cell   = parse(Int,atoms[2])
            τ      = parse(Int,atoms[3])
            x₀     = parse(T1, atoms[4])

            # record phonon field
            x[τ,cell,phonon] = x₀
        end
    end

    # construct exponentiated interaction matrix
    update_model!(ssh)

    return nothing
end

"Write Hamiltonian K[τ] matrix for given time slice τ to file. Includes on-site energies."
function write_K_matrix!(filename::String, model::SSHModel{T1,T2,T3}, τ::Int=1) where {T1,T2,T3}

    open(filename,"w") do file

        # write header
        write(file,"col row val\n")

        # iterate over on-site energies
        for i in 1:model.Nsites
            ϵ = -model.μ[i]
            write(file,"$i $i $ϵ\n")
        end

        # iterate over bonds
        for bond in 1:model.Nbonds
            # matrix element based on hopping amplitude
            val = -model.t′[τ,bond]
            # neighboring sites associated with bond
            indx = model.checkerboard_perm[bond]
            i = model.neighbor_table[1,indx]
            j = model.neighbor_table[2,indx]
            # write to file
            write(file,"$i $j $val\n")
            write(file,"$j $i $(conj(val))\n")
        end
    end

    return nothing
end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
fp = dirname(@__FILE__)

@testset "VAR impulse responses" begin
    data = load(joinpath(fp, "../../reference/var_irf_out.jld2"))
    β = data["bhat_full"]
    Σ = data["sigmahat_full"]
    irfs1 = data["irfs_full_choleskyLR"]
    irfs2 = data["irfs_full_maxBC"]
    irfs3 = data["irfs_full_cholesky"]
    out1 = impulse_responses(β, Σ, 1, 20;
                             method = :choleskyLR,
                             flip_shocks = true)
    out2 = impulse_responses(β, Σ, 1, 20;
                             method = :maxBC,
                             flip_shocks = true)
    out3 = impulse_responses(β, Σ, 1, 20;
                             method = :cholesky,
                             flip_shocks = true)
    out4 = impulse_responses(β, Σ, 1, 20;
                             method = :cholesky_long_run,
                             flip_shocks = true)
    out5 = impulse_responses(β, Σ, 1, 20;
                             method = :maximum_business_cycle_variance,
                             flip_shocks = true)
    out6 = impulse_responses(β, Σ, 1, 20;
                             method = :choleskyLR,
                             flip_shocks = false)
    out7 = impulse_responses(β, Σ, 1, 20;
                             method = :maxBC,
                             flip_shocks = false)
    out8 = impulse_responses(β, Σ, 1, 20;
                             method = :cholesky,
                             flip_shocks = false)

    @test @test_matrix_approx_eq irfs1' out1
    @test @test_matrix_approx_eq irfs2' out2
    @test @test_matrix_approx_eq irfs3' out3
    @test @test_matrix_approx_eq irfs1' out4
    @test @test_matrix_approx_eq irfs2' out5
    @test @test_matrix_approx_eq irfs1[1, :] -out6[:, 1]
    @test @test_matrix_approx_eq irfs2[1, :] -out7[:, 1]
    @test @test_matrix_approx_eq irfs3[1, :] -out8[:, 1]
end

@testset "VECM impulse responses" begin
    data = load(joinpath(fp, "../../reference/dsgevecm_lambda_irfs.jld2"))
    expout = load(joinpath(fp, "../../reference/vecm_irfs.jld2"), "expout")
    β = data["cct_sim"]
    Σ = data["sig_sim"]
    coint_mat = data["cointvec"]
    out1 = impulse_responses(β, Σ, coint_mat, 1, 5;
                             method = :cholesky,
                             flip_shocks = true,
                             use_intercept = true)
    out2 = impulse_responses(β, Σ, coint_mat, 1, 5;
                             method = :cholesky,
                             flip_shocks = false,
                             use_intercept = true)
    out3 = impulse_responses(β[vcat(1:3, 5:32), :], Σ, coint_mat, 1, 5;
                             method = :cholesky, use_intercept = false,
                             flip_shocks = true)
    out4 = impulse_responses(β[vcat(1:3, 5:32), :], Σ, coint_mat, 1, 5;
                             method = :cholesky, use_intercept = false,
                             flip_shocks = false)

    @test_throws ErrorException impulse_responses(β, Σ, coint_mat, 1, 5;
                                                  method = :not_a_method,
                                                  use_intercept = true)

    @test @test_matrix_approx_eq expout out1
    @test @test_matrix_approx_eq expout -out2
    @test @test_matrix_approx_eq expout out3
    @test @test_matrix_approx_eq expout -out4
end

nothing


text:
using Distributed
addprocs(3)
@everywhere import Primes
ref = @spawn Primes.factor(21883298135690819)
factors = fetch(ref)
@assert collect(factors) == [93234719=>1, 234711901=>1]

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


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

name = "GSL"
version = v"2.7.1" # <--- This version number is a lie to build for experimental platforms

# Collection of sources required to build GSL
sources = [
    ArchiveSource("http://ftp.gnu.org/gnu/gsl/gsl-$(version.major).$(version.minor).tar.gz",
                  "efbbf3785da0e53038be7907500628b466152dbc3c173a87de1b5eba2e23602b"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gsl-*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgslcblas", :libgslcblas),
    LibraryProduct("libgsl", :libgsl),
    ExecutableProduct("gsl-histogram", :gsl_histogram),
    ExecutableProduct("gsl-randist", :gsl_randist),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Plots
using SolidStateDetectors

T = Float32
sim = Simulation(SSD_examples[:InvertedCoax])

calculate_electric_potential!(sim, max_refinements = 0)
plot(sim.ρ)

# The charge density model (actually it should be named impurity density model, we have to change this at some point)
# is part of the semiconductor object of the detector:
sim.detector.semiconductors[1].charge_density_model

# It can be overwritten by an custom density model.
# This must be subtype of AbstractChargeDensity model:

struct CustomChargeDensity{T} <: SolidStateDetectors.AbstractChargeDensity{T}
    # Define all the parameters for your own model here
    # e.g.:
    offset::T
    z_gradient::T
    r_gradient::T
end

# Create a model with arbitray parameters
ccdm = CustomChargeDensity{Float32}( -10^16, -10^18, -10^18) # meter not mm

# overwrite the predefined model
sim.detector.semiconductors[1].charge_density_model = ccdm

# We also need to define two methods for this density model:
# (So actually the model)
function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::CylindricalPoint{T}) where {T}
    # This is the actually the impurity density: it will be multiplied by the elementary charge (in SSD)
    # Just an example:
    charge_density::T = cdm.offset
    charge_density += cdm.z_gradient * pt[3] # z coordinate
    charge_density += cdm.r_gradient * pt[1] # r coordinate
    # Note, that pt will be parsed to this function in SI units -> meter
    return charge_density
end
# In case you use cartesian coordinates for the simulation just add this method:
function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::CartesianPoint{T}) where {T}
    return SolidStateDetectors.get_charge_density(cdm, CylindricalPoint(pt))
end

calculate_electric_potential!(sim, max_refinements = 2, depletion_handling=true)

plot(sim.point_types)
plot(sim.ρ)

# Plot the profile:
# We want to add this as a plot recipe in the future but is not yet implemented.
rs = 0.0:1e-5:4e-2
charge_vec = similar(rs)
for ir in eachindex(rs)
    pt = CylindricalPoint{T}(rs[ir], 0, 1e-3)
    charge_vec[ir] = pt in sim.detector.semiconductors[1] ?
        SolidStateDetectors.get_charge_density(sim.detector.semiconductors[1], pt) :
        0
end
plot(rs, charge_vec, ylabel = "impurity density [1/m^3]", xlabel = "r [m]")

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt
using Test
using KCenters

include("kcenters.jl")


text:
using DCELs
using Test

@testset "Unit tests" begin
    include("connecting.jl")
    include("errors.jl")
    include("geometry.jl")
    include("iterators.jl")
    include("operations.jl")
    include("printing.jl")
    include("sorting.jl")
    include("squeezing.jl")
end


text:
module FEPreProcessing
    __precompile__()
    using LinearAlgebra
    #
    using StructArrays
    using StaticArrays
    using SparseArrays
    using Parameters
    #
    using Materials

    const IF64 = Union{Int64, Float64}

    include("FEPrepBase.jl")
    include("mesh_data.jl")
    include("assembly.jl")
    include("prepare_material_data.jl")
    #
    include("UnitCellFEModel.jl")
    export make_unit_cell_FEA_model
    #
end


text:
push!(LOAD_PATH, "./models")
push!(LOAD_PATH, ".")
module Utils
using SparseInverseProblemsMod
using SuperResModels
threshold_weight = 1e-1
export run_simulation, generate_and_reconstruct_all

using PyCall
@pyimport numpy as np

function phi_noise(model, thetas, weights, sigma)
    return sum([weights[i] * psi_noise(model, vec(thetas[:, i]), sigma) for i in 1:size(thetas, 2)])
end
function psi_noise(model :: DynamicSuperRes, theta :: Vector{Float64}, dsec :: Float64)
    # This function computes the direct problem for a single point theta, with
    # and added acceleration that is proportional to the speed by dsec value.
    # aceleration * tau K / speed = dsec.
    d = dim(model)
    tauK = maximum(model.times)
    return vcat([psi(model.static, theta[1:d] + t * theta[d+1:end] + t^2*dsec*theta[d+1:end]/tauK/2) for t in model.times]...)
end

function run_simulation_target(model, thetas, weights, target, algIter)
# We simulate the reconstruction algorithm on the data given by target, and we return
# the error values from the background truth.
    function callback(old_thetas, thetas, weights, output, old_obj_val)
       #evalute current OV
       new_obj_val,t = loss(LSLoss(), output - target)
       #println("gap = $(old_obj_val - new_obj_val)")
       if old_obj_val - new_obj_val < 1E-4
           return true
       end
    return false
    end
    # Inverse problem
    (thetas_est,weights_est) = ADCG(model, LSLoss(), target,sum(weights), callback=callback, max_iters=algIter)
    return (thetas_est, weights_est)
end

function generate_target(model, thetas, weights, noise_level = 0.0, noise_position = 0.0)
    # Noise is added to the simulations. Noise_level is a percentage of the total
    # signal in L2 norm, in other words, it is the 1/signal-to-noise.
    if noise_position > 0
        target =  phi_noise(model, thetas, weights, noise_position)
    else
        target = phi(model, thetas, weights)
    end
    # Noise is bounded by 1 in modulus
    noise = randn(size(target))*noise_level
    # Each particle's signal is bounded by it's weights, we set them to 1 ideally
    if (maximum(weights)>2 || minimum(weights)<0.1) && noise_level > 0
	warn(" The proportionality of noise is not accurate as the weights are too high or low ")
    end
    # So a porcentage is obtained by just ponderating the noise.
    return target + noise
end

function match_points(theta_1, theta_2)
    # Method to match the reconstruction's particles with the given ones.
    n_points = size(theta_1, 2)
    corres = zeros(Int, n_points)
    distances = zeros(n_points,n_points)
    # to get from one dimensional array to the indexes in two dimensions.
    function iindex(x,N)
	return mod(x-1,N)+1
    end
    function jindex(x,N)
	return div(x-iindex(x,N),N)+1
    end
    # distances matrix.
    for i = 1:n_points
        for j = 1:n_points
            distances[i,j] = norm(theta_1[:,i] - theta_2[:,j])
        end
    end
    # We take the closest ones and match them
    for k = 1:n_points
	x = indmin(distances)
	(i,j) = (iindex(x,n_points), jindex(x,n_points))
	corres[i]=j
	distances[i,:] = Inf
	distances[:,j] = Inf
    end
    return corres
end

function to_static(thetas, t, x_max)
    # For a moving particle and a time, gives its static position.
    d = div(size(thetas, 1),2)
    pts = thetas[1:d,:]
    velocities = thetas[d+1:2*d,:]
    pts_t = pts + velocities * t
    pts_t = mod.(pts_t, x_max)
    return pts_t
end

function target_to_static(target_dynamic, times)
    # For time dependent measurements, return a container with each static measurement.
    static_size = div(length(target_dynamic), length(times))
    static_target = []
    for k = 1:length(times)
	push!(static_target, target_dynamic[static_size*(k-1) + 1: static_size*k])
    end
    return static_target
end

function separation_val(thetas,K,tau,x_max)
# function that evaluates the separation between the particles
# We project into the time measurements, and we take the third best
# separation.
separations = zeros(2*K+1)
# We also consider toroidal modulus distance on the interval
toroidalSeparation = zeros(2*K+1)
for k in 1:2*K+1
   separations[k] = minimum([abs(thetas[1,i]-thetas[1,j]+tau*(k-K-1)*(thetas[2,i]-thetas[2,j])) + (i==j)*x_max for i in 1:size(thetas,2), j in 1:size(thetas,2)])
   toroidalSeparation[k] = minimum([abs(thetas[1,i]+tau*(k-K-1)*(thetas[2,i]) - (thetas[1,j]+tau*(k-K-1)*(thetas[2,j])-x_max)) + (i==j)*x_max for i in 1:size(thetas,2), j in 1:size(thetas,2)])
end
separations = [min(separations[k],toroidalSeparation[k]) for k in 1:length(separations)]
sort!(separations)
return separations[end-2]
end

function separation_norm(thetas,K,tau,x_max)
# function that evaluates the separation between the particles, but in terms
# of the dynamic norm.
delta = K*tau
norm = (x,v) -> abs(x) + delta*abs(v)
return minimum([norm(thetas[1,i] - thetas[1,j], thetas[2,i] - thetas[2,j]) + (i==j)*Inf for i in size(thetas,2), j in 1:size(thetas,2) ])
end

function Rejection_sampling(test_case, bins, density, K, tau, x_max)
    # We use a rejection sampling algorithm to generate test_cases whose
    # separation is distributed uniformly. To save time in Montecarlo simulations.
    # Since we have no formula for the distribution of separations for random particles,
    # the distribution is estimated with separationDistribution.jl script
    # WARNING: This works only for cloud_1d_full in TestCases.jl.
    dx = bins[2]-bins[1]
    minDensity = minimum(density)
    densityFunc = y -> density[floor(Int,y/dx)+1]
    densityLength = length(density)
    safetyInd = 1
    (thetas, weights) = test_case()
    y = separation_val(thetas,K,tau,x_max)
    # Reject all the particules outside the desired interval.
    while floor(Int,y/dx)+1 >= densityLength
        (thetas, weights) = test_case()
        y = separation_val(thetas,K,tau,x_max)
    end
    # rejection-aceptation method
    reject_var = rand()
    while safetyInd<1e6 && reject_var >= minDensity/densityFunc(y)
        (thetas, weights) = test_case()
        y = separation_val(thetas,K,tau,x_max)
        while floor(Int,y/dx)+1 >= densityLength
            (thetas, weights) = test_case()
            y = separation_val(thetas,K,tau,x_max)
        end
        reject_var = rand()
	safetyInd = safetyInd + 1
    end
    if safetyInd == 1e6
        error(" Something is going wrong with the rejection algorithm ! We are
        not finding any appropiate sample ! ")
    end
    return (thetas, weights)
end

function generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, target, algIter)
    d = dim(model_dynamic)
    # run the algorithm on the points thetas with the respective noises.
    (thetas_est, weights_est) = try run_simulation_target(model_dynamic, thetas, weights, target, algIter) catch
	warn("fail")
	([0], [0])
	end
    # Discard any reconstruction with smaller weight than some threshold.
    thetas_est = thetas_est[:, weights_est .> threshold_weight]
    weights_est = weights_est[weights_est .> threshold_weight]
    if (length(thetas) == length(thetas_est))
        # The output of the minimization algorithm returns an unsorted vector with the
        # estimated positions ans weights of particles. match_points is a function that
        # returns a vector that tries to match the best possible the corresponding particles
        # to the estimated particles.
        corres = match_points(thetas, thetas_est)
        dist_x = norm(thetas[1:d, :] - thetas_est[1:d, corres], Inf)
        dist_v = norm(thetas[d+1:end, :] - thetas_est[d+1:end, corres], Inf)
        dist_w = norm(weights[:] - weights_est[corres], Inf)
        # returns the infinite distance between the matched particles, in space and velocity.
        return (dist_x, dist_v, dist_w)
    end
    # if the reconstructed particles don't match in size, we call this a failure and the
    # distance is the maximal given one.
    return model_dynamic.static.x_max, model_dynamic.v_max, maximum(weights)
end

function generate_and_reconstruct_static(model_static, thetas, weights, target_static, algIter)
    # Description is similar to the one in generate_and_reconstruct_dynamic.
    d = dim(model_static)
    (thetas_est, weights_est) = try
        run_simulation_target(model_static, thetas, weights, target_static, algIter)
    catch
        warn("fail")
        ([0], [0])
    end
    thetas_est = thetas_est[:, weights_est .> threshold_weight]
    weights_est = weights_est[weights_est .> threshold_weight]
    if (length(thetas) == length(thetas_est))
        corres = match_points(thetas, thetas_est)
        dist_x = norm(thetas[1:d, :] - thetas_est[1:d, corres], Inf)
        dist_w = norm(weights[:] - weights_est[corres], Inf)
        return dist_x, dist_w
    end
    return model_static.x_max, maximum(weights)
end

function generate_and_reconstuct_static_best(model_dynamic, model_static, thetas, weights, target_dynamic, algIter)
    # Given a dynamic set of particles, we obtain the respective static reconstructions
    # at each time sample, and afterwards we return the error in space and weight
    # for the case of the best one, and the third best one.
    dist_x = zeros(length(model_dynamic.times))
    dist_w = zeros(length(model_dynamic.times))
    static_target = target_to_static(target_dynamic, model_dynamic.times)
    for k = 1:length(model_dynamic.times)
        thetas_t = to_static(thetas, model_dynamic.times[k], model_static.x_max)
        (dist_x[k], dist_w[k]) = generate_and_reconstruct_static(model_static, thetas_t, weights, static_target[k], algIter)
    end
    sortOrder = sortperm(dist_x)
    dist_x = dist_x[sortOrder]
    dist_w = dist_w[sortOrder]
    # The first two outputs are the best case and its weight
    # the second two outpust are the third best case with the worse weight.
    return dist_x[1], dist_w[1], dist_x[3], maximum(dist_w[1:3])
end

function generate_and_reconstruct_all(model_static, model_dynamic, bins, density, algIter, test_case, noises_data, noises_position, cases)
    ### This function generates a test case and returns the distances between the generated particles
    ### given by dx, dv and norm. Furthermore it returns a result matrix that contains
    ### the first 3 colums correspond to the dynamic results, (space x velocity x weights)
    ### whereas the last 4 columns correspond to the static results-: best measurement location
    ### and weight, third best location, worse weight difference among the best 3.
    ### The rows are respectively no noise, noises in data, noises in position.
    (noiseless_dynamic, noiseless_static, noise_dynamic, noise_static, curvature_dynamic) = cases
    K = div(length(model_dynamic.times)-1,2)
    tau = maximum(model_dynamic.times)/K
    (thetas, weights) =  Rejection_sampling(test_case, bins, density, K, tau, model_static.x_max)
    d = dim(model_dynamic)
    results = zeros(1 + length(noises_data) + length(noises_position), 7)
    # Obtain the separation of the configuration of particles
    separation = separation_val(thetas,K,tau,model_static.x_max)
    # And the separationg given by the dynamic norm.
    separation_dyn = separation_norm(thetas,K,tau,model_static.x_max)
    # Dynamic case, no noise
    target = generate_target(model_dynamic, thetas, weights)
    println("######### STARTING ########")
    # We reconsctruct the location of the particles for the target data.
    if noiseless_dynamic == true
	println(" trying noiseless dynamical ###### ")
       (results[1, 1], results[1, 2], results[1,3]) = generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, target, algIter)
    else
        (results[1,1]. results[1,2], results[1,3]) = (0,0,0)
    end
    # Static case, no noise
    if noiseless_static == true
	println(" trying noiseless static  ###### ")
        (results[1,4], results[1,5], results[1,6], results[1,7]) = generate_and_reconstuct_static_best(model_dynamic, model_static, thetas, weights, target, algIter)
    else
        (results[1,4], results[1,5], results[1,6], results[1,7]) = (0,0,0,0)
    end
    i = 2
    # recomputes, adding differents levels of noise, for both dynamic and static.
    for noise_data in noises_data
        target = generate_target(model_dynamic,thetas,weights, noise_data, 0.0)
        # Dynamic case
        if noise_dynamic == true
	   println(" trying dynamical noise data ", noise_data, "####" )
           (results[i, 1], results[i, 2], results[i,3]) = generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, target, algIter)
	else
	   (results[i, 1], results[i, 2], results[i,3]) = (0,0,0)
	end
        # Static case
	if noise_static == true
	   println(" trying static noise data ", noise_data, "####" )
           (results[i,4], results[i,5], results[i,6], results[i,7]) = generate_and_reconstuct_static_best(model_dynamic, model_static, thetas, weights, target, algIter)
	else
	   (results[i,4], results[i,5], results[i,6], results[i,7]) = (0,0,0,0)
	end
        i = i + 1
    end
    for noise_position in noises_position
        target = generate_target(model_dynamic,thetas, weights, 0.0, noise_position)
        # This result is only analyzed in the Dynamic case for the moment
	if curvature_dynamic == true
	    println(" trying noise position", noise_position, "#####" )
            (results[i, 1], results[i, 2], results[i,3]) = generate_and_reconstruct_dynamic(model_dynamic, thetas, weights, target, algIter)
	else
	    (results[i, 1], results[i, 2], results[i,3]) = (0,0,0)
	end
        results[i, 4] = 0.0
        results[i, 5] = 0.0
        results[i, 6] = 0.0
        results[i, 7] = 0.0
        i = i + 1
    end
    println(" ################# FINISHED #######################")
    return separation, separation_dyn, results
end
end


text:
using Random;
using CPUTime;
using LinearAlgebra;
include("pep_linear.jl");
include("../expfam.jl");
include("samplingrules.jl");

# Run the learning algorithm, paramterised by a sampling rule
# The stopping and recommendation rules are common
#
# βs must be a list of thresholds *in increasing order*

function play!(i, k, rng, pep, µ, S, N, P, Vinv)
    arm = pep.arms[k]
    Y = sample(rng, getexpfam(pep, 1), arm'µ)
    S .+= Y .* arm
    Vinv .= sherman_morrison(Vinv, arm)
    N[k] += 1
end

function runit(seed, sr, μs, pep::Union{LinearBestArm,LinearThreshold}, βs, δs)
    # seed: random seed. UInt.
    # sr: sampling rule.
    # µs: mean vector.
    # pep: pure exploration problem.
    # βs: list of thresholds.
    convex_sr = typeof(sr) == ConvexGame  # test if P is needed.

    βs = collect(βs) # mutable copy

    rng = MersenneTwister(seed)

    K = narms(pep, µs)
    nb_I = nanswers(pep, µs)
    dim = length(μs)
    convex_sr ? P = ones(Int64, (nb_I, K)) : P = ones(Int64, (1, 1))  # counts detailed by answer
    N = zeros(Int64, K)              # counts
    S = zeros(dim)                     # sum of samples
    Vinv = Matrix{Float64}(I, dim, dim)  # inverse of the design matrix

    baseline = CPUtime_us()

    # pull each arm once. TODO: pull instead d arms spanning the space.
    for k = 1:K
        play!(1, k, rng, pep, µs, S, N, P, Vinv)
    end

    state = start(sr, N, P)

    R = Tuple{Int64,Array{Int64,1},Array{Int64,2},UInt64}[] # collect return values

    while true
        t = sum(N)

        hµ = Vinv * S  # emp. estimates
        #println("hµ : $(round.(hµ, digits=3)); Vinv: $Vinv")

        # test stopping criterion
        if convex_sr
            Z, _, (star, ξ) = glrt(pep, P, hμ)
        else
            Z, (_, _), (star, ξ) = glrt(pep, N, hμ)
        end
        #println("star : $star; Z : $Z; beta : $(βs[1](t))")

        while Z > βs[1](t)
        #println("Z big")
            popfirst!(βs)
            push!(R, (star, copy(N), copy(P), CPUtime_us() - baseline))
            if isempty(βs)
                return R
            end
        end

        # invoke sampling rule
        i, k = nextsample(state, pep, star, ξ, N, P, S, Vinv)
        play!(i, k, rng, pep, µ, S, N, P, Vinv)
        convex_sr ? P[i, k] += 1 : nothing
        t += 1
    end
end


text:
using miniWeather

nx_glob = 100
nz_glob = 50
sim_time = 1000.0
output_freq = 10.0
# possible options: collision, thermal, mountain_waves, turbulence, density_current, injection
#weather_type = "thermal"
weather_type = "density_current"

model, grid = init(nx_glob, nz_glob, sim_time, weather_type);

mass0, te0 = reductions(model, grid)

run!(model, grid, output_freq, "output1.nc")

mass, te = reductions(model, grid)
println("d_mass: ", (mass - mass0) / mass0)
println("d_te:   ", (te - te0) / te0)

#run!(model, grid, output_freq, "output2.nc")
#mass2, te2 = reductions(model, grid)
#println("d_mass: ", (mass2 - mass) / mass)
#println("d_te:   ", (te2 - te) / te)


text:
import DataFrames: nrow, ncol

"""
    nrow(a::Array)

Returns size of the first dimenstion of Array a (== nrow() for Arrays).
"""
function nrow(a::Array)

    return size(a)[1]
end

"""
    ncol(a::Array)

Returns size of the second dimenstion of Array a (== ncol() for Arrays).
"""
function ncol(a::Array)

    return size(a)[2]
end

"""
    rowSample(df::DataFrame, n::Int)

Take n random rows from a DataFrame with replacement.
"""
function rowSample(df::DataFrame, n::Int)

    return df[rand(1:nrow(df), n),:]
end

"""
    rowSample(df::Array, n::Int)

Take n random rows from an Array with replacement and return it as a row vector.
"""
function rowSample(a::Array, n::Int)

    return a[rand(1:nrow(a), n),:]
end

"""
    rowSample(df::Array)

Take 1 random row from an Array and return it as a 1-d-vector (== column vector).
"""
function rowSample(a::Array)

    return vec(a[rand(1:nrow(a), 1),:])
end

"""
    initCodes(num::Int, x::Array)

Return num rows from Array x as initial codes.
"""
function initCodes(num::Int, x::Array{Float64}, colNames)

    codes = rowSample(x, num)
    return codes
end

"""
    findWinner(cod, sampl)

Return index of the winner neuron for sample sampl.
"""
function findWinner(cod, sampl)

    dist = floatmax()
    winner = 1
    n = nrow(cod)

    for i in 1:n

        d = euclidean(sampl, cod[i,:])
        if (d < dist)
            dist = d
            winner = i
        end
    end

    return winner
end

"""
    normTrainData(train::DataFrame, normParams::DataFrame)

Normalise every column of training data with the params.

# Arguments
- `train`: DataFrame with training Data
- `normParams`: Shift and scale parameters for each attribute column.
"""
function normTrainData(x::Array{Float64,2}, normParams)

    for i in 1:ncol(x)
        x[:,i] = (x[:,i] .- normParams[1,i]) ./ normParams[2,i]
    end

    return x
end

"""
    normTrainData(train::DataFrame, norm::Symbol)

Normalise every column of training data.

# Arguments
- `train`: DataFrame with training Data
- `norm`: type of normalisation; one of `minmax, zscore, none`
"""
function normTrainData(train::Array{Float64,2}, norm::Symbol)

    normParams = zeros(2, ncol(train))

    if  norm == :minmax
        for i in 1:ncol(train)
            normParams[1,i] = minimum(train[:,i])
            normParams[2,i] = maximum(train[:,i]) - minimum(train[:,i])
        end
    elseif norm == :zscore
        for i in 1:ncol(train)
            normParams[1,i] = mean(train[:,i])
            normParams[2,i] = std(train[:,i])
        end
    else
        for i in 1:ncol(train)
            normParams[1,i] = 0.0  # shift
            normParams[2,i] = 1.0  # scale
        end
    end

    # do the scaling:
    if norm == :none
        x = train
    else
        x = normTrainData(train, normParams)
    end

    return x, normParams
end

function convertTrainingData(data)::Array{Float64,2}

    if typeof(data) == DataFrame
        train = convert(Matrix{Float64}, data)

    elseif typeof(data) != Matrix{Float64}
        try
            train = convert(Matrix{Float64}, data)
        catch ex
            Base.showerror(STDERR, ex, backtrace())
            error("Unable to convert training data to Array{Float64,2}!")
        end
    else
        train = data
    end

    return train
end

prettyPrintArray(arr) = println("$(show(IOContext(STDOUT, limit=true), "text/plain", arr))")


text:
include("../examples/rijke.jl")
using OrdinaryDiffEq
using PyPlot
using JLD
function plot_attractor(beta)
	n_runup = 400000
	d = N
	u0 = rand(d)
	tau = 0.2
	s = [beta, tau]
	u = Rijke_ODE(u0, s, n_runup)
	n_samples = 400000
	uf, qdot = zeros(n_samples), zeros(n_samples)
	u0 = copy(u)
	for k = 1:n_samples
		u .= Rijke_ODE(u0, s, 1)[:, end]
		uf[k] = dot(u[1:Ng], cjpixf)
		qdot[k] = beta*u[tNg]
		u0 .= u
	end
	save(string("../data/attractor/ufqdot_beta", beta, "_.jld"),
		 "uf", uf, "qdot", qdot)
	fig = figure(figsize=(10,6))
	ax = fig.add_subplot(111)
	ax.plot(qdot, uf, "b.", ms=0.2)
	ax.set_xlabel(L"\dot{q}",fontsize=46)
	ax.set_ylabel(L"u_f",fontsize=46)
	ax.xaxis.set_tick_params(labelsize=46)
	ax.yaxis.set_tick_params(labelsize=46)
	ax.set_title("\$\\beta = $(beta) \$", fontsize=46)
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
# Compute the Token Error Rate between to transcript files.
#

using ArgParse

function getargs()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--analysis", "-a"
            help = "output analysis file"
            default = ""
        "ref"
            help = "reference transcription"
            arg_type = String
            required = true
        "hyp"
            help = "hypothesis transcription"
            arg_type = String
            required = true
    end
    args = parse_args(s)
end

function loadtrans(fname::AbstractString)
    trans = Dict{String, String}()
    open(fname, "r") do f
        for line in eachline(f)
            tokens = split(line)
            trans[tokens[1]] = join(tokens[2:end], " ")
        end
    end
    trans
end

function editdistance(reftokens, hyptokens)
    distmat = zeros(Int64, length(reftokens) + 1, length(hyptokens) + 1)

    for i in 1:length(reftokens)+1 distmat[i, 1] = i-1 end
    for j in 1:length(hyptokens)+1 distmat[1, j] = j-1 end

    for i in 2:length(reftokens)+1
        for j in 2:length(hyptokens)+1

            if reftokens[i-1] == hyptokens[j-1]
                distmat[i, j] = distmat[i-1, j-1]
            else
                sub = distmat[i-1, j-1] + 1
                insert = distmat[i, j-1] + 1
                delete = distmat[i-1, j] + 1
                distmat[i, j] = min(sub, insert, delete)
            end
        end
    end
    distmat
end

function analyze(reftokens, hyptokens, editmat)
    x, y = length(reftokens) + 1, length(hyptokens) + 1

    subs = Dict{Tuple{AbstractString, AbstractString}, Int64}()
    inserts = Dict{AbstractString, Int64}()
    deletions = Dict{AbstractString, Int64}()

    while x > 1 || y > 1
        if x >= 2 && y >= 2 && editmat[x, y] == editmat[x-1, y-1] && reftokens[x-1] == hyptokens[y-1]

            x -= 1
            y -= 1
        elseif y >= 2 && editmat[x, y] == editmat[x, y-1] + 1

            key = hyptokens[y-1]
            v = get(inserts, key, 0)
            inserts[key] = v + 1

            y -= 1
        elseif x >= 2 && y >= 2 && editmat[x, y] == editmat[x-1, y-1] + 1

            key = (reftokens[x-1], hyptokens[y-1])
            v = get(subs, key, 0)
            subs[key] = v + 1

            x -= 1
            y -= 1
        else
            key = reftokens[x-1]
            v = get(deletions, key, 0)
            deletions[key] = v + 1

            x -= 1
        end
    end
    subs, inserts, deletions
end

const args = getargs()
const analysisfile = args["analysis"]
const reffile = args["ref"]
const hypfile = args["hyp"]

@info "loading the transcriptions"
const refs = loadtrans(reffile)
const hyps = loadtrans(hypfile)

const subs = Dict{Tuple{String, String}, Int64}()
const inserts = Dict{String, Int64}()
const deletions = Dict{String, Int64}()
count_utt = 0
count_token = 0
for uttid in keys(hyps)
    ref = split(refs[uttid])
    hyp = split(hyps[uttid])

    mat = editdistance(ref, hyp)
    s, i, d = analyze(ref, hyp, mat)

    for k in keys(s)
        v = get(subs, k, 0)
        subs[k] = v + 1
    end
    for k in keys(d)
        v = get(deletions, k, 0)
        deletions[k] = v + 1
    end
    for k in keys(i)
        v = get(inserts, k, 0)
        inserts[k] = v + 1
    end

    global count_utt += 1
    global count_token += length(ref)
end

# Number of substitutions
S = reduce(+, values(subs), init = 0)
D = reduce(+, values(deletions), init = 0)
I = reduce(+, values(inserts), init = 0)

# Token Error Rate
TER = round(100 * (S + D + I) / count_token, digits = 3)

@info "TER = $TER % (substitutions: $S, deletions: $D, insertions: $I) evaluated from $count_utt utterances"

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM