text:
using LogRoller, Test, Logging, Dates, JSON

rolledfile(path, n) = string(path, "_", n, ".gz")

function test_filewriter()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)
        @test !isfile(filepath)
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # initialize
        io = RollingFileWriter(filepath, 1000, 3)
        @test isfile(filepath)
        logstr = "-"^100

        # not rolled yet
        println(io, logstr)
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # roll once
        for count in 1:10
            println(io, logstr)
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed

        # roll twice
        for count in 1:10
            println(io, logstr)
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed
        @test stat(rolledfile(filepath, 2)).size > 0
        @test stat(rolledfile(filepath, 2)).size < 1000  # compressed

        # roll 4 times
        for count in 1:20
            println(io, logstr)
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test isfile(rolledfile(filepath, 2))
        @test isfile(rolledfile(filepath, 3))
        @test !isfile(rolledfile(filepath, 4)) # max 3 rolled files
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed
        @test stat(rolledfile(filepath, 2)).size > 0
        @test stat(rolledfile(filepath, 2)).size < 1000  # compressed
        @test stat(rolledfile(filepath, 3)).size > 0
        @test stat(rolledfile(filepath, 3)).size < 1000  # compressed

        close(io)
        @test !isopen(io.stream)
    end
end

function test_pipelined_tee()
    mktempdir() do logdir
        filename1 = "test1.log"
        filename2 = "test2.log"
        filepath1 = joinpath(logdir, filename1)
        filepath2 = joinpath(logdir, filename2)
        @test !isfile(filepath1)
        @test !isfile(filepath2)

        logger1io = open(filepath1, "w+")
        logger1 = SimpleLogger(logger1io)
        logger2 = RollingFileWriterTee(filepath2, 1000, 3, logger1, Logging.Info)

        julia = joinpath(Sys.BINDIR, "julia")
        cmd = pipeline(`$julia -e 'for i in 1:5 println(string("hello",i)) end; flush(stdout)'`; stdout=logger2, stderr=logger2)
        run(cmd)

        close(logger2)
        close(logger1io)

        @test isfile(filepath1)
        @test isfile(filepath2)
        @test (2*length(readlines(filepath2))) == length(readlines(filepath1)) # Julia logger entry will be two lines for every raw message
        if !Sys.iswindows()    # streams do not seem to be flushed cleanly on Windows
            @test length(readlines(filepath2)) == 5
        end
    end
end

function test_logger()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)
        @test !isfile(filepath)
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # initialize
        logger = RollingLogger(filepath, 1000, 3)
        @test isfile(filepath)
        logstr = "-"^40 # account for headers added by logger

        # not rolled yet
        with_logger(logger) do
            @info(logstr)
        end
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # roll once
        with_logger(logger) do
            for count in 1:10
                @info(logstr)
            end
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed

        # roll 4 times
        with_logger(logger) do
            for count in 1:40
                @info(logstr)
            end
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test isfile(rolledfile(filepath, 2))
        @test isfile(rolledfile(filepath, 3))
        @test !isfile(rolledfile(filepath, 4)) # max 3 rolled files
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed
        @test stat(rolledfile(filepath, 2)).size > 0
        @test stat(rolledfile(filepath, 2)).size < 1000  # compressed
        @test stat(rolledfile(filepath, 3)).size > 0
        @test stat(rolledfile(filepath, 3)).size < 1000  # compressed

        rm(rolledfile(filepath, 1))
        logger_roi = RollingLogger(filepath, 1000, 3, rotateOnInit=true)
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))

        close(logger)
        @test !isopen(logger.stream.stream)
    end
end

function test_process_streams()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)
        @test !isfile(filepath)
        @test !isfile(rolledfile(filepath, 1))

        io = RollingFileWriter(filepath, 1000, 3)
        @test isfile(filepath)

        julia = joinpath(Sys.BINDIR, "julia")
        cmd = pipeline(`$julia -e 'println("-"^100)'`; stdout=io, stderr=io)
        run(cmd)

        @test !isfile(rolledfile(filepath, 1))

        for count in 1:10
            run(cmd)
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))

        @test io.procstream !== nothing
        @test io.procstreamer !== nothing
        @test !istaskdone(io.procstreamer)

        close(io)
        @test io.procstream === nothing
        @test io.procstreamer === nothing
    end
end

function test_postrotate()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)
        @test !isfile(filepath)
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # initialize
        logger = RollingLogger(filepath, 1000, 3)
        rotatedfiles = Vector{String}()
        postrotate(logger) do rotatedfilename
            push!(rotatedfiles, rotatedfilename)
        end
        @test isfile(filepath)
        logstr = "-"^40 # account for headers added by logger

        # not rolled yet
        with_logger(logger) do
            @info(logstr)
        end
        @test !isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))

        # roll once
        with_logger(logger) do
            for count in 1:10
                @info(logstr)
            end
        end
        @test isfile(filepath)
        @test isfile(rolledfile(filepath, 1))
        @test !isfile(rolledfile(filepath, 2))
        @test !isfile(rolledfile(filepath, 3))
        @test stat(filepath).size > 0
        @test stat(filepath).size < 1000
        @test stat(rolledfile(filepath, 1)).size > 0
        @test stat(rolledfile(filepath, 1)).size < 1000  # compressed

        @test length(rotatedfiles) == 1
        @test rotatedfiles[1] == rolledfile(filepath, 1)
    end
end

function test_timestamp_handling()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)

        logger = RollingLogger(filepath, 1000, 3)
        testdt = DateTime(2000, 1, 1, 1, 1, 1, 1)
        timenow = time()
        dt_from_time = Dates.unix2datetime(timenow)
        with_logger(logger) do
            @info("without_timestamp")
            @info("without_timestamp", testarg=1)
            @info("with_timestamp", time=timenow)
            @info("with_timestamp", time=testdt)
            @info("with_timestamp", time=testdt, testarg=1)
            @test isfile(filepath)
        end
        close(logger)

        filecontents = readlines(filepath)
        @test length(filecontents) == 12
        dates = DateTime[]
        for line in filecontents
            if startswith(line, "┌ Info: ")
                parts = split(line, ' ')
                @test length(parts) == 4
                push!(dates, DateTime(strip(parts[3])[1:(end-1)]))
            end
        end
        @test length(dates) == 5
        @test dates[3] == dt_from_time
        @test dates[4] == testdt
        @test dates[5] == testdt
    end
end

function test_json_format()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)

        logger = RollingLogger(filepath, 2000, 3; format=:json)
        testdt = DateTime(2000, 1, 1, 1, 1, 1, 1)
        timenow = time()
        dt_from_time = Dates.unix2datetime(timenow)
        with_logger(logger) do
            @info("without_timestamp")
            @info("without_timestamp", testarg=1)
            @info("with_timestamp", time=timenow)
            @info("with_timestamp", time=testdt)
            @info("with_timestamp", time=testdt, testarg=1)
            @test isfile(filepath)
        end
        close(logger)

        dates = DateTime[]
        open(filepath) do readio
            while !eof(readio)
                try
                    entry = JSON.parse(readio)
                    push!(dates, DateTime(entry["metadata"]["time"]))
                catch ex
                    eof(readio) || rethrow()
                end
            end
        end

        @test length(dates) == 5
        @test dates[3] == dt_from_time
        @test dates[4] == testdt
        @test dates[5] == testdt
    end

    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)

        logger = RollingLogger(filepath, 20000, 3; format=:json)
        with_logger(logger) do
            for i in 1:4
                @info("log message $i", time=time(), randval=rand(Int))
            end
            with_logger(current_logger()) do
                @info("test nested with_logger", time=time())
            end
            @info rand(1000, 1000)
            @info(Vector{Bool})
            try
                error("test exception")
            catch ex
                @error("caught an exception", ex)
                @error("this is the exception with backtrace", exception=(ex,catch_backtrace()))
            end
            @warn("test other types", ptr1=Ptr{Nothing}(), ptr2=Ptr{Int}(10), sv=Core.svec(1,2,3), typ=String)
        end
        close(logger)

        open(filepath) do readio
            for idx in 1:4
                entry = JSON.parse(readio)
                @test entry["metadata"]["level"] == "Info"
                @test entry["message"] == "log message $idx"
            end
            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Info"
            @test entry["message"] == "test nested with_logger"

            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Info"
            @test length(entry["message"]) == (4*1024 + 3)

            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Info"
            @test entry["message"] == "Array{Bool,1}"

            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Error"
            @test entry["message"] == "caught an exception"
            @test haskey(entry["keywords"], "ex")
            @test haskey(entry["keywords"]["ex"], "msg")
            @test entry["keywords"]["ex"]["msg"] == "test exception"

            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Error"
            @test entry["message"] == "this is the exception with backtrace"
            @test haskey(entry["keywords"], "exception")
            @test startswith(entry["keywords"]["exception"], "test exception\nStacktrace:")

            entry = JSON.parse(readio)
            @test entry["metadata"]["level"] == "Warn"
            @test entry["message"] == "test other types"
            @test haskey(entry["keywords"], "ptr1")
            @test haskey(entry["keywords"], "ptr2")
            @test haskey(entry["keywords"], "sv")
            @test haskey(entry["keywords"], "typ")
            @test entry["keywords"]["sv"] == [1,2,3]
            @test entry["keywords"]["typ"] == "String"
        end
    end
end

function test_size_limits()
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)

        logger = RollingLogger(filepath, 2000, 3; format=:json, entry_size_limit=250+length(filepath))
        with_logger(logger) do
            @info("short msg")
            @info("short msg", v=1)
            @info("long msg", v="-"^100)
            @info("long msg " * "-"^100)
        end
        close(logger)

        open(filepath) do readio
            entry = JSON.parse(readio)
            @test entry["message"] == "short msg"
            @test !haskey(entry, "keywords")
            entry = JSON.parse(readio)
            @test entry["message"] == "short msg"
            @test entry["keywords"]["v"] == 1
            entry = JSON.parse(readio)
            @test entry["message"] == "long msg"
            @test !haskey(entry, "keywords")
            entry = JSON.parse(readio)
            @test !haskey(entry, "message")
            @test !haskey(entry, "keywords")
        end
    end
    mktempdir() do logdir
        filename = "test.log"
        filepath = joinpath(logdir, filename)

        logger = RollingLogger(filepath, 2000, 3; format=:console, entry_size_limit=100)
        with_logger(logger) do
            @info("short msg")
            @info("short msg", v=1)
            @info("long msg", v="-"^100)
            @info("long msg " * "-"^100)
        end
        close(logger)
        filecontents = readlines(filepath)
        @test length(filecontents) == 8
    end
end

@testset "file writer" begin
    test_filewriter()
end
@testset "pipelined tee" begin
    test_pipelined_tee()
end
@testset "process streams" begin
    test_process_streams()
end
@testset "logger" begin
    test_logger()
    test_timestamp_handling()
    test_postrotate()
    test_json_format()
    test_size_limits()
end


text:
#####
##### Outer functions for setting velocity on boundary and filling halo beyond boundary.
#####

@kernel function set_east_west_u_velocity!(u, i_boundary, bc, grid, clock, state)
    j, k = @index(Global, NTuple)
    @inbounds u[i_boundary, j, k] = getbc(bc, j, k, grid, clock, state)
end

@kernel function set_north_south_v_velocity!(v, j_boundary, bc, grid, clock, state)
    i, k = @index(Global, NTuple)
    @inbounds v[i, j_boundary, k] = getbc(bc, i, k, grid, clock, state)
end

@kernel function set_top_bottom_w_velocity!(w, k_boundary, bc, grid, clock, state)
    i, j = @index(Global, NTuple)
    @inbounds w[i, j, k_boundary] = getbc(bc, i, j, grid, clock, state)
end

  fill_west_halo!(u, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :yz, set_east_west_u_velocity!,   u,           1, bc, grid, clock, state, dependencies=dep)
  fill_east_halo!(u, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :yz, set_east_west_u_velocity!,   u, grid.Nx + 1, bc, grid, clock, state, dependencies=dep)
 fill_south_halo!(v, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :xz, set_north_south_v_velocity!, v,           1, bc, grid, clock, state, dependencies=dep)
 fill_north_halo!(v, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :xz, set_north_south_v_velocity!, v, grid.Ny + 1, bc, grid, clock, state, dependencies=dep)
fill_bottom_halo!(w, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :xy, set_top_bottom_w_velocity!,  w,           1, bc, grid, clock, state, dependencies=dep)
   fill_top_halo!(w, bc::NFBC, arch, dep, grid, clock, state) = launch!(arch, grid, :xy, set_top_bottom_w_velocity!,  w, grid.Nz + 1, bc, grid, clock, state, dependencies=dep)

ext: jl


text:
# Supported domains: Real, Boolean, Ordinal, Periodic, Count

# The purpose of domains is to be able to sample over different possible values of `a` regardless of
# the loss that was used in the GLRM. The reason for doing this is to evaluate the performance of GLRMS.
# For instance, let's say we use PCA (QuadLoss losses) to model a binary data frame (not the best idea).
# In order to override the standard imputation with `sample(QuadLoss(), u)`, which assumes imputation over the reals,
# we can use `sample(BoolDomain(), QuadLoss(), u)` and see which of {-1,1} is best. The reason we want to be able to
# do this is to compare a baseline model (e.g. PCA) with a more logical model using heterogenous losses,
# yet still give each model the same amount of information regarding how imputation should be done.

# The domains themselves are defined in domains.jl

# In order to accomplish this we define a series of domains that describe how imputation should be performed over
# them. Each combination of domain and loss must have the following:
#	Methods:
#     `sample(D::my_Domain, l::my_loss_type, u::Float64) ::Float64`
#           Samples aᵤ from among the range of possible values of a. The range of
#			possible values of a should be implicitly or explicitly provided by `D`.
#			There should be an sample method for every combination of datatype and loss.

# DataTypes are assigned to each column of the data and are not part of the low-rank model itself, they just serve
# as a way to evaluate the performance of the low-rank model.

import StatsBase: sample, WeightVec
export sample, sample_missing

########################################## REALS ##########################################
# Real data can take values from ℜ

# l.scale should be 1/var
sample(D::RealDomain, l::QuadLoss, u::Float64) = u + randn()/sqrt(l.scale)

########################################## BOOLS ##########################################
# Boolean data should take values from {true, false}

function sample(D::BoolDomain, l::LogisticLoss, u::Float64)
	rand()<=(1/(1+exp(-u))) ? true : false
end

# generic method
# Evaluate w/ a=-1 and a=1 and see which is better according to that loss.
# This is fast and works for any loss.
function sample(D::BoolDomain, l::Loss, u::AbstractArray)
    prob = exp.(-[evaluate(l, u, i) for i in (true, false)])
		return sample(WeightVec(prob))
end

########################################## ORDINALS ##########################################
# Ordinal data should take integer values ranging from `min` to `max`

# a DiffLoss is one in which l(u,a) = f(u-a) AND argmin f(x) = 0
# for example, QuadLoss(u,a)=(u-a)² and we can write f(x)=x² and x=u-a
function sample(D::OrdinalDomain, l::DiffLoss, u::Float64)
		uint = round(Int, u)
		uclip = max(D.min, min(D.max, uint))
		return uclip
end

# generic method
function sample(D::OrdinalDomain, l::Loss, u::AbstractArray)
    prob = exp.(-[evaluate(l, u, i) for i in D.min:D.max])
		return sample(WeightVec(prob))
end

########################################## CATEGORICALS ##########################################
# Categorical data should take integer values ranging from 1 to `max`

function sample(D::CategoricalDomain, l::MultinomialLoss, u::Array{Float64})
	return sample(WeightVec(exp.(u)))
end

# sample(D::CategoricalDomain, l::OvALoss, u::Array{Float64}) = ??

# generic method
function sample(D::CategoricalDomain, l::Loss, u::AbstractArray)
    prob = exp.(-[evaluate(l, u, i) for i in D.min:D.max])
		return sample(WeightVec(prob))
end

########################################## PERIODIC ##########################################
# Periodic data can take values from ℜ, but given a period T, we should have error_metric(a,a+T) = 0

# Since periodic data can take any real value, we can use the real-valued imputation methods
# sample(D::PeriodicDomain, l::Loss, u::Float64) = ??

########################################## COUNTS ##########################################
# Count data can take values over ℕ, which we approximate as {0, 1, 2 ... `max_count`}

# Our approximation of ℕ is really an ordinal
sample(D::CountDomain, l::Loss, u::Float64) = sample(OrdinalDomain(0,D.max_count), l, u)

####################################################################################
# Use impute and error_metric over arrays
function sample{DomainSubtype<:Domain,LossSubtype<:Loss}(
			domains::Array{DomainSubtype,1},
			losses::Array{LossSubtype,1},
			U::Array{Float64,2}) # U = X'*Y
	m, d = size(U)
	n = length(losses)
	yidxs = get_yidxs(losses)
	A_sampled = Array(Number, (m, n));
	for f in 1:n
		for i in 1:m
			if length(yidxs[f]) > 1
				A_sampled[i,f] = sample(domains[f], losses[f], vec(U[i,yidxs[f]]))
			else
				A_sampled[i,f] = sample(domains[f], losses[f], U[i,yidxs[f]])
			end
		end
	end
	return A_sampled
end

# sample missing entries in A according to the fit model (X,Y)
function sample_missing(glrm::GLRM)
		do_sample(e::Int, f::Int) = !(e in glrm.observed_examples[f])
		return sample(glrm, do_sample)
end

all_entries(e::Int,f::Int) = true

# sample all entries in A according to the fit model (X,Y)
# do_sample is a function that takes an example-feature pair (e,f)
# and returns true if that entry should be replaced by a sample from the model
# is_dense controls whether the output should be a dense matrix
# it's true by default because we sample all entries by default
function sample(glrm::GLRM, do_sample::Function=all_entries, is_dense::Bool=true)
	U = glrm.X'*glrm.Y
	m, d = size(U)
	n = length(glrm.losses)
	yidxs = get_yidxs(glrm.losses)
	domains = Domain[domain(l) for l in glrm.losses]
	# make sure we don't mutate the type of the array A
	# even if all data for some real loss take integer values
	for j=1:n
		if isa(domains[j], RealDomain) && isa(glrm.A[j], DataArray{Int64,1})
			domains[j] = OrdinalDomain(minimum(dropna(glrm.A[j])), maximum(dropna(glrm.A[j])))
		end
	end
	A_sampled = copy(glrm.A);
	if is_dense && isa(A_sampled, SparseMatrixCSC)
		A_sampled = Matrix(A_sampled)
	end
	for f in 1:n
		for e in 1:m
			if do_sample(e,f)
				A_sampled[e,f] = sample(domains[f], glrm.losses[f], U[e,yidxs[f]])
			end
		end
	end
	return A_sampled
end

function sample{LossSubtype<:Loss}(losses::Array{LossSubtype,1}, U::Array{Float64,2})
	domains = Domain[domain(l) for l in losses]
	sample(domains, losses, U)
end

### Hack to sample from non-probabilistic losses
sample(D::Domain, l::Loss, u) = impute(D, l, u)


text:
using SigmoidNumbers
using Base.Test

#/((ULP{4,0}(0b0001) → ULP{4,0}(0b0001)), (Exact{4,0}(0b0000) → Exact{4,0}(0b1000))) |> println
#exit()

#include("h-layer-test.jl")
#include("breakdown-test.jl")
#
#include("Posits5-test.jl")
#include("Posits7-test.jl")
#
#include("misc-test.jl")
#include("fdp-test.jl")
#
##include("fdp-matrix-test.jl")
#
##tests for the valids
#include("warlpiri-test.jl")

include("blas-test.jl")

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

"""
    YAActL

Yet Another Actor Library (in Julia `VERSION ≥ v"1.3"`).

The current stable, registered version is installed with
```julia
pkg> add YAActL
```

The development version is installed with:
```julia
pkg> add("https://github.com/juliaactors/YAActL.jl")
```
"""
module YAActL

"Gives the package version."
const version = v"0.2.2"

using Distributed

include("types.jl")
include("messages.jl")
include("links.jl")
include("com.jl")
include("actors.jl")
include("api.jl")
include("registry.jl")
include("diag.jl")
include("init.jl")

export  Message, Response, Request, Timeout, Func, Args,
        Link, LinkParams, parallel, USR, islocal,
         send!, request!, receive!,
        Actor, become, stop,
        Dispatch, full, state,
        become!, call!, cast!, exec!, exit!, init!,
        query!, self, set!, term!, update!,
        register, unregister, whereis, registered,
        register!, info

end


text:
module ModelParametersTest

using Base.Test
using ModelParameters

const x = parameter(name = "A", id = :a, value = 10)

@test value(x) == 10
@test name(x) == "A"
@test id(x) === :a
@test units(x) == "-"
@test description(x) == ""
@test transformer(x) === identity

@test transform(x) == 10

const mass = parameter(
	name = "Mass",
	id = :mass,
	value = 300.0,
	units = "g",
	transformer = x -> x / 1000.0)

@test transform(mass) ≈ 0.3

const children = [
	parameter(
		name = "Temperature",
		id = :temperature,
		value = 30.0,
		units = "°C",
		transformer = x -> x + 273.15),
	parameter(
		name = "Pressure",
		id = :pressure,
		value = 2.0,
		units = "bar",
		transformer = x -> x * 1e5)]

const idealgas = parameter(
	name = "Ideal gas",
	id = :idealgas,
	children = children)

@test id(idealgas) === :idealgas
@test value(idealgas) == Dict(:temperature => 30.0, :pressure => 2.0)

setvalue!(idealgas, Dict(:pressure => 1.5))

@test value(get(idealgas, :pressure)) == 1.5

let expected = Dict(:pressure => 1.5 * 1e5, :temperature => 30.0 + 273.15)
	for (k,v) in transform(idealgas)
		@test_approx_eq v expected[k]
	end
end

const temps = parameter(
	name = "Temperature",
	id = :temperature,
	children = [
		parameter(
			name = "T (°C)",
			id = :celcius,
			value = 30.0,
			units = "°C",
			transformer = x -> x + 273.15),
		parameter(
			name = "T (K)",
			id = :kelvin,
			value = 293.15,
			units = "K")],
	selection = :celcius)

@test typeof(temps) == ParameterOptions

@test value(temps) == 30.0

# transform(temps)

@test transform(temps) ≈ 303.15

@test_throws KeyError setselection!(temps, :blah)

setselection!(temps, :kelvin)

@test value(temps) ≈ 293.15

@test transform(temps) ≈ 293.15

const velocity = parameter(
	name = "Velocity",
	id = :velocity,
	value = 3.0,
	units = "m/s",
	perturbation = 0.1)

@test value(velocity) == 3.0

let err = 0.0, larger = 0, less = 0
	for i=1:10000
		let v = transform(velocity)
			if v > 3.0
				larger += 1
			else
				less += 1
			end
			err += abs(transform(velocity) - 3.0)
		end
	end
	@test 1450.0 < err < 1650.0
	@test 4500 < less < 5500
	@test 4500 < larger < 5500
end

const vols = parameter(
	name = "Volumes",
	id = :volumes,
	size = 1000,
	children = [
		parameter(
			name = "Volume",
			id = :volume,
			value = 20.0,
			units = "L",
			transformer = x -> x / 1000.0,
			perturbation = 0.1)],
	transformer = (; args...) -> Dict{Symbol, Any}(args))

@test value(vols) == Dict(:broadcastsize => 1000, :volume => 20.0)

@test typeof(transform(vols)) == Vector{Dict{Symbol, Any}}

let v = transform(vols), err = 0.0, larger = 0, n = broadcastsize(vols)
	for i=1:n
		if v[i][:volume] > 0.02
			larger += 1
		end
		err += abs(v[i][:volume] - 0.02)
	end
	err /= 0.02   # relative error
	err /= n      # average relative error
	@test 0.045 < err < 0.055
	@test (n/2 * 0.95) < larger < (n/2 * 1.05)
end

end # module


text:
using Documenter
using Plots
using SpatialDatasets
using SpatialDependence
using StableRNGs

# Workaround for GR errors when plotting
ENV["GKSwstype"] = "100"

# Load datasets to not display Downloading Artifact message in the documentation
begin
    guerry = sdataset("Guerry")
    boston = sdataset("Bostonhsg")
end

DocMeta.setdocmeta!(SpatialDependence, :DocTestSetup, :(using SpatialDependence); recursive=true)

makedocs(;
    modules=[SpatialDependence],
    authors="Javier Barbero and contributors",
    repo="https://github.com/javierbarbero/SpatialDependence.jl/blob/{commit}{path}#{line}",
    sitename="SpatialDependence.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://javierbarbero.github.io/SpatialDependence.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Spatial Weight Matrices" => "sweightmatrices.md",
        "Spatial Lag" => "spatiallag.md",
        "Global Spatial Autocorrelation" => "gspatialautocorrelation.md",
        "Local Spatial Autocorrelation" => "lspatialautocorrelation.md",
        "Choropleth Maps" => "choropleth.md",
        "Parallel Computing" => "parallelcomputing.md",
        "Bibliography" => "bibliography.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/javierbarbero/SpatialDependence.jl",
    devbranch = "main"
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ThriftJuliaCompiler_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ThriftJuliaCompiler")
JLLWrappers.@generate_main_file("ThriftJuliaCompiler", UUID("815b9798-8dd0-5549-95cc-3cf7d01bce66"))
end  # module ThriftJuliaCompiler_jll


text:
module ContinuousInference

using LPSjl

"""
```julia
initial_net(bayesianNetwork::BayesianNetwork, A::Array{String}, a::Array{Int})
```
# Arguments
- `bayesianNetwork`: A singly-connected BN(G, P) containing continuous variables, where G=(V, E)
- `A`: a subset Α ⊆ V
- `a`: a set of values for A

# Returns
- The BN(G, P) updated according to the values in Α. All expectations and variances, including those in messages, are considered part of the network.

# References
- Neapolitan, p. 187, Algorithm 4.1 - Inference With Continuous Variables
"""
function initial_net(bayesianNetwork::BayesianNetwork, A::Array{String}, a::Array{Int})

end

function update_tree(bayesianNetwork::BayesianNetwork, A::Array{String}, a::Array{String},
                    V::String, v::Int)

end

function send_lambda_message(Y::String, X::String)

end

function send_pi_message(Z::String, X::String)

end

export initial_net

end # module


text:
function combinationSum2(cands, target)
    sort!(cands) # sort all cands
    listOfSets = [Set() for _ in 1:(target+1)]

    push!(listOfSets[1], []) # index 1 actually point to 0, because julia index start from 1

    for
        cand in cands,
        t in range(target, cand, step=-1),
        prev in listOfSets[t+1-cand];

        cache = copy(prev)
        push!(cache, cand)
        push!(listOfSets[t+1], cache)
    end

    last(listOfSets)
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
using .Kronecker

export kernelkronmat

function kernelkronmat(
    κ::Kernel,
    X::AbstractVector,
    dims::Int
    )
    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices (see `iskroncompatible()`)"
    k = kernelmatrix(κ, X)
    kronecker(k, dims)
end

function kernelkronmat(
    κ::Kernel,
    X::AbstractVector{<:AbstractVector};
    obsdim::Int=defaultobs
    )
    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices"
    Ks = kernelmatrix.(κ, X)
    K = reduce(⊗, Ks)
end

"""
    To be compatible with kroenecker constructions the kernel must satisfy
    the property : for x,x' ∈ ℜᴰ
    k(x,x') = ∏ᵢᴰ k(xᵢ,x'ᵢ)
"""
@inline iskroncompatible(κ::Kernel) = false # Default return for kernels


text:
# To be removed after running this test individually in the REPL successfully
@everywhere using DSGE
@everywhere using HDF5, JLD2, Random, DelimitedFiles
@everywhere import Test: @test, @testset

write_test_output = false

path = dirname(@__FILE__)

###################################################################
# An Schorfheide
###################################################################
m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = h5read("$path/../../reference/smc.h5", "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, false)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:tempering_target, 0.95)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:use_fixed_schedule, true)

####################################################################
test_init_cloud = ParticleCloud(m, get_setting(m, :n_particles))
@everywhere Random.seed!(42)
DSGE.initial_draw!(m, data, test_init_cloud)

if write_test_output
    #=JLD.jldopen("$path/../../reference/initial_draw.jld", "w") do file
    write(file, "cloud", test_init_cloud)
    end =#
    JLD2.jldopen("$path/../../reference/initial_draw.jld2", "w") do file
        write(file, "cloud", test_init_cloud)
    end
end

saved_init_cloud = load("$path/../../reference/initial_draw.jld2", "cloud")

@testset "Initial draw" begin
    @test @test_matrix_approx_eq DSGE.get_vals(test_init_cloud) DSGE.get_vals(saved_init_cloud)
    @test @test_matrix_approx_eq DSGE.get_loglh(test_init_cloud) DSGE.get_loglh(saved_init_cloud)
end

#=###################################################################
# Smets Wouters
###################################################################
m = SmetsWoutersOrig()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = readdlm("$path/../../reference/YY.txt")
data = Matrix{Float64}(data')

m <= Setting(:n_particles, 12000)
m <= Setting(:n_Φ, 500)
m <= Setting(:λ, 2.1)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, false)
m <= Setting(:step_size_smc, 0.4)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, 0.9)
m <= Setting(:tempering_target, 0.95)
m <= Setting(:resampling_threshold, 0.5)
m <= Setting(:use_fixed_schedule, true)
m <= Setting(:use_chand_recursion, true)

####################################################################
test_init_cloud = ParticleCloud(m, get_setting(m, :n_particles))

@everywhere Random.seed!(42)
DSGE.initial_draw!(m, data, test_init_cloud)

if write_test_output
    JLD2.jldopen("$path/../../reference/initial_draw_sw.jld2", true, true, true, IOStream) do file
        write(file, "cloud", test_init_cloud)
    end
end

saved_init_cloud = load("$path/../../reference/initial_draw_sw.jld2", "cloud")

@testset "Initial draw" begin
    @test @test_matrix_approx_eq DSGE.get_vals(test_init_cloud) DSGE.get_vals(saved_init_cloud)
    @test @test_matrix_approx_eq DSGE.get_loglh(test_init_cloud) DSGE.get_loglh(saved_init_cloud)
end
=#


text:
module TestHessian
using Test
using CovarianceFunctions
using CovarianceFunctions: EQ, RQ, Dot, ExponentialDot, NN
using CovarianceFunctions: HessianKernel, ValueGradientHessianKernel,
                           input_trait, BlockFactorization,
                           IsotropicHessianKernelElement,
                           DotProductHessianKernelElement,
                           ValueGradientHessianKernelElement
using LinearAlgebra
const AbstractMatOrFac = Union{AbstractMatrix, Factorization}

@testset "Hessian kernels" begin
    # nd test
    n = 2
    d = 3
    X = randn(d, n) / sqrt(d)
    ε = 1e2eps()
    @testset "HessianKernelElements" begin
        kernels = [EQ(), RQ(1.), Dot()^4, ExponentialDot()] #, NN()]
        for k in kernels
            K = HessianKernel(k)(X[:, 1], X[:, 2])
            @test K isa Union{IsotropicHessianKernelElement, DotProductHessianKernelElement}
            H = HessianKernel((x, y) -> k(x, y))(X[:, 1], X[:, 2])
            @test H isa Matrix{<:Real}
            @test size(K) == (d^2, d^2)
            @test Matrix(K) ≈ H
            a, b = randn(d^2), randn(d^2)
            α, β = randn(2)
            c = α * H * a + β * b
            mul!(b, K, a, α, β)
            @test b ≈ c
        end
    end

    @testset "HessianKernel" begin
        aa = randn(n*d^2)
        bb = randn(n*d^2)
        kernels = [EQ(), RQ(1.), Dot()^3, Dot()^4] #, NN()]
        for k in kernels
            G = HessianKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            @test K isa BlockFactorization{<:Real}

            MK = Matrix(K)
            @test maximum(abs, MK - MK') < ε
            MK = Symmetric(MK)
            @test all(≥(-10maximum(MK)*eps()), eigvals(MK)) # positive semidefinite
            @test size(MK) == (d^2*n, d^2*n)

            # create anonymous wrapper around kernel to trigger generic fallback
            k2 = (x, y)->k(x, y)
            G2 = HessianKernel(k2)
            # G2(x[1], x[2]) # this takes centuries to precompile, maybe because of nested differentiation?
            # not really important because it is a fallback anyway
            K2 = CovarianceFunctions.gramian(G2, X)
            MK2 = Matrix(K2)
            @test MK ≈ MK2 # compare generic and specialized evaluation of gradient kernel

            # testing matrix multiply
            Kab = deepcopy(bb)
            α, β = randn(2) # trying out two-argument mul
            mul!(Kab, K, aa, α, β) # saves some memory, Woodbury still allocates a lot
            MKab = @. α * $*(MK, aa) + β * bb
            @test Kab ≈ MKab
        end
    end

    dd = d^2 + d + 1
    @testset "ValueGradientHessianKernelElement" begin
        kernels = [EQ(), RQ(1.)] #, Dot()^2, Dot()^4]# Dot()^5] #, NN()]
        for k in kernels
            K = ValueGradientHessianKernel(k)(X[:, 1], X[:, 2])
            @test K isa ValueGradientHessianKernelElement
            H = ValueGradientHessianKernel((x, y) -> k(x, y))(X[:, 1], X[:, 2])
            @test H isa Matrix{<:Real}
            @test size(K) == (dd, dd)
            @test Matrix(K) ≈ H

            # block indices
            vi = 1 # value index
            gi = (1+1):(d+1) # gradient index
            hi = (d+1+1):dd # hessian index

            a = randn(dd)
            # a[vi] = 0 # used to debug individual components
            # @. a[gi] = 0
            # @. a[hi] = 0
            b = randn(dd)
            α, β = randn(2)
            c = α * H * a + β * b
            mul!(b, K, a, α, β)
            @test b[vi] ≈ c[vi]
            @test b[gi] ≈ c[gi]
            @test b[hi] ≈ c[hi]
        end
    end

    @testset "ValueGradientHessianKernel" begin
        kernels = [EQ()] #, RQ(1.), Dot(), Dot()^3, NN()]
        ε = 1e4eps()
        aa = randn(n*dd)
        bb = randn(n*dd)
        for k in kernels
            G = ValueGradientHessianKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            @test K isa BlockFactorization{<:Real}

            MK = Matrix(K)
            @test maximum(abs, MK - MK') < ε
            MK = Symmetric(MK)
            @test all(≥(-1e4eps()), eigvals(MK)) # positive semidefinite
            @test size(MK) == (dd*n, dd*n)
            # # create anonymous wrapper around kernel to trigger generic fallback
            # k2 = (x, y)->k(x, y)
            # G2 = ValueGradientKernel(k2)
            # # G2(x[1], x[2]) # this takes centuries to precompile, maybe because of nested differentiation?
            # # not really important because it is a fallback anyway
            # K2 = CovarianceFunctions.gramian(G2, X)
            # MK2 = Matrix(K2)
            # @test MK ≈ MK2 # compare generic and specialized evaluation of gradient kernel
            #
            # testing matrix multiply
            Kab = deepcopy(bb)
            α, β = randn(2) # trying out two-argument mul
            mul!(Kab, K, aa, α, β) # saves some memory, Woodbury still allocates a lot
            MKab = @. α * $*(MK, aa) + β * bb
            @test Kab ≈ MKab
        end
        # testing matrix solve
        for k in [EQ(), RQ(1.)]
            G = ValueGradientHessianKernel(k)
            K = CovarianceFunctions.gramian(G, X)
            a = randn(n*dd)
            Ka = K*a
            as = K\Ka
            @test norm(K*as-Ka) < 1e-6
        end
    end
end

end # module


text:
using Printf
using StaticArrays

export Atom

mutable struct Atom
    atno::Int
    xyz::MVector{3,Float64}
end

atom(sym::String,xyz::MVector{3,Float64},units=:Angstrom) = atom(sym2no[sym],xyz,units)
atom(atno::Int,x::Float64,y::Float64,z::Float64,units=:Angstrom) = atom(atno,MVector(x,y,z),units)
atom(sym::String,x::Float64,y::Float64,z::Float64,units=:Angstrom)= atom(sym2no[sym],MVector(x,y,z),units)
function atom(atno::Int,xyz::MVector{3,Float64},units=:Angstrom)
    if units == :Angstrom xyz /= 0.52918 end
    return Atom(atno,xyz)
end

nuclear_repulsion(a::Atom,b::Atom)= a.atno*b.atno/sqrt(dist2(a.xyz,b.xyz))
nuclear_repulsion(atoms::Vector{Atom}) = sum(nuclear_repulsion(atoms[i],atoms[j]) for (i,j) in spairs(nat(atoms)))

nel(atoms::Vector{Atom}) = sum([at.atno for at in atoms])
nat(atoms::Vector{Atom}) = length(atoms)

# Array of symbols, masses
symbol = [
    "H","He",
    "Li","Be","B","C","N","O","F","Ne",
    "Na","Mg","Al","Si","P","S","Cl","Ar",
    "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe",
    "Co", "Ni", "Cu", "Zn",
    "Ga", "Ge", "As", "Se", "Br", "Kr",
    "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru",
    "Rh", "Pd", "Ag", "Cd",
    "In", "Sn", "Sb", "Te", "I", "Xe",
    "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm",  "Eu",
    "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu",
    "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
    "Tl","Pb","Bi","Po","At","Rn"];
sym2no = Dict((sym => i) for (i,sym) in enumerate(symbol));
mass = [
    1.0008, 4.0026,
    6.941,9.0122,
    10.811,12.011,14.007,15.999,18.998,20.179,
    22.990,24.305,
    26.982,28.086,30.974,32.066,35.453,39.948,
    39.098, 40.078,
    44.9559, 47.867, 50.9415, 51.9961, 54.938, 55.845,
    58.9332, 58.6934, 63.546,65.39,
    69.723, 72.61, 74.9216, 78.96, 79.904, 83.80,
    85.4678, 87.62,
    88.90686, 91.224, 92.90638, 95.94, 98, 101.07,
    102.90550, 106.42, 107.8682, 112.411,
    114.818, 118.710, 121.760, 127.60, 126.90447, 131.29,
    132.90545, 137.327, 138.9055, 140.11, 140.90765, 144.24,
    145.0, 150.36, 151.964,
    157.25, 158.92534, 162.5, 164.93, 167.259, 168.934, 173.04, 174.967,
    178.49, 180.9479, 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655,
    200.59];

    function read_xyz(lines)
        nat = parse(Int,strip(lines[1]))
        title = strip(lines[2])
        atoms = Atom[]
        for i in 1:nat
            words = split(lines[2+i])
            sym = words[1]
            atno = sym2no[sym]
            xyz = MVector(parse(Float64,w) for w in words[2:4])
            push!(atoms,Atom(atno,xyz))
        end
        return atoms
    end
    function write_xyz(atoms,fname,title="Written by MolecularIntegrals.jl")
        f = open(fname,"w")
        write(f,"$(length(atoms))\n")
        write(f,"$title\n")
        for atom in atoms
            @printf(f,"%-10s %10.6f %10.6f %10.6f\n",
                symbol[atom.atno],atom.xyz[1],atom.xyz[2],atom.xyz[3])
        end
        close(f)
    end
    calcmass(atoms) = sum(mass[atom.atno] for atom in atoms)
    com(atoms) = sum(mass[atom.atno]*[atom.xyz] for atom in atoms)/calcmass(atoms)
    function translate!(atoms,xyz)
        for atom in atoms
            atom.xyz += xyz
        end
        return nothing
    end
    # using Bio3DView, able to view("mol.xyz") for a three.js view of molecule.


text:
"""
# PiAVRDude.jl

AVRDude interface for Raspberry Pi GPIO ISP.
"""
module PiAVRDude

export AVRDude

using UnixIO

struct AVRDude

    name::String
    device::String
    config::String
    usb_port::Union{Nothing, String}

    function AVRDude(;device=nothing,
                      sck=nothing, miso=nothing, mosi=nothing, reset=nothing,
                      usb_port=nothing)
        if usb_port != nothing
            new("avrisp2", device, """
                default_safemode = no;
                """,
                usb_port)
        else
            name = "gpio$reset"
            new(name, device, """
                default_safemode = no;
                programmer
                  id    = "$name";
                  desc  = "GPIO ISP";
                  type  = "linuxgpio";
                  sck   = $sck;
                  miso  = $miso;
                  mosi  = $mosi;
                  reset = $reset;
                ;
                """,
                nothing)
        end
    end
end

config(isp) = isp.config

function avrdude(isp, cmd)
    mktempdir() do d
        conf=joinpath(d, "avrdude.conf")
        write(conf, config(isp))
        options = `-p $(isp.device) -C +$conf -c $(isp.name)`
        if isp.usb_port != nothing
            options = `$options -D -P $(isp.usb_port)`
        end
        cmd = `avrdude $options $cmd`
        @info cmd
        UnixIO.system(cmd)
    end
end

status(isp) = avrdude(isp, `-v`)

flash_hex(isp, file) = avrdude(isp, `-U flash:w:$file:i`)
flash_elf(isp, file) = avrdude(isp, `-U flash:w:$file:e`)
eeprom_hex(isp, file) = avrdude(isp, `-U eeprom:w:$file:i`)
eeprom_elf(isp, file) = avrdude(isp, `-U eeprom:w:$file:e`)
flash(i, f) = endswith(f, ".elf") ? flash_elf(i, f) : flash_hex(i, f)
eeprom(i, f) = endswith(f, ".elf") ? eeprom_elf(i, f) : eeprom_hex(i, f)

read(isp, type, file) = avrdude(isp, `-U $type:r:$file:i`)

fuse(isp, fuse) = avrdude(isp, `-U $fuse`)

end # module


text:
using LinearAlgebra
using BenchmarkTools
using Statistics
using PyPlot
# pygui(true);

function mathieu(z, zτ, t, κ, δ, ϵ, b, τ, ω)
        dzdt=zeros(ComplexF64,size(z))
                dzdt[1]=z[2]
                dzdt[2]=-(δ+(ϵ * cos(ω * t)))*z[1] - κ * z[2] + b * zτ[1]
                return dzdt
end

#integration

function Sinterp(Sloc::Array{ComplexF64,2},indf::Float64)
        i_sub=Int64(floor(indf))::Int64
        remi=(indf-i_sub)::Float64
        # if abs(remi)<(eps(typeof(remi))*5.0)
        #         return S[:,i_sub] .* (1-remi)
        # else
        return Sloc[:,i_sub] .* (1.0-remi) .+ Sloc[:,i_sub+1]*remi
        # end
end

function Sinterp(S::Array{ComplexF64,2},i::Integer)
        return S[:,i]
end

# S=rand(ComplexF64,1,n+1+m)
# @time Sinterp(S,4.3)
# @time Sinterp(S,4)
# @time S[:,4]
# figure(1)
# plot(1:n+1+m,S[1,:])
# v=1:0.1:n+1+m
# scatter(v,[Sinterp(S[:,:,1],i) for i in v])

function mathiestab( κ, δ, ϵ, b, τ, ω, n,m,mult,iter,RKtype::String)

dt=τ/n
S=rand(ComplexF64,2,n+1+m,mult)
λj=Array{ComplexF64}(undef,mult)
for kinter in 1:iter
        #Threads.@threads
        for mi in 1:mult
                for kt in n+1:n+m
                if RKtype=="EE"
                        k1=mathieu(S[:,kt,mi], Sinterp(S[:,:,mi],kt-n),(kt-(n+1))*dt, κ, δ, ϵ, b, τ, ω)
                        k=k1 ./ 1.0
                elseif RKtype=="RK2-fast"
                         k1=mathieu(S[:,kt,mi]      ,S[:,kt-n  ,mi],(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                         k2=mathieu(S[:,kt,mi]+k1*dt,S[:,kt-n+1,mi],(kt-(n+1)+1)*dt, κ, δ, ϵ, b, τ, ω)
                         k=(k1+k2) ./ 2.0
                elseif RKtype=="RK2"
                        # RK2
                        k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        k2=mathieu(S[:,kt,mi]+k1*dt,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+k2) ./ 2.0
                elseif RKtype=="RK3"
                        # RK3
                        k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        k2=mathieu(S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                        k3=mathieu(S[:,kt,mi]-k1*dt*1.0+k2*dt*2.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+4*k2+k3) ./ 6.0
                elseif RKtype=="RK4"
                         # RK4
                         k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                         k2=mathieu(S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                         k3=mathieu(S[:,kt,mi]+k2*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                         k4=mathieu(S[:,kt,mi]+k3*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                         k=(k1+2*k2+2*k3+k4) ./ 6.0
                else
                        error("nincs ilyen integrator")
                end
                S[:,kt+1,mi]=S[:,kt,mi]+dt*k;
                end
        end

        Sj=reshape(S[:,1:n+1,:],(2*(n+1),mult))
        Vj=reshape(S[:,end-n:end,:],(2*(n+1),mult))
        Hj=Sj\Vj
        (λj,Gj)=eigen(Hj)

        S[:,1:n+1,:] = (Vj*Gj)[:]

        # figure(3)
        # clf()
        # plot((1:n+1+m).*dt,S[2,1:n+1+m,:])
        # plot((1:n+1).*dt,S[2,1:n+1,:])
        # pause(0.2)
        #
        # figure(4)
        # scatter(kinter .* ones(mult),abs.(λj))
end

return λj
end

println("----------=======<<<<<<<<<<<<<################>>>>>>>>>>>>>>========----------")
κ=0.05
δ=5
ϵ=1.0
b=0.1
τ=2pi
ω=1
#
mult=20 # number of computed eigenvalues
iter=30 # number of iteration in issi

# @debug
# @elapsed mu_ref=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω,10^6,m,mult,iter,"RK4")))
@time  mu_ref=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω, 2^4, Int64(2^4*(1 / ω)),mult,iter,"RK4")))
@time  mu_ref=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω, 2^4, Int64(2^4*(1 / ω)),mult,iter,"RK4")))

n=2^7#16
m=Int64(n*(1 / ω))
@time mu_ref=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω, n,m,mult,iter,"RK4")))

integrators=["EE","RK2","RK3","RK4"]
ni=2 .^ (3:6)
mumax=zeros(length(ni),length(integrators))
timelength=zeros(length(ni),length(integrators))
timestd=zeros(length(ni),length(integrators))

figure(1)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("n / error")
xlabel("n")
ylabel("error")
#
figure(2)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("n / time")
xlabel("n")
ylabel("time")
#
figure(3)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("time / error")
xlabel("time")
ylabel("error")

Nsample=10

for kn in 1:length(ni)

        Threads.@threads for ksample in 1:Nsample
                println("<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>")
                println([ksample,kn,ni[kn]])
                for kRKi in 1:length(integrators)
                        integratortype=integrators[kRKi]
                        print(integratortype)
                        print("  -  ")
        n=ni[kn]
        m=Int64(n*(1 / ω))
        # timeing = @benchmark (mumax[kn]=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω,n,m,mult,iter,"RK4"))))
        # timelength[kn]=mean(timeing.times)
        # timestd[kn]=std(timeing.times)
        timeing = @elapsed (mumax[kn,kRKi]=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω,n,m,mult,iter,integratortype))))
        timelength[kn,kRKi]+=timeing
        println(timeing)

        figure(1)
        scatter(n,abs(mumax[kn,kRKi] - mu_ref))
        figure(2)
        scatter(n,timeing)
        figure(3)
        scatter(timeing,abs(mumax[kn,kRKi] - mu_ref))
                end
        end

figure(1)
plot(ni[1:kn,:],abs.(mumax[1:kn,:] .- mu_ref))
savefig("1_nt")
# savefig("1_t")

figure(2)
plot(ni[1:kn,:],timelength[1:kn,:]./Nsample)
savefig("2_nt")
# savefig("2_t")

figure(3)
plot(timelength[1:kn,:]./Nsample,abs.(mumax[1:kn,:] .- mu_ref))
savefig("3_nt")
# savefig("3_t")
end

println("-------------------")
                # figure(1)
                # plot(log.(ni)./log(2),log.(abs.(mumax .- mu_ref))./log(2))
                # figure(2)
                # plot(log.(ni)./log(2),log.(timelength)./log(2))
                # # plot(log.(timelength .+ timestd))
                # # plot(log.(timelength .- timestd))

# #----------------------------
# using MDBM
#
# foo=(x,y)->(maximum(abs.(mathiestab( κ, x, ϵ, y, τ, ω,n,m,mult,iter,"RK4")))-1)
#
# δv=-0.1:0.2:10
# # ϵv=-3:0.1:3
# bv=-1:0.2:1
#
#
# mymdbm=MDBM_Problem(foo,Axis.([δv,bv]))
# MDBM.solve!(mymdbm,3)
#
# x_eval,y_eval=getevaluatedpoints(mymdbm)
# x_sol,y_sol=getinterpolatedsolution(mymdbm)
#
# fig = figure(11);clf()
# PyPlot.scatter(x_eval,y_eval,s=5)
# PyPlot.scatter(x_sol,y_sol)
#
# println("cuccma")


text:
# This file defines the CompoundState type. This type defines a state which
#   is a collection of several base states.

# The CompoundState type requires no other types.
requiredTypes = []

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export CompoundState
"""
This type defines a compound state which groups a number of states in a manpower
simulation. This is used to generate a hierarchy of levels that personnel
members can occupy. For example: assume that there are states for each
combination of service branch (Land, Air, Marine, Medical) and grade category
(Volunteer, Non-com, Officer). The compound states can then group all personnel
members by grade category, so each compound state corresponds to a service
branch.

The type contains the following fields:
* `name::String`: the name of the compound state.
* `stateTarget::Int`: the target number of personnel members in this state.
* `stateList::Vector{String}`: the names of the base states the compound state
  consists of. Remark that the simulation must ensure that all the component
  base states actually exist.
"""
type CompoundState

    name::String
    stateTarget::Int
    stateList::Vector{String}

    function CompoundState( name::String )

        newCompState = new()
        newCompState.name = name
        newCompState.stateTarget = -1
        newCompState.stateList = Vector{String}()
        return newCompState

    end  # CompoundState( name )

end  # CompoundState


text:
@testset "events" begin

let
    start = CuEvent()
    stop = CuEvent()
    @test start != stop

    record(start)
    record(stop)
    synchronize(stop)

    @test elapsed(start, stop) > 0
end

@test (CUDAdrv.@elapsed begin end) > 0

@test (CUDAdrv.@elapsed CuDefaultStream() begin end) > 0

CuEvent(CUDAdrv.EVENT_BLOCKING_SYNC)
CuEvent(CUDAdrv.EVENT_BLOCKING_SYNC | CUDAdrv.EVENT_DISABLE_TIMING)

@testset "stream wait" begin
    event  = CuEvent()
    stream = CuStream()

    CUDAdrv.record(event, stream)

    CUDAdrv.wait(event)
    synchronize()
end

@testset "elapsed stream" begin
    stream = CuStream()
    @test (CUDAdrv.@elapsed stream begin end) > 0
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    NearestMapping

A mapping method in which data points are assigned to their nearest
point in the domain.
"""
struct NearestMapping <: MappingMethod end

function map(sdata, sdomain, targetvars, ::NearestMapping)
  N = ncoords(sdata)
  T = coordtype(sdata)

  @assert targetvars ⊆ name.(variables(sdata)) "target variables must be present in spatial data"

  # dictionary with mappings
  mappings = Dict(var => Dict{Int,Int}() for var in targetvars)

  # pre-allocate memory for coordinates
  coords = MVector{N,T}(undef)

  # nearest neighbor search method
  neighbor = Vector{Int}(undef, 1)
  searcher = KNearestSearch(sdomain, 1)

  for ind in 1:nelms(sdata)
    # update datum coordinates
    coordinates!(coords, sdata, ind)

    # find nearest location in the domain
    search!(neighbor, coords, searcher)

    # save pair if there is data for variable
    for var in targetvars
      if !ismissing(sdata[var][ind])
        push!(mappings[var], neighbor[1] => ind)
      end
    end
  end

  mappings
end


text:
t(t0, tf, dt) = t0:dt:tf
t(dom::Domain) = t(dom.initial_year,dom.final_year,dom.dt)
t(params::ClimateModelParameters) = t(params.domain)
t(m::ClimateModel) = t(m.domain)

future_mask(t, present_year) = t .<= present_year
future_mask(model) = future_mask(t(model), model.domain.present_year)


text:
module CachedArrays

#below deps are just for utility functions (may move those?)
using Images, Unitful
using AxisArrays
const axes = Base.axes

import Base: size, getindex, setindex!, show, eltype, parent, axes

export CachedArray,
        cached_axes,
        noncached_axes,
        axisspacing,      #utility
        match_axisspacing #utility

include("util.jl")
include("abstract_cached_array.jl")
include("cached_array.jl")

end # module


text:
using Documenter, DocumenterTools, DocumenterLaTeX

using TSML
using TSML.DecisionTreeLearners

using Test

const ROOT = joinpath(@__DIR__, "..")

# Documenter package docs
doc = makedocs(
    debug = true,
    root = ROOT,
    build = "pdf/build",
    modules = [TSML, DecisionTreeLearners],
    clean = false,
	format = LaTeX(platform = "docker"),
	sitename = "TSML Documentation",
	authors = "Paulito P. Palmes",
	pages = Any[
				"HOME" => "index.md",
				"Tutorial" => Any[
								  "tutorial/aggregators.md",
								  "tutorial/pipeline.md",
								  "tutorial/statistics.md",
								  "tutorial/tsdetectors.md"
								 ],
				"Manual" => Any[
								"Date Processing" => "man/dateproc.md",
								"Value Processing" => "man/valueproc.md",
								"Aggregation" => "man/aggregation.md",
								"Imputation" => "man/imputation.md",
								"Monotonic Detection" => "man/monotonic.md",
								"TS Classification" => "man/tsclassification.md",
								"CLI Wrappers" => "man/cli.md"
							   ],
				"Library" => Any[
								 "Decision Tree" => "lib/decisiontree.md"
								 #"Scikit Learners" => "lib/sklearn.md",
								 #"Caret Learners" => "lib/caretlearn.md"
								]
	]
);

# hack to only deploy the actual pdf-file
mkpath(joinpath(ROOT, "pdf", "build", "pdfdir"))
let files = readdir(joinpath(ROOT, "pdf", "build"))
    for f in files
        if startswith(f, "Documenter.jl") && endswith(f, ".pdf")
            mv(joinpath(ROOT, "pdf", "build", f),
               joinpath(ROOT, "pdf", "build", "pdfdir", f))
        end
    end
end

deploydocs(
    repo = "github.com/IBM/TSML.jl.git",
    root = ROOT,
    target = "pdf/build/pdfdir",
    branch = "master",
    forcepush = true,
)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_tab

"""
    dcc_tab(;kwargs...)
    dcc_tab(children::Any;kwargs...)
    dcc_tab(children_maker::Function;kwargs...)

A Tab component.
Part of dcc.Tabs - this is the child Tab component used to render a tabbed page.
Its children will be set as the content of that tab, which if clicked will become visible.
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `label` (String; optional): The tab's label
- `value` (String; optional): Value for determining which Tab is currently selected
- `disabled` (Bool; optional): Determines if tab is disabled or not - defaults to false
- `disabled_style` (Dict; optional): Overrides the default (inline) styles when disabled
- `disabled_className` (String; optional): Appends a class to the Tab component when it is disabled.
- `className` (String; optional): Appends a class to the Tab component.
- `selected_className` (String; optional): Appends a class to the Tab component when it is selected.
- `style` (Dict; optional): Overrides the default (inline) styles for the Tab component.
- `selected_style` (Dict; optional): Overrides the default (inline) styles for the Tab component when it is selected.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function dcc_tab(; kwargs...)
        available_props = Symbol[:children, :id, :label, :value, :disabled, :disabled_style, :disabled_className, :className, :selected_className, :style, :selected_style, :loading_state]
        wild_props = Symbol[]
        return Component("dcc_tab", "Tab", "dash_core_components", available_props, wild_props; kwargs...)
end

dcc_tab(children::Any; kwargs...) = dcc_tab(;kwargs..., children = children)
dcc_tab(children_maker::Function; kwargs...) = dcc_tab(children_maker(); kwargs...)

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
path = dirname(@__FILE__)

# Set up arguments
global m = AnSchorfheide(testing = true)
m <= Setting(:impulse_response_horizons, 10)
m[:σ_z].value = sqrt(.0201332) # need to change parameters or you get a PositiveDefiniteException
m[:σ_g].value = sqrt(1.09725)  # in the long run identification
m[:σ_R].value = sqrt(.0590387)
m[:ρ_z].value = 0.99
m[:ρ_g].value = 0.99
m[:ρ_R].value = 0.99
θ = map(x -> x.value, m.parameters)
DSGE.update!(m, θ)
system = compute_system(m)
θmat = convert(Matrix{Float64}, hcat(θ, θ)')
obs_shock = zeros(n_observables(m))
obs_shock[1] = 1.

@testset "Wrapper for Cholesky-identified DSGE impulse responses" begin
    for do_parallel in [false, true]
        for do_flip in [false, true]
            for method in [:cholesky, :choleskyLR, :cholesky_long_run]
                states_chol, obs_chol, pseudo_chol =
                    impulse_responses(system, impulse_response_horizons(m),
                                      Matrix{Float64}(I, n_observables(m), n_observables(m)),
                                      obs_shock, flip_shocks = do_flip, restriction = method)

                states1, obs1, pseudo1 = impulse_responses(m, θ, :mode, method, 1;
                                                           flip_shocks = do_flip,
                                                           parallel = do_parallel)
                states2, obs2, pseudo2 = impulse_responses(m, θmat, :full, method, 1,
                                                           flip_shocks = do_flip,
                                                           parallel = do_parallel)
                @test @test_matrix_approx_eq states_chol states1[:, :, 1]
                @test @test_matrix_approx_eq obs_chol obs1[:, :, 1]
                @test @test_matrix_approx_eq pseudo_chol pseudo1[:, :, 1]
                @test @test_matrix_approx_eq states_chol states2[:, :, 1]
                @test @test_matrix_approx_eq obs_chol obs2[:, :, 1]
                @test @test_matrix_approx_eq pseudo_chol pseudo2[:, :, 1]
                @test @test_matrix_approx_eq states_chol states2[:, :, 2]
                @test @test_matrix_approx_eq obs_chol obs2[:, :, 2]
                @test @test_matrix_approx_eq pseudo_chol pseudo2[:, :, 2]
            end
        end
    end

    # No tests but check that there is no error if we create a MeansBands
    mb = impulse_responses(m, θmat, :full, :cholesky, 1,
                           create_meansbands = true,
                           test_meansbands = true)
    mb = impulse_responses(m, θmat, :full, :choleskyLR, 1,
                           create_meansbands = true,
                           test_meansbands = true)
    mb = impulse_responses(m, θmat, :full, :maxBC, 1,
                           create_meansbands = true,
                           test_meansbands = true)
end

@testset "Wrapper for `maxBC`-identified DSGE impulse responses" begin
    for do_parallel in [false, true]
        for do_flip in [false, true]
            for method in [:maxBC, :maximum_business_cycle_variance]
                states_maxBC, obs_maxBC, pseudo_maxBC =
                    impulse_responses(system, impulse_response_horizons(m),
                                      (2 * π / 32, 2 * π / 6), 1;
                                      flip_shocks = do_flip)

                states1, obs1, pseudo1 = impulse_responses(m, θ, :mode, method, 1,
                                                           flip_shocks = do_flip,
                                                           parallel = do_parallel)
                states2, obs2, pseudo2 = impulse_responses(m, θmat, :full, method, 1,
                                                           flip_shocks = do_flip,
                                                           parallel = do_parallel)
                @test @test_matrix_approx_eq states_maxBC states1[:, :, 1]
                @test @test_matrix_approx_eq obs_maxBC obs1[:, :, 1]
                @test @test_matrix_approx_eq pseudo_maxBC pseudo1[:, :, 1]
                @test @test_matrix_approx_eq states_maxBC states2[:, :, 1]
                @test @test_matrix_approx_eq obs_maxBC obs2[:, :, 1]
                @test @test_matrix_approx_eq pseudo_maxBC pseudo2[:, :, 1]
                @test @test_matrix_approx_eq states_maxBC states2[:, :, 2]
                @test @test_matrix_approx_eq obs_maxBC obs2[:, :, 2]
                @test @test_matrix_approx_eq pseudo_maxBC pseudo2[:, :, 2]
            end
        end
    end

    # No tests but check that there is no error if we create a MeansBands
    mb = impulse_responses(m, θmat, :full, :cholesky, 1,
                           create_meansbands = true,
                           test_meansbands = true)
    mb = impulse_responses(m, θmat, :full, :choleskyLR, 1,
                           create_meansbands = true,
                           test_meansbands = true)
    mb = impulse_responses(m, θmat, :full, :maxBC, 1,
                           create_meansbands = true,
                           test_meansbands = true)
end

nothing


text:
module Quantum
	#####################################
	#constants###########################
	#####################################
	const lang = "\u27E8"
	const rang = "\u27E9"
	const otimes = "\u2297"
	const vdots ="\u205E"

	#####################################
	#imports#############################
	#####################################
	import Base:
			length,
			size,
			endof,
			isequal,
			==,
			copy,
			getindex,
			ctranspose,
			show,
			zero,
			one,
			conj,
			eltype,
			hash,
		    showcompact,
		    eltype,
		    get,
		    in,
		    setdiff,
		    +,
		    -,
		    *,
		    /,
		    ^,
		    (.*),
		    (.+),
		    (.-),
		    (./),
		    (.^),
		    abs,
		    exp,
		    map,
		    filter,
		    find,
		    ndims,
		    findn,
		    findnz,
		    nnz,
		    countnz,
		    ndims,
		    norm,
		    summary,
		    kron,
		    trace,
		    promote_rule,
		    convert

	#####################################
	#abstract types######################
	#####################################
	abstract Dirac
	abstract AbstractScalar <: Dirac
	abstract AbstractOperator <: Dirac
	abstract State{S} <: Dirac

	typealias DiracCoeff Union(Number, AbstractScalar)

	#####################################
	#includes############################
	#####################################
	include("crossjoin.jl")
	include("state.jl")
	include("products.jl")
	include("basis.jl")
	include("scalar.jl")
	include("diracvector.jl")
	include("diracmatrix.jl")
	include("misc.jl")
	include("fock.jl")

	#####################################
	#exports#############################
	#####################################
	export Dirac,
		   State,
		   Bra,
		   Ket,
		   Tensor,
		   AbstractBasis,
		   AbstractScalar,
		   InnerProduct,
		   ScalarExpr,
		   Basis,
		   DiracVector,
		   DiracMatrix,
		   tensor,
		   basis,
		   dual,
		   label,
		   bsym,
		   separate,
		   isdual,
		   labeldelta,
		   inner,
		   samebasis,
		   kind,
		   bcat,
		   findstates,
		   filtercoeffs,
		   filterstates,
		   getpos,
		   normalize,
		   qeval,
		   dvec,
		   svec,
		   dmat,
		   bra,
		   ket,
		   fbasis,
		   fvec,
		   fcreate,
		   fdestroy,
		   fnum,
		   feye,
		   commutator,
		   ptrace,
		   scalar,
		   actop,
		   filternz,
		   showsp
end


text:
#
# Calculations for the Circular Restricted Three Body Problem.
#

"""
Normalizes a CR3BP orbit in the rotating reference frame.
"""
function LinearAlgebra.normalize(r::AbstractVector{<:Number}, v::AbstractVector{<:Number}, t::Number, a::Number, μs::Tuple{<:Number, <:Number}; lengthunit=u"km", timeunit=u"s")

    rₙ = r ./ a
    Tₛ = period(a, sum(μs))
    vₙ = v ./ (a / Tₛ)
    tₙ = t / Tₛ
    μₙ = min(μs[1], μs[2]) / (μs[1]+μs[2])

    DU = a  isa Unitful.Length ? a : a * lengthunit
    DT = Tₛ isa Unitful.Time ? Tₛ : Tₛ * timeunit

    return rₙ, vₙ, tₙ, μₙ, DU, DT

end

"""
Redimensionalizes a CR3BP orbit in the rotating reference frame.
"""
function redimension(rₙ::AbstractVector{<:Real}, vₙ::AbstractVector{<:Real}, tₙ::Real, μₙ::Real, DU::Number, TU::Number)

    r = rₙ .* DU
    v = vₙ .* DU ./ TU
    t = tₙ * TU

    sum_μs = DU^3 / ((TU / 2π)^2)
    μ₂ = μₙ * sum_μs
    μ₁ = sum_μ - μ₂

    return r, v, t, DU, (μ₁,μ₂)

end

"""
Returns the spacecraft's nondimensional position w.r.t. body 1 (or 2).
"""
nondimensional_radius(r, xᵢ) = √( (r[1]-xᵢ)^2 + r[2]^2 + r[3]^2 )

"""
Returns synodic distance to primary body.
"""
distance_to_primary(orbit::CR3BPOrbit) = norm(OrbitalStates.get_r(state(orbit)) .- SVector(-massparameter(orbit), 0, 0))

"""
Returns synodic distance to secondary body.
"""
distance_to_secondary(orbit::CR3BPOrbit) = norm(OrbitalStates.get_r(state(orbit)) .- SVector(1-massparameter(orbit), 0, 0))

"""
Returns the potential energy `U` in the Synodic frame with Normalized units.
"""
potential_energy(r, μ) = (r[1]^2 + r[2]^2)/2 + ((1-μ)/nondimensional_radius(r,-μ)) + (μ/nondimensional_radius(r,1-μ))

"""
Returns the potential energy `U` in the Synodic frame with Normalized units.
"""
function potential_energy(orbit::CR3BPOrbit)
    return potential_energy(OrbitalStates.get_r(state(orbit)), massparameter(system(orbit)))
end

"""
Returns the Jacobi Constant, `C` in the Synodic frame with Normalized units.
"""
jacobi_constant(r, v, μ) = 2*potential_energy(r, μ) - (v⋅v)

"""
Returns the Jacobi Constant, `C` in the Synodic frame with Normalized units.
"""
function jacobi_constant(orbit::CR3BPOrbit)
    return jacobi_constant(OrbitalStates.get_r(state(orbit)), OrbitalStates.get_v(state(orbit)) massparameter(system(orbit)))
end

"""
Given the Synodic frame vector, returns the vector in the barycenteric-inertial reference frame.
"""
function inertial(vecₛ::AbstractVector, t, ω::Unitful.AbstractQuantity=1.0u"rad"/unit(t))

    θ = ω*t
    ᴵTₛ = @SMatrix [
        cos(θ) sin(θ) 0
       -sin(θ) cos(θ) 0
        0      0      1
    ]

    return  ᴵTₛ * vecₛ

end

"""
Given an `InertialCartesianState`, returns the state in the synodic (rotating) reference frame.
"""
function synodic(state::AbstractVector, t)
    θ = ω*t
    ˢTᵢ = inv(SMatrix{3,3}([
        cos(θ) sin(θ) 0
       -sin(θ) cos(θ) 0
        0      0      1
    ]))

    return ˢTᵢ * state
end

"""
Position of primary body.
"""
function primary_position(orbit::CR3BPOrbit)
    μ = massparameter(orbit)
    SVector{3}(
        -μ,
        zero(μ),
        zero(μ)
    )
end

"""
Position of primary body.
"""
function primary_position(orbit::CR3BPOrbit)
    μ = massparameter(orbit)
    SVector{3}(
        one(μ)-μ,
        zero(μ),
        zero(μ)
    )
end

"""
Returns the lagrange points for a CR3BP system.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `L`: Langrange points requested, must be in range [1,5]

__Outputs:__
- Tuple of Lagrange points
- Throws `ArgumentError` if L is out of range [1,5]

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/)
"""
function lagrange(μ::Real, L=1:5)

    if !all(L[i] ∈ (1,2,3,4,5) for i ∈ 1:length(L))
        throw(ArgumentError("Requested lagrange points must be in range [1,5]"))
    end

	expressions = @SVector [
		x -> x - (1-μ)/(x+μ)^2 + μ/(x+μ-1)^2,
		x -> x - (1-μ)/(x+μ)^2 - μ/(x+μ-1)^2,
		x -> x + (1-μ)/(x+μ)^2 + μ/(x+μ+1)^2
	]

	return  (map(f->[find_zero(f, (-3,3)), 0, 0], expressions)...,
			[(1/2) - μ, √(3)/2, 0], [(1/2) - μ, -√(3)/2, 0])[L]

end

"""
Returns the lagrange points for a CR3BP system.

__Arguments:__
- `sys`: Non-dimensional mass parameter for the CR3BP system.
- `L`: Langrange points requested, must be in range [1,5]

__Outputs:__
- Tuple of Lagrange points
- Throws `ArgumentError` if L is out of range [1,5]

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/)
"""
function lagrange(sys::CR3BPParameters, L=1:5)

    return lagrange(massparameter(sys), L)
end

"""
Returns an analytical solution for a Halo orbit about `L`.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `Az`: Desired non-dimensional Z-amplitude for Halo orbit.
- `ϕ`: Desired Halo orbit phase.
- `steps`: Number of non-dimensional timepoints in returned state.
- `L`: Lagrange point to orbit (L1 or L2).
- `hemisphere`: Specifies northern or southern Halo orbit.

__Outputs:__
- Synodic position vector `r::Array{<:AbstractFloat}`
- Synodic velocity vector `v::Array{<:Abstractfloat}`.
- Halo orbit period `Τ`.
- Throws `ArgumentError` if L is not `1` or `2`.

__References:__
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/).
"""
function analyticalhalo(μ; Az=0.00, ϕ=0.0, steps=1,
                       L=1, hemisphere=:northern)

    if L == 1
        point = first(lagrange(μ, 1))
        γ = abs(one(μ) - μ - point)
        n = collect(1:4) .* one(μ)
        c = @. (μ + (-one(1))^n * (one(μ)-μ)γ^(n+1)) / (γ^3 * (one(μ) - γ^(n+1)))
    elseif L == 2
        point = first(lagrange(μ, 2))
        γ = abs(point - one(μ) + μ)
        n = collect(1:4) .* one(μ)
        c = @. ((-one(μ))^n * μ + (-one(μ))^n * (one(μ)-μ)γ^(n+1)) / (γ^3 * (one(μ) + γ^(n+1)))
    else
        throw(ArgumentError("Only Halo orbits about L1 or L2 are supported."))
    end

    ωₚ  = √((2 - c[2] + √((9c[2]^2 - 8c[2])))/2)
    k   = (ωₚ^2 + 1 + 2c[2]) / (2ωₚ)

    d₁  = (3ωₚ^2 / k) * (k*(6ωₚ^2 - 1) - 2ωₚ)
    d₂  = (8ωₚ^2 / k) * (k*(11ωₚ^2 - 1) - 2ωₚ)
    a₂₁ = (3c[3] * (k^2 - 2)) / (4(1 + 2c[2]))
    a₂₂ = (3c[3]) / (4(1 + 2c[2]))
    a₂₃ = (-3c[3]ωₚ / (4k*d₁)) * (3k^3 * ωₚ - 6k*(k-ωₚ) + 4)
    a₂₄ = (-3c[3]ωₚ / (4k*d₁)) * (2 + 3k*ωₚ)
    b₂₁ = (-3c[3]ωₚ / (2d₁)) * (3k*ωₚ - 4)
    b₂₂ = -3c[3]*ωₚ / d₁
    d₂₁ = -c[3] / (2ωₚ^2)
    a₃₁ = (-9ωₚ / (4d₂)) * (4c[3] * (k*a₂₃-b₂₁) + k*c[4]*(4+k^2)) +
          ((9ωₚ^2 + 1 - c[2]) / (2d₂)) * (3c[3]*(2a₂₃-k*b₂₁) + c[4]*(2+3k^2))
    a₃₂ = (-9ωₚ / (4d₂)) * (4c[3] * (3k*a₂₄-b₂₂) + k*c[4]) -
          (3 / (2d₂)) * (9ωₚ^2 + 1 - c[2]) * (c[3]*(k*b₂₂+d₂₁-2a₂₄) - c[4])
    b₃₁ = (3 / (8d₂)) * 8ωₚ * (3c[3] * (k*b₂₁ - 2a₂₃) - c[4]*(2+3k^2)) +
          (3/(8d₂)) * (9ωₚ^2 + 1 + 2c[2]) * (4c[3]*(k*a₂₃-b₂₁) + k*c[4]*(4+k^2))
    b₃₂ = (9ωₚ/d₂)*(c[3]*(k*b₂₂+d₂₁-2a₂₄)-c[4]) +
          (3(9ωₚ^2 + 1 + 2c[2]) / (8d₂) * (4c[3]*(k*a₂₄-b₂₂)+k*c[4]))
    d₃₁ = (3 / (64ωₚ^2)) * (4c[3]*a₂₄ + c[4])
    d₃₂ = (3 / (64 + ωₚ^2)) * (4c[3]*(a₂₃ - d₂₁) + c[4]*(4+k^2))

    s₁  = (1 / (2ωₚ*(ωₚ*(1+k^2) - 2k))) *
          (3c[3]/2 * (2a₂₁*(k^2 - 2) - a₂₃*(k^2 + 2) - 2k*b₂₁) -
            (3c[4]/8) * (3k^4 - 8k^2 + 8))
    s₂  = (1 / (2ωₚ*(ωₚ*(1+k^2) - 2k))) *
          (3c[3]/2 * (2a₂₂*(k^2-2) + a₂₄*(k^2 + 2) + 2k*b₂₂ + 5d₂₁) +
            (3c[4]/8) * (12 - k^2))
    l₁  = (-3c[3] / 2) * (2a₂₁ + a₂₃ + 5d₂₁) - (3c[4]/8)*(12 - k^2) + 2ωₚ^2 * s₁
    l₂  = (3c[3]/2) * (a₂₄ - 2a₂₂) + (9c[4]/8) + 2ωₚ^2 * s₂
    Δ   = ωₚ^2 - c[2]

    Aᵧ  = Az / γ
    Aₓ  = √((-l₂*Aᵧ^2 - Δ) / l₁)

    ν   = 1 + s₁*Aₓ^2 + s₂*Aᵧ^2
    Τ   = 2π / (ωₚ*ν)
    τ   = ν .* (steps > 1 ? range(0, stop=Τ, length=steps) : range(0, stop=Τ, length=1000))

    if hemisphere == :northern
        m = 1.0
    elseif hemisphere == :southern
        m = 3.0
    else
        throw(ArgumentError("`hemisphere` must be `:northern` or `:southern`."))
    end

    δₘ = 2 - m
    τ₁ = @. ωₚ*τ + ϕ

    x = @. γ * (a₂₁*Aₓ^2 + a₂₂*Aᵧ^2 - Aₓ*cos(τ₁) + (a₂₃*Aₓ^2 -
                    a₂₄*Aᵧ^2)*cos(2τ₁) + (a₃₁*Aₓ^3 - a₃₂*Aₓ*Aᵧ^2)*cos(3τ₁)) + 1 - μ - (L == 1 ? γ : -γ)
    y = @. γ * (k*Aₓ*sin(τ₁) + (b₂₁*Aₓ^2 - b₂₂*Aᵧ^2)*sin(2τ₁) +
                    (b₃₁*Aₓ^3 - b₃₂*Aₓ*Aᵧ^2)*sin(3τ₁))
    z = @. γ * (δₘ*Aᵧ*cos(τ₁) + δₘ*d₂₁*Aₓ*Aᵧ*(cos(2τ₁)-3) +
                    δₘ*(d₃₂*Aᵧ*Aₓ^2 - d₃₁*Aᵧ^3)*cos(3τ₁))

    ẋ = @. γ * (ωₚ*ν*Aₓ*sin(τ₁) - 2ωₚ*ν*(a₂₃*Aₓ^2 - a₂₄*Aᵧ^2)*sin(2τ₁) -
                    3ωₚ*ν*(a₃₁*Aₓ^3 - a₃₂*Aₓ*Aᵧ^2)*sin(3τ₁))
    ẏ = @. γ * (ωₚ*ν*k*Aₓ*cos(τ₁) + 2ωₚ*ν*(b₂₁*Aₓ^2 - b₂₂*Aᵧ^2)*cos(2τ₁) +
                    3ωₚ*ν*(b₃₁*Aₓ^3 - b₃₂*Aₓ*Aᵧ^2)*cos(3τ₁))
    ż = @. γ * (-ωₚ*ν*δₘ*Aᵧ*sin(τ₁) - 2ωₚ*ν*δₘ*d₂₁*Aₓ*Aᵧ*sin(2τ₁) -
                    3ωₚ*ν*δₘ*(d₃₂*Aᵧ*Aₓ^2 - d₃₁*Aᵧ^2)*sin(3τ₁))

    return hcat(x, y, z)[1:steps, :], hcat(ẋ, ẏ, ż)[1:steps, :], Τ

end

"""
Returns a `Vector` of `Matrix` values.
Each `Matrix` contains a 3-column nondimensional
position trajectory in the Synodic frame which
represents a Zero Velocity Curve.
"""
function zerovelocity_curves(r::AbstractVector, v::AbstractVector, μ::Real;
                             nondimensional_range = range(-2; stop=2, length=1000))

    x = nondimensional_range
    y = nondimensional_range

    r = orbit |> synodic |> normalize |> position_vector
    v = orbit |> synodic |> normalize |> velocity_vector
    μ  = normalized_mass_parameter(orbit.system)

    z = [
        2 * potential_energy([xs, ys, 0.0], μ) - jacobi_constant(r, v, μ)
        for xs ∈ x, ys ∈ y
    ]

    zvc_contour = contours(x,y,z,[0.0])

    curves = Vector{Matrix{Float64}}()

    for zvc_level ∈ Contour.levels(zvc_contour)

        x = Float64[]
        y = Float64[]
        for zvc_line ∈ Contour.lines(zvc_level)
            xs, ys = coordinates(zvc_line)
            if length(x) == length(y) == 0
                x = xs
                y = ys
            else
                x = vcat(x, xs)
                y = vcat(y, ys)
            end
        end

        if length(curves) == 0
            curves = [hcat(x,y)]
        else
            curves = vcat(curves, [hcat(x, y)])
        end
    end

    return curves

end


text:
#This file has the definitions for the gray Flux-Limited Diffusion Solver
#The code uses units of
# cm for length
# sh for time (1 sh = 10^-8 s)
# GJ (gigajoules) for energy; 1 GJ is also known as a jerk
# keV for temperature (1 keV is 1.1605 * 10^7 K)

using PyCall
using JLD, HDF5

#radiation constant in GJ/keV^4/cm^3
const a = 0.01372
#speed of light in cm/sh; 1 sh is 10^-8 s
const c = 299.98
const Planck_int_const = 15*a*c/(4*pi^5)

"This function maps an Nx*Ny 1-D array to a i,j on an Nx by Ny grid"
function vect_to_xy(l,Nx,Ny)
    j = div(l,Nx) + 1
    i = mod(l-j,Nx) + 1
    (i,j)
end

"This function maps i,j on an Nx by Ny grid to an Nx*Ny 1-D array"
function xy_to_vect(i,j,Nx,Ny)
l = Nx*(j-1) + i
l
end
function nm(f, fp, x; tol=sqrt(eps()))

    ctr, max_steps = 0, 100

    while (abs(f(x)) > tol) && ctr < max_steps
        x = x - f(x) / fp(x)
        ctr = ctr + 1
    end

    ctr >= max_steps ? error("Method did not converge.") : return (x, ctr)

end

"This function computes the integral of the planck function from E1 = h*nu [keV] to infinity at T[keV]"
function B1(E1,T)
    r = (E1./T)

    num_terms = convert(Int32,round(minimum([2+20/minimum(r),1024])))
    output = 0*T
    for n in 1:num_terms
        invn = 1.0/n
        output += invn*exp.(-n*r).*(r.^3 + 3*r.^2*invn .+ 6*r * invn^2 .+ 6*invn^3)
    end
    return output*Planck_int_const
end

"This function gives Bg(T)"
function Bg(E1,E2,T)
    return B1(E1,T) - B2(E2,T)
end

"This function computes the derivative of B1 w.r.t. T"
function DB1DT(y,T)
    r = y./T
    num_terms = convert(Int32,round(minimum([2+20/minimum(r),1024])))
    output = 0*T
    for n in 1:num_terms
        output += exp.(-n*r).*(24*T.^4 .+ 24*n*T.^3 .*y + 12*n^2*T.^2 .*y.^2 + 4*n^3*T.*y.^3 .+ n^4*y.^4)./(n^4*T)
    end
    return output*Planck_int_const
end
"This function computes the derivative of Bg w.r.t. T"
function DBgDT(E1,E2,T)
    return DB1DT(E1,T) .- DB1DT(E2,T)
end
"""
    create_A(D,sigma,Nr,Nz,Lr,Lz; lower_z = "refl",upper_z="refl",upper_r="refl")

Compute the matrix and RHS for an RZ diffusion problem
# Arguments
- `D::Float64(Nr,Nz)`: matrix of diffusion coefficients
- `sigma::Float64(Nr,Nz)`: matrix of absorption coefficients
- `Nr::Integer`: number of cells in r
- `Nz::Integer`: number of cells in z
- `Lr::Float64`: size of domain in r
- `Lz::Float64`: size of domain in z
...
"""
function create_A(Dleft,Dright, Dtop,Dbottom, sigma,Nr,Nz,Lr,Lz; lower_z = "refl",upper_z="refl",upper_r="refl")

    redges = linspace(0,Lr,Nr+1)
    zedges = linspace(0,Lz,Nz+1)

    A = spzeros(Nr*Nz,Nr*Nz)
    b = zeros(Nr*Nz)
    dz = Lz/Nz
    dr = Lr/Nr
    idz = 1/dz
    idz2 = idz*idz
    idr = 1/dr

    #loop over cells
    for i in 1:Nr
        for j in 1:Nz
            here = xy_to_vect(i,j,Nr,Nz)
            Vij = pi*(redges[i+1]^2 - redges[i]^2)
            iVij = 1/Vij
            Siplus = 2*pi*redges[i+1]
            Siminus = 2*pi*redges[i]
            #add in diag
            A[here,here] = sigma[i,j]*c
            #add in i+1
            if (i< Nr)
                Diplus = 2*Dright[i,j]*Dleft[i+1,j]/(Dright[i,j] + Dleft[i+1,j])
                tmp_zone = xy_to_vect(i+1,j,Nr,Nz)
                A[here,tmp_zone] = -Siplus*iVij*idr*Diplus
                A[here,here] += Siplus*iVij*idr*Diplus
            elseif (upper_r == "vacuum")
                A[here,here] += Siplus*iVij*0.5*c #*idr*D[i,j]
            elseif (isa(upper_r,Float64))
                A[here,here] += Siplus*iVij*0.5*c #*idr*D[i,j]
                b[here] += Siplus*iVij*0.5*upper_r*c
            end
            #add in i-1
            if (i>1)
                Diminus = 2*Dleft[i,j]*Dright[i-1,j]/(Dleft[i,j] + Dright[i-1,j])
                tmp_zone = xy_to_vect(i-1,j,Nr,Nz)
                A[here,tmp_zone] = -Siminus*iVij*idr*Diminus
                A[here,here] += Siminus*iVij*idr*Diminus
            end
            #add in j+1
            if (j < Nz)
                Djplus = 2*Dtop[i,j]*Dbottom[i,j+1]/(Dtop[i,j] + Dbottom[i,j+1])
                tmp_zone = xy_to_vect(i,j+1,Nr,Nz)
                #println("l = $(tmp_zone) i = $(i), j = $(j+1)")
                A[here,tmp_zone] = -idz2*Djplus
                A[here,here] += idz2*Djplus
            elseif (upper_z == "vacuum")
                A[here,here] += idz*0.5*c   #idz2*D[i,j]
            elseif (isa(upper_z,Float64))
                A[here,here] += idz*0.5*c #idz2*D[i,j]
                b[here] += idz*upper_z*0.5*c #idz2*D[i,j]*upper_z
            end
            #add in j-1
            if (j>1)
                Djminus = 2*Dbottom[i,j]*Dtop[i,j-1]/(Dbottom[i,j] + Dtop[i,j-1])
                tmp_zone = xy_to_vect(i,j-1,Nr,Nz)
                A[here,tmp_zone] = -idz2*Djminus
                A[here,here] += idz2*Djminus
            elseif (lower_z == "vacuum")
                A[here,here] += idz*0.5*c #idz2*D[i,j]
            elseif (isa(lower_z,Float64))
                A[here,here] += idz*c*0.5 #idz2*D[i,j]
                b[here] += idz*lower_z*c #idz2*D[i,j]*lower_z
            end
        end #for j
    end #for i
    A,b
end #function

"This function solves equation 19 from Brunner's paper. It gives E_star_g"
function Estar_g(delta_t,T,Eg,Eg_l,D_func,sigma_func,Q_func,RHS,Nr,Nz,Lr,Lz;
    lower_z = "refl",upper_z="refl",upper_r="refl", LOUD=0)

    sigma_star = c*sigma_func(T) + 1.0/(dt)
    dz = Lz/Nz
    dr = Lr/Nr
    A,b = create_A(c*Dleft,c*Dright, c*Dtop,c*Dbottom, sigma, Nr, Nz, Lr, Lz, lower_z=lower_z,upper_z=upper_z, upper_r=upper_r);
    b += RHS
    return A\b
end #function

"One step of Equation 17"
function Estep(delta_t,T,T_prev,E,E_prev,E_l,D_func_arr,sigma_func_arr,Q_func_arr,group_bounds,EOS,Inv_EOS,Cv,rho,Nr,Nz,Lr,Lz,G; lower_z = "refl",upper_z="refl",upper_r="refl", LOUD=0)
    idt = 1.0/delta_t
    Jee = zeros(Nr*Nz)
    lagged_Emission = 0*Eg
    for g in 1:G
        Jee += sigma_func_arr[g](T).*DBgDT(group_bounds[g],group_bounds[g+1],T)
        lagged_Emission += reshape(sigma_func_arr[g](T),Nr*Nz).*(reshape(Bg(group_bounds[g],group_bounds[g+1],T),Nr*Nz) - E[:,g])*c
    end
    Jee *= c./reshape(Cv(T))
    Jee +=  reshape(rho,Nr,Nz)*idt

    Esteps = zeros(Nr*Nz,G)
    sigma_input = sigma_func_arr[g](T)
    e = EOS(T)
    e_prev = EOS(T_prev)
    for g in 1:G
        RHS = -c./reshape(Cv(T)).*reshape(sigma_func_arr[g](T).*DBgDT(group_bounds[g],group_bounds[g+1],T),Nr,Nz)./Jee
        RHS *= (reshape(rho,Nr,Nz)*idt*(e-e_prev) + reshape(lagged_Emission,Nr,Nz))
        RHS +=  reshape(sigma_func_arr[g](T),Nr,Nz).*reshape(E[:,g],Nr,Nz)*c #this group emission (2nd term on RHS)
        RHS +=  idt*reshape(E[:,g],Nr,Nz) #first term on RHS of 17
        RHS +=  reshape(Q_func(T),Nr,Nz) #source term on RHS of 17
        sigma_input *= (1-c./reshape(Cv(T)).*reshape(sigma_func_arr[g](T).*
                DBgDT(group_bounds[g],group_bounds[g+1],T),Nr,Nz)./Jee)
        sigma_input += idt

        #define Tleft
        tmp = ones(Nr,Nz)
        tmp[2:Nr,:] = 0.5*(Trect[1:(Nr-1),:] + Trect[2:Nr,:])
        sigma_left = sigma_func(t,tmp,Nr,Nz,Lr,Lz)
        Dleft = D_func_arr[g](t,tmp,Nr,Nz,Lr,Lz,Er,sigma_left)
        #define Tright
        tmp = ones(Nr,Nz)
        tmp[1:(Nr-1),:] = 0.5*(Trect[1:(Nr-1),:] + Trect[2:Nr,:])
        sigma_right = sigma_func(t,tmp,Nr,Nz,Lr,Lz)
        Dright = D_func_arr[g](t,tmp,Nr,Nz,Lr,Lz,Er,sigma_right)
        #define Ttop
        tmp = ones(Nr,Nz)
        tmp[:,1:(Nz-1)] = 0.5*(Trect[:,1:(Nz-1)] + Trect[:,2:Nz])
        sigma_top = sigma_func(t,tmp,Nr,Nz,Lr,Lz)
        Dtop = D_func_arr[g](t,tmp,Nr,Nz,Lr,Lz,Er,sigma_top)
        #define Tbottom
        tmp = ones(Nr,Nz)
        tmp[:,2:Nz] = 0.5*(Trect[:,1:(Nz-1)] + Trect[:,2:Nz])
        sigma_bottom = sigma_func(t,tmp,Nr,Nz,Lr,Lz)
        Dbottom = D_func_arr[g](t,tmp,Nr,Nz,Lr,Lz,Er,sigma_bottom)

        A,b = create_A(c*Dleft,c*Dright, c*Dtop,c*Dbottom, sigma_input, Nr, Nz, Lr, Lz, lower_z=lower_z,upper_z=upper_z, upper_r=upper_r);
        b += RHS

        Esteps[:,g] = A\b

    end

    return Esteps
end #function

max_stars_repo_name: DrRyanMc/HEDPRadDiffusion

max_issues_repo_name: DrRyanMc/HEDPRadDiffusion

max_forks_repo_name: DrRyanMc/HEDPRadDiffusion


text:
macro def_discrete_fit_sepp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    n = length(l)
    x0_expr = Expr[]
    fitted_params = :(setproperty!(sepp, $(l[1]), Optim.minimizer(res)[1]))

    h_it = Expr[]

    kw_expr_it = Expr[]

    for i in 1:n
        push!(x0_expr, :(θ[$(l[i])]))
        push!(kw_expr_it, :($(l[i]) => x[$i]))
        push!(h_it, :(getproperty(sepp, $(l[i]))))
    end

    for i in 2:n
        fitted_params = :($fitted_params ; $(:(setproperty!(sepp, $(l[i]), Optim.minimizer(res)[$i]))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, :x), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return discrete_negloglik(ts; final_kw...)
    end)

    h = :([$(h_it...)])

    x0_expr = :([$(x0_expr...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:discrete_fit!))), Expr(Symbol("::"), :sepp, :SEPPExpKern), m...),quote

        ts = sepp.data
        isnothing(ts) && error("No data in model, can't fit")

        θ = params(sepp)

        $to_min_expr

        lower = zeros($n)
        upper = fill(Inf, $n)

        res = optimize(to_min, lower, upper, $x0_expr)

        $fitted_params

        hess = ForwardDiff.hessian(to_min, $h)

        sepp.cov_mat = inv(hess)

        return res

    end)
end

#=
macro def_discrete_fit_sepp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    x = Symbol[]
    variables = :(@variable(model, $(Symbol("y1")) >= 0, start = θ[$(l[1])]))

    fitted_params = :(setproperty!(sepp, $(l[1]), value($(Symbol("y1")))))

    h_it = Expr[]

    kw_expr_it = Expr[]

    opt_var = Symbol[]

    for i in 1:length(l)
        push!(x, Symbol("x$i"))
        push!(opt_var, Symbol("y$i"))
        push!(kw_expr_it, :($(l[i]) => $(Symbol("x$i"))))
        push!(h_it, :(getproperty(sepp, $(l[i]))))
    end

    for i in 2:length(l)
        variables = :($variables ; @variable(model, $(Symbol("y$i")) >= 0, start = θ[$(l[i])]))

        fitted_params = :($fitted_params ; $(:(setproperty!(sepp, $(l[i]), value($(Symbol("y$i")))))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, x...), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return discrete_negloglik(ts; final_kw...)
    end)

    h = :([$(h_it...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:discrete_fit!))), Expr(Symbol("::"), :sepp, :SEPPExpKern), m...),quote

        ts = sepp.data
        isnothing(ts) && error("No data in model, can't fit")

        model = Model(Ipopt.Optimizer)
        θ = params(sepp)

        $to_min_expr

        JuMP.register(model, :to_min, length($params_to_fit) , to_min, autodiff=true)

        $variables

        @NLobjective(model, Min, to_min($(opt_var...)))

        optimize!(model)

        $fitted_params

        to_hess(θ) = to_min(θ...)

        hess = ForwardDiff.hessian(to_hess, $h)

        sepp.cov_mat = inv(hess)

        return objective_value(model)

    end)
end

=#

macro def_discrete_fit_sempp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    n = length(l)
    x0_expr = Expr[]
    fitted_params = :(setproperty!(sempp, $(l[1]), Optim.minimizer(res)[1]))

    h_it = Expr[]

    kw_expr_it = Expr[]

    for i in 1:n
        push!(x0_expr, :(θ[$(l[i])]))
        push!(kw_expr_it, :($(l[i]) => x[$i]))
        push!(h_it, :(getproperty(sempp, $(l[i]))))
    end

    for i in 2:n
        fitted_params = :($fitted_params ; $(:(setproperty!(sempp, $(l[i]), Optim.minimizer(res)[$i]))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, :x), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ], :δ => θ[:δ], :ξ => θ[:ξ], :β => θ[:β], :α => θ[:α], :κ => θ[:κ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return discrete_negloglik(mts, markdens, impact_func; final_kw...)
    end)

    h = :([$(h_it...)])

    x0_expr = :([$(x0_expr...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:discrete_fit!))), Expr(:parameters, Expr(:kw, :lower, zeros(n)), Expr(:kw, :upper, fill(Inf, n))), Expr(Symbol("::"), :sempp, :SEMPPExpKern), m...),quote

        mts = sempp.data
        isnothing(mts) && error("No data in model, can't fit")

        θ = params(sempp)
        markdens = θ[:markdens]
        impact_func = θ[:impact_function]

        $to_min_expr

        res = optimize(to_min, lower, upper, $x0_expr)

        $fitted_params

        hess = ForwardDiff.hessian(to_min, $h)

        sempp.cov_mat = inv(hess)

        return res

    end)
end

#=
macro def_discrete_fit_sempp(params_to_fit)
    l = (qn -> qn.value).(params_to_fit.args)
    x = Symbol[]
    variables = :(@variable(model, $(Symbol("y1")) >= 0, start = θ[$(l[1])]))
    constraints = :(@constraint(model, $(Symbol("y1")) >= 0))

    fitted_params = :(setproperty!(sempp, $(l[1]), value($(Symbol("y1")))))

    h_it = Expr[]

    kw_expr_it = Expr[]

    opt_var = Symbol[]

    for i in 1:length(l)
        push!(x, Symbol("x$i"))
        push!(opt_var, Symbol("y$i"))
        push!(kw_expr_it, :($(l[i]) => $(Symbol("x$i"))))
        push!(h_it, :(getproperty(sempp, $(l[i]))))
    end

    for i in 2:length(l)
        variables = :($variables ; @variable(model, $(Symbol("y$i")) >= 0, start = θ[$(l[i])]))

        constraints = :($constraints ;  @constraint(model, $(Symbol("y$i")) >= 0))

        fitted_params = :($fitted_params ; $(:(setproperty!(sempp, $(l[i]), value($(Symbol("y$i")))))))
    end

    kw_expr = :(kw = Dict($(kw_expr_it...)))

    to_min_expr = Expr(:function, Expr(:call, :to_min, x...), quote
        default_kw = Dict{Symbol, Real}(:μ => θ[:μ], :ϕ => θ[:ϕ], :γ => θ[:γ], :δ => θ[:δ], :ξ => θ[:ξ], :β => θ[:β], :α => θ[:α], :κ => θ[:κ])
        $kw_expr
        final_kw = merge(default_kw, kw)
        return discrete_negloglik(mts, markdens, impact_func; final_kw...)
    end)

    h = :([$(h_it...)])

    m = esc.(l)

    return Expr(:function, Expr(:call, :($(esc(:discrete_fit!))), Expr(:parameters, Expr(:kw, :tol, 10^(-8)), Expr(:kw, :max_iter, 1000)), Expr(Symbol("::"), :sempp, :SEMPPExpKern), m...),quote

        mts = sempp.data
        isnothing(mts) && error("No data in model, can't fit")

        model = Model(Ipopt.Optimizer)
        set_optimizer_attribute(model, "tol", tol)
        set_optimizer_attribute(model, "max_iter", max_iter)

        θ = params(sempp)
        markdens = θ[:markdens]
        impact_func = θ[:impact_function]

        $to_min_expr

        JuMP.register(model, :to_min, length($params_to_fit) , to_min, autodiff=true)

        $variables

        @NLobjective(model, Min, to_min($(opt_var...)))

        optimize!(model)

        $fitted_params

        to_hess(θ) = to_min(θ...)

        hess = ForwardDiff.hessian(to_hess, $h)

        sempp.cov_mat = inv(hess)

        return objective_value(model)

    end)
end
=#


text:
isfsal(alg::OrdinaryDiffEqAlgorithm) = false
isfsal(alg::DP5) = true
isfsal(alg::DP5Threaded) = true
isfsal(alg::DP8) = true
isfsal(alg::BS3) = true
isfsal(alg::BS5) = true
isfsal(alg::Tsit5) = true
isfsal(alg::Vern6) = true
isfsal(alg::Rosenbrock23) = true
isfsal(alg::Rosenbrock32) = true
isfsal(alg::Euler) = true
isfsal(alg::Midpoint) = true
isfsal(alg::RK4) = true
isfsal(alg::Feagin10) = true
isfsal(alg::Feagin12) = true
isfsal(alg::Feagin14) = true
isfsal(alg::TanYam7) = true
isfsal(alg::TsitPap8) = true
isfsal(alg::Trapezoid) = true
isfsal(alg::ImplicitEuler) = true
isfsal(alg::ExplicitRK) = true
isfsal{MType,VType,fsal}(tab::ExplicitRKTableau{MType,VType,fsal}) = fsal
#isfsal(tab::ImplicitRKTableau) = false
isfsal(alg::CompositeAlgorithm) = true # Every algorithm is assumed FSAL. Good assumption?

isimplicit(alg::OrdinaryDiffEqAlgorithm) = false
isimplicit(alg::ImplicitEuler) = true
isimplicit(alg::Trapezoid) = true

isdtchangeable(alg::OrdinaryDiffEqAlgorithm) = true

ismultistep(alg::OrdinaryDiffEqAlgorithm) = false

isadaptive(alg::OrdinaryDiffEqAlgorithm) = false
isadaptive(alg::OrdinaryDiffEqAdaptiveAlgorithm) = true
isadaptive(alg::OrdinaryDiffEqCompositeAlgorithm) = isadaptive(alg.algs[1])

qmin_default(alg::OrdinaryDiffEqAlgorithm) = 1//5
qmin_default(alg::DP8) = 1//3

qmax_default(alg::OrdinaryDiffEqAlgorithm) = 10
qmax_default(alg::DP8) = 6

get_chunksize(alg::OrdinaryDiffEqAlgorithm) = error("This algorithm does not have a chunk size defined.")
get_chunksize{CS,AD}(alg::ImplicitEuler{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Trapezoid{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Rosenbrock23{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Rosenbrock32{CS,AD}) = CS

alg_extrapolates(alg::OrdinaryDiffEqAlgorithm) = false
alg_extrapolates(alg::ImplicitEuler) = true
alg_extrapolates(alg::Trapezoid) = true

alg_autodiff(alg::OrdinaryDiffEqAlgorithm) = error("This algorithm does not have an autodifferentiation option defined.")
alg_autodiff{CS,AD}(alg::ImplicitEuler{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Trapezoid{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Rosenbrock23{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Rosenbrock32{CS,AD}) = AD

alg_order(alg::OrdinaryDiffEqAlgorithm) = error("Order is not defined for this algorithm")
alg_adaptive_order(alg::OrdinaryDiffEqAdaptiveAlgorithm) = error("Algorithm is adaptive with no order")

alg_order(alg::Discrete) = 0
alg_order(alg::Euler) = 1
alg_order(alg::Midpoint) = 2
alg_order(alg::RK4) = 4
alg_order(alg::ExplicitRK) = alg.tableau.order
alg_order(alg::BS3) = 3
alg_order(alg::BS5) = 5
alg_order(alg::DP5) = 5
alg_order(alg::DP5Threaded) = 5
alg_order(alg::Tsit5) = 5
alg_order(alg::DP8) = 8
alg_order(alg::Vern6) = 6
alg_order(alg::Vern7) = 7
alg_order(alg::Vern8) = 8
alg_order(alg::Vern9) = 9
alg_order(alg::TanYam7) = 7
alg_order(alg::TsitPap8) = 8
alg_order(alg::ImplicitEuler) = 1
alg_order(alg::Trapezoid) = 2
alg_order(alg::Rosenbrock23) = 2
alg_order(alg::Rosenbrock32) = 3
alg_order(alg::Feagin10) = 10
alg_order(alg::Feagin12) = 12
alg_order(alg::Feagin14) = 14

alg_order(alg::CompositeAlgorithm) = alg_order(alg.algs[1])

alg_adaptive_order(alg::ExplicitRK) = alg.tableau.adaptiveorder
alg_adaptive_order(alg::BS3) = 2
alg_adaptive_order(alg::BS5) = 4
alg_adaptive_order(alg::DP5) = 4
alg_adaptive_order(alg::DP5Threaded) = 4
alg_adaptive_order(alg::Tsit5) = 4
alg_adaptive_order(alg::DP8) = 6
alg_adaptive_order(alg::Vern6) = 5
alg_adaptive_order(alg::Vern7) = 6
alg_adaptive_order(alg::Vern8) = 7
alg_adaptive_order(alg::Vern9) = 8
alg_adaptive_order(alg::TanYam7) = 6
alg_adaptive_order(alg::TsitPap8) = 7
alg_adaptive_order(alg::Rosenbrock23) = 3
alg_adaptive_order(alg::Rosenbrock32) = 2
alg_adaptive_order(alg::Feagin10) = 8
alg_adaptive_order(alg::Feagin12) = 10
alg_adaptive_order(alg::Feagin14) = 12

beta2_default(alg::OrdinaryDiffEqAlgorithm) = 2//(5alg_order(alg))
beta2_default(alg::Discrete) = 0
beta2_default(alg::DP8) = 0//1
beta2_default(alg::DP5) = 4//100
beta2_default(alg::DP5Threaded) = 4//100

beta1_default(alg::OrdinaryDiffEqAlgorithm,beta2) = 7//(10alg_order(alg))
beta1_default(alg::Discrete,beta2) = 0
beta1_default(alg::DP8,beta2) = typeof(beta2)(1//alg_order(alg)) - beta2/5
beta1_default(alg::DP5,beta2) = typeof(beta2)(1//alg_order(alg)) - 3beta2/4

discrete_apply_map{apply_map,scale_by_time}(alg::Discrete{apply_map,scale_by_time}) = apply_map
discrete_scale_by_time{apply_map,scale_by_time}(alg::Discrete{apply_map,scale_by_time}) = scale_by_time