text:
export sheets_client

# The Google Sheets API has access limits.
# The default limits are 100 requests per 100 seconds per user.
# Limits for reads and writes are tracked separately.
# Higher limits are available in the Google Cloud console.
# https://developers.google.com/sheets/api/limits
#
# Reality seems to be more complex.
# The tokens per second seems to match the described limit.
# If the max tokens are too large, the burstiness triggers undocumented
# limits, without clearly indicating that burstiness is the problem.
#
# These defaults seem to work.
default_rate_limiter_tokens_per_sec = 0.9
default_rate_limiter_max_tokens = 4
default_rate_limiter_read = TokenBucketRateLimiter(default_rate_limiter_tokens_per_sec, default_rate_limiter_max_tokens, default_rate_limiter_max_tokens)
default_rate_limiter_write = TokenBucketRateLimiter(default_rate_limiter_tokens_per_sec, default_rate_limiter_max_tokens, default_rate_limiter_max_tokens)

"""
    update_default_rate_limiter(rate_limiter_tokens_per_sec::Float64; rate_limiter_max_tokens::Float64=5)

Update the default rate limiter.

# Arguments
- `rate_limiter_tokens_per_sec::Float64`: rate limiter number of tokens per second.
- `rate_limiter_max_tokens::Float64=5`: rate limiter max number of tokens.
"""
function update_default_rate_limiter(rate_limiter_tokens_per_sec::Float64; rate_limiter_max_tokens::Float64=5)
    global default_rate_limiter_read = TokenBucketRateLimiter(rate_limiter_tokens_per_sec, rate_limiter_max_tokens, rate_limiter_max_tokens)
    global default_rate_limiter_write = TokenBucketRateLimiter(rate_limiter_tokens_per_sec, rate_limiter_max_tokens, rate_limiter_max_tokens)
end

"""
    update_default_rate_limiter(rate_limiter_read::AbstractRateLimiter, rate_limiter_write::AbstractRateLimiter)

Update the default rate limiter.

# Arguments
- `rate_limiter_read::AbstractRateLimiter`: rate limiter for reading.
- `rate_limiter_write::AbstractRateLimiter`: rate limiter for writing.
"""
function update_default_rate_limiter(rate_limiter_read::AbstractRateLimiter, rate_limiter_write::AbstractRateLimiter)
    global default_rate_limiter_read = rate_limiter_read
    global default_rate_limiter_write = rate_limiter_write
end

"""
Print details on a python exception.
"""
macro _print_python_exception(ex)
    # MacroTools.@q is used instead of quote so that the returned stacktrace
    # has line numbers from the calling function and not the macro.
    return esc(MacroTools.@q begin
        try
            $ex
        catch e
            if hasfield(typeof(e), :traceback)
                println("Python error:")
                println(e)
                println("Python stacktrace:")
                tb = pyimport("traceback")
                tb.print_exception(e.traceback)
                tb.print_tb(e.traceback)
            end
            rethrow(e)
        end
    end)
end

"""
Directory containing configuration files.
"""
config_dir = joinpath(homedir(),".julia/config/google_sheets/")

const _permission_urls = Dict(
    AUTH_SCOPE_READONLY => "https://www.googleapis.com/auth/spreadsheets.readonly",
    AUTH_SCOPE_READWRITE => "https://www.googleapis.com/auth/spreadsheets",
)

"""
Maps authorization scopes to the appropriate permission URLs.
"""
function _scope_urls(scopes::AuthScope)::Array{String,1}
    return [_permission_urls[scopes]]
end

"""
Maps authorization scopes to the appropriate permission URLs.
"""
function _scope_urls(scopes::Array{AuthScope,1})::Array{String,1}
    return [_permission_urls[scope] for scope in scopes]
end

"""
    credentials_file()::String

Gets the credentials file needed to log into Google.
The file is loaded from the GOOGLE_SHEETS_CREDENTIALS environment variable
if it is present; otherwise, it is loaded from the configuration directory,
which defaults to ~/.julia/config/google_sheets/.

See the python quick start reference for a link to generate credentials.
https://developers.google.com/sheets/api/quickstart/python

See the API credentials page to create or revoke credentials.
https://console.developers.google.com/apis/credentials
"""
function credentials_file()::String
    if haskey(ENV, "GOOGLE_SHEETS_CREDENTIALS")
        return ENV["GOOGLE_SHEETS_CREDENTIALS"]
    end

    file = "credentials.json"
    return joinpath(config_dir, file)
end

"""
Gets a cached token file which allows access to Google Sheets with
specific authorization scopes.
"""
function _token_file(scopes::AuthScope)::String
    return _token_file([scopes])
end

"""
Gets a cached token file which allows access to Google Sheets with
specific authorization scopes.
"""
function _token_file(scopes::Array{AuthScope,1})::String
    s = sort(unique(copy(scopes)))

    id = 0

    for i in 1:length(s)
        id += Int(s[i]) * 10^i
    end

    file = "google_sheets_token.$id.pickle"
    return joinpath(config_dir, file)
end

"""
    sheets_client(scopes::Union{AuthScope,Array{AuthScope,1}};
        rate_limiter_read::AbstractRateLimiter=default_rate_limiter_read,
        rate_limiter_write::AbstractRateLimiter=default_rate_limiter_write)::GoogleSheetsClient

Creates a client for accessing Google Sheets.

# Arguments
- `scopes::Union{AuthScope,Array{AuthScope,1}}`: authorization scopes.
- `rate_limiter_read::AbstractRateLimiter=default_rate_limiter_read`: rate limiter for reading.
- `rate_limiter_write::AbstractRateLimiter=default_rate_limiter_write`: rate limiter for writing.
"""
function sheets_client(scopes::Union{AuthScope,Array{AuthScope,1}};
        rate_limiter_read::AbstractRateLimiter=default_rate_limiter_read,
        rate_limiter_write::AbstractRateLimiter=default_rate_limiter_write)::GoogleSheetsClient

    pickle = pyimport("pickle")
    os_path = pyimport("os.path")
    build = pyimport("googleapiclient.discovery").build
    InstalledAppFlow = pyimport("google_auth_oauthlib.flow").InstalledAppFlow
    Request = pyimport("google.auth.transport.requests").Request
    open = pybuiltin("open")

    credentialsFile = credentials_file()
    tokenFile = _token_file(scopes)
    scopeUrls = _scope_urls(scopes)

    @_print_python_exception begin
        creds = nothing

        # The file token.pickle stores the user's access and refresh tokens, and is
        # created automatically when the authorization flow completes for the first
        # time.
        if os_path.exists(tokenFile)
            @pywith open(tokenFile, "rb") as token begin
                creds = pickle.load(token)
            end
        end

        # If there are no (valid) credentials available, let the user log in.
        if isnothing(creds) || !creds.valid
            if !isnothing(creds) && creds.expired && !isnothing(creds.refresh_token)
                creds.refresh(Request())
            else
                if !isfile(credentialsFile)
                    error("Credentials file not found: $credentialsFile")
                end

                flow = InstalledAppFlow.from_client_secrets_file(credentialsFile, scopeUrls)
                creds = flow.run_local_server(port=0)
            end

            # Save the credentials for the next run
            @pywith open(tokenFile, "wb") as token begin
                pickle.dump(creds, token)
            end
        end

        return GoogleSheetsClient(build("sheets", "v4", credentials=creds), rate_limiter_read, rate_limiter_write)
    end
end

"""
    gsheet_api_spreadsheet_get(client::GoogleSheetsClient; kwargs...)

GoogleSheets spreadsheet-get API.

# Arguments
- `client::GoogleSheetsClient`: client
- `kwargs`: keyword arguments.
"""
function gsheet_api_spreadsheet_get(client::GoogleSheetsClient; kwargs...)
    @rate_limit client.rate_limiter_read 1 @_print_python_exception begin
        return client.client.spreadsheets().get(;kwargs...).execute()
    end
end

"""
    gsheet_api_sheet_get(client::GoogleSheetsClient; kwargs...)

GoogleSheets sheet-get API.

# Arguments
- `client::GoogleSheetsClient`: client
- `kwargs`: keyword arguments.
"""
function gsheet_api_sheet_get(client::GoogleSheetsClient; kwargs...)
    @rate_limit client.rate_limiter_read 1 @_print_python_exception begin
        return client.client.spreadsheets().values().get(;kwargs...).execute()
    end
end

"""
    gsheet_api_sheet_batchget(client::GoogleSheetsClient; kwargs...)

GoogleSheets sheet-batchGet API.

# Arguments
- `client::GoogleSheetsClient`: client
- `kwargs`: keyword arguments.
"""
function gsheet_api_sheet_batchget(client::GoogleSheetsClient; kwargs...)
    @rate_limit client.rate_limiter_read 1 @_print_python_exception begin
        return client.client.spreadsheets().values().batchGet(;kwargs...).execute()
    end
end

"""
    gsheet_api_sheet_update(client::GoogleSheetsClient; kwargs...)

GoogleSheets sheet-update API.

# Arguments
- `client::GoogleSheetsClient`: client
- `kwargs`: keyword arguments.
"""
function gsheet_api_sheet_update(client::GoogleSheetsClient; kwargs...)
    @rate_limit client.rate_limiter_write 1 @_print_python_exception begin
        return client.client.spreadsheets().values().update(;kwargs...).execute()
    end
end

"""
    gsheet_api_speadsheet_batchupdate(client::GoogleSheetsClient; kwargs...)

GoogleSheets spreadsheet-batchUpdate API.

# Arguments
- `client::GoogleSheetsClient`: client
- `kwargs`: keyword arguments.
"""
function gsheet_api_speadsheet_batchupdate(client::GoogleSheetsClient; kwargs...)
    @rate_limit client.rate_limiter_write 1 @_print_python_exception begin
        return client.client.spreadsheets().batchUpdate(;kwargs...).execute()
    end
end

ext: jl


text:
using PackageCompiler

create_sysimage([:Gridap,:PartitionedArrays,:GridapMHD],
  sysimage_path=joinpath(@__DIR__,"..","GridapMHD.so"),
  precompile_execution_file=joinpath(@__DIR__,"warmup.jl"))


text:
# Autogenerated wrapper script for gmmlib_jll for x86_64-linux-gnu
export libigdgmm

JLLWrappers.@generate_wrapper_header("gmmlib")
JLLWrappers.@declare_library_product(libigdgmm, "libigdgmm.so.11")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libigdgmm,
        "lib/libigdgmm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
Clang.isdecl(x) = false
function _convert(ctx::ConverterContext, decl::CLTypedefDecl)
	name = _convertType(ctx, decl)
	(o, c, typ) = _convertTypedefType(ctx, underlying_type(decl))
	typ = isdecl(typ) || typedecl(typ) isa CLNoDeclFound ? typ : typedecl(typ)
	kind = :atcompile

	# if typ is the definition and it was just converted then we need to remove that conversion in favor of this one
	if (typ isa CLEnumDecl || typ isa CLUnionDecl || typ isa CLStructDecl) && !isempty(ctx.converted) && getdef(typ) == last(ctx.converted).decl
		pop!(ctx.converted)
		typ = typ isa CLEnumDecl ? _convertEnum(ctx, typ) : _convertAggregate(ctx, typ)
	else
		typ = _convertType(ctx, typ)
		ctx.typedefs[name] = "$(o)$(typ)$(c)"
		kind = :atcompile_typedefs
	end

	_export(ctx, name)
	push!(ctx.converted, JuliaizedC(
		decl,
		"CBinding.@ctypedef $(name) $(o)$(typ)$(c)",
		kind,
	))
end

function _convert(ctx::ConverterContext, decl::CLEnumDecl)
	expr = _convertEnum(ctx, decl)
	push!(ctx.converted, JuliaizedC(
		decl,
		expr,
		:atcompile,
	))
end

function _convert(ctx::ConverterContext, decl::Union{CLUnionDecl, CLStructDecl})
	expr = _convertAggregate(ctx, decl)
	isnothing(expr) || push!(ctx.converted, JuliaizedC(
		decl,
		expr,
		:atcompile,
	))
end

function _convertEnum(ctx, decl::CLEnumDecl; indent::Int = 0)
	name = _convertType(ctx, decl)

	vals = []
	for c in children(decl)
		vname = _convertName(ctx, c)
		vtype = _convertType(ctx, type(c))
		vval  = value(c)
		_export(ctx, vname)
		push!(vals, "$(vname) = $(vtype)($(vval))")
	end

	tabs = repeat('\t', indent)
	body = "{"*(isempty(vals) ? "" : "\n$(tabs)\t$(join(vals, "\n$(tabs)\t"))")*"\n$(tabs)}"

	_export(ctx, name)
	return "CBinding.@cenum $(name) $(body)"
end

function _convertAggregate(ctx, decl::Union{CLUnionDecl, CLStructDecl}; indent::Int = 0)
	kind = decl isa CLUnionDecl ? "CBinding.@cunion" : "CBinding.@cstruct"
	name = _convertType(ctx, decl)
	packing = nothing

	if getdef(decl) isa CLInvalidFile  # this is an opaque type, don't convert, use _convertType(ctx, decl) instead
		return nothing
	elseif getdef(decl) != decl  # this is a forward decl, so wait for actual decl to convert
		return nothing
	else
		fields = []
		lastAgg = nothing
		for c in children(decl)
			#@info c
			if c isa CLFieldDecl
				fname = _convertName(ctx, c)
				ftype = type(c)
				if ftype isa CLElaborated
					ftype = typedecl(ftype)
					if lastAgg == ftype && !isnothing(last(last(fields)))
						push!(last(last(fields)), fname)
					else
						ftypename = _convertType(ctx, ftype)
						push!(fields, ftypename => fname)
						lastAgg = nothing
					end
				else
					if isbit(c)
						fname = "($(fname):$(bitwidth(c)))"
						ftypename = _convertType(ctx, ftype)
						push!(fields, ftypename => fname)
						lastAgg = nothing
					else
						(augment, atype) = _convertAugmentedType(ctx, ftype)
						if lastAgg == atype && !isnothing(last(last(fields))) && (atype isa CLEnumDecl || atype isa CLUnionDecl || atype isa CLStructDecl)
							push!(last(last(fields)), fname*"::"*augment)
						else
							ftypename = _convertType(ctx, ftype)
							push!(fields, ftypename => fname)
							lastAgg = nothing
						end
					end
				end
			elseif c isa CLEnumDecl || c isa CLUnionDecl || c isa CLStructDecl
				ftype = c isa CLEnumDecl ? _convertEnum(ctx, c, indent = indent+1) : _convertAggregate(ctx, c, indent = indent+1)
				ftype = isnothing(ftype) ? _convertType(ctx, c) : ftype
				push!(fields, ftype => (isempty(spelling(c)) ? [] : nothing))
				lastAgg = c
			elseif c isa CLPackedAttr && isnothing(packing)
				packing = spelling(c)
			else
				error("Unable to convert nested aggregate type at $(Interval(c)) `$(c)`")
			end
		end

		fields = map(fields) do (ftype, fname)
			if isnothing(fname) || (fname isa Array && isempty(fname))
				return ftype
			elseif fname isa Array && !isempty(fname)
				return "($(join(fname, ", ")))::$(ftype)"
			else
				return "$(fname)::$(ftype)"
			end
		end

		tabs = repeat('\t', indent)
		body = "{"*(isempty(fields) ? "" : "\n$(tabs)\t$(join(fields, "\n$(tabs)\t"))")*"\n$(tabs)}"
		body = isnothing(packing) ? body : "$(body) __$(packing)__"
	end

	_export(ctx, name)
	return "$(kind) $(name) $(body)"
end

_convertTypedefType(ctx, decl) = ("", "", decl)
_convertTypedefType(ctx, decl::CLTypedef) = typedecl(decl) isa CLNoDeclFound ? ("", "", decl) : _convertTypedefType(ctx, typedecl(decl))
function _convertTypedefType(ctx, decl::CLConstantArray)
	(o, c, typ) = _convertTypedefType(ctx, element_type(decl))
	return ("(CBinding.@carray ("*o,c*")[$(element_num(decl))])", typ)
end

_convertAugmentedType(ctx, decl) = ("{}", decl)
 _convertAugmentedType(ctx, decl::CLElaborated) = _convertAugmentedType(ctx, typedecl(decl))
function _convertAugmentedType(ctx, decl::CLConstantArray)
	(augment, typ) = _convertAugmentedType(ctx, element_type(decl))
	return ("(CBinding.@carray ($(augment))[$(element_num(decl))])", typ)
end
todo"determine the best way to represent incomplete array types, as pointers, or as empty arrays"
function _convertAugmentedType(ctx, decl::CLIncompleteArray)
	(augment, typ) = _convertAugmentedType(ctx, element_type(decl))
	augment = augment == "{}" ? "" : augment
	return ("Base.Ptr{$(augment)}", typ)
end
function _convertAugmentedType(ctx, decl::CLPointer)
	(augment, typ) = _convertAugmentedType(ctx, pointee_type(decl))
	augment = augment == "{}" ? "" : augment
	return ("Base.Ptr{$(augment)}", typ)
end

_convertType(ctx, typ::CLChar_S) = "Base.Cchar"
_convertType(ctx, typ::CLChar_U) = "Base.Cuchar"
_convertType(ctx, typ::CLSChar) = "Base.Cchar"
_convertType(ctx, typ::CLUChar) = "Base.Cuchar"

_convertType(ctx, typ::CLShort) = "Base.Cshort"
_convertType(ctx, typ::CLInt) = "Base.Cint"
_convertType(ctx, typ::CLLong) = "Base.Clong"
_convertType(ctx, typ::CLLongLong) = "Base.Clonglong"

_convertType(ctx, typ::CLUShort) = "Base.Cushort"
_convertType(ctx, typ::CLUInt) = "Base.Cuint"
_convertType(ctx, typ::CLULong) = "Base.Culong"
_convertType(ctx, typ::CLULongLong) = "Base.Culonglong"

_convertType(ctx, typ::CLFloat) = "Base.Cfloat"
_convertType(ctx, typ::CLDouble) = "Base.Cdouble"
_convertType(ctx, typ::CLLongDouble) = "CBinding.Clongdouble"
_convertType(ctx, typ::CLComplex) = "Base.Complex{$(_convertType(ctx, element_type(typ)))}"

_convertType(ctx, typ::CLBool) = "Base.Cbool"
_convertType(ctx, typ::CLVoid) = "Base.Cvoid"
_convertType(ctx, typ::CLFirstBuiltin) = "Base.Cvoid"   # TODO:  this decision needs to be evaluated!
_convertType(ctx, typ::CLVector) = error("$(typ) $(element_num(typ))")
function _convertType(ctx, typ::Union{CLTypedef, CLTypedefDecl})
	name = typ isa CLTypedefDecl || typedecl(typ) isa CLNoDeclFound ? _convertName(ctx, typ) : _convertType(ctx, typedecl(typ))
	return haskey(ctx.typedefs, name) ? ctx.typedefs[name] : name
end
_convertType(ctx, typ::CLPointer) = "Base.Ptr{$(_convertType(ctx, pointee_type(typ)))}"
_convertType(ctx, typ::CLConstantArray) = "(CBinding.@carray ($(_convertType(ctx, element_type(typ))))[$(element_num(typ))])"
_convertType(ctx, typ::CLIncompleteArray) = "Base.Ptr{$(_convertType(ctx, element_type(typ)))}"
_convertType(ctx, typ::CLElaborated) = _convertType(ctx, typedecl(typ))
_convertType(ctx, typ::Union{CLEnumDecl, CLUnionDecl, CLStructDecl}) = _convertName(ctx, typ)

Clang.calling_conv(t::CLType) = Clang.calling_conv(t.type)
Clang.argnum(t::CLType) = argnum(t.type)
Clang.argtype(t::CLType, i) = CLType(argtype(t.type, UInt(i)))
Clang.result_type(t::CLType) = CLType(result_type(t.type))
function _convertType(ctx, typ::CLUnexposed)
	num = argnum(typ)
	if num >= 0   # this is a function signature
		ret = _convertType(ctx, result_type(typ))
		args = map(i -> _convertType(ctx, argtype(typ, i-1)), 1:num)
		isvariadic(typ) && push!(args, "Base.Vararg")
		return "CBinding.Cfunction{$(ret), Base.Tuple{$(join(args, ", "))}}"
	else
		error("A non-function `unexposed` type is not yet supported")
	end
end

_convertName(ctx, decl::CLCursor) = _convertName(ctx, spelling(decl))
_convertName(ctx, name::String) = (name in Clang.RESERVED_WORDS) ? name*"_" : name


text:
abstract type AbstractGridWorld end

#####
##### game logic
#####

reset!(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")
act!(env::AbstractGridWorld, action) = error("Method not implemented for $(typeof(env))")

#####
##### optional (pretty printing, playing, etc...)
#####

get_pretty_tile_map(env::AbstractGridWorld, position::CartesianIndex{2}) = error("Method not implemented for $(typeof(env))")
get_pretty_sub_tile_map(env::AbstractGridWorld, window_size, position::CartesianIndex{2}) = error("Method not implemented for $(typeof(env))")
get_action_keys(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")
get_action_names(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")
get_object_names(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")
get_height(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")
get_width(env::AbstractGridWorld) = error("Method not implemented for $(typeof(env))")

function get_pretty_tile_map(env::AbstractGridWorld)
    height = get_height(env)
    width = get_width(env)

    str = ""

    for i in 1:height
        for j in 1:width
            str = str * get_pretty_tile_map(env, CartesianIndex(i, j))
        end
        if i < height
            str = str * "\n"
        end
    end

    return str
end

function get_pretty_sub_tile_map(env::AbstractGridWorld, window_size)
    height, width = window_size

    str = ""

    for i in 1:height
        for j in 1:width
            str = str * get_pretty_sub_tile_map(env, window_size, CartesianIndex(i, j))
        end
        if i < height
            str = str * "\n"
        end
    end

    return str
end

#####
##### sub tile map
#####

function get_window_size(env::AbstractGridWorld)
    height = get_height(env)
    width = get_width(env)
    return (2 * (height ÷ 4) + 1, 2 * (width ÷ 4) + 1)
end

function get_window_region((i, j), (m, n))
    temp1 = m ÷ 2
    temp2 = i - temp1
    temp3 = n ÷ 2
    temp4 = j - temp3
    return CartesianIndices((temp2 : temp2 + m - 1, temp4 : temp4 + n - 1))
end

function get_window_region((i, j), (m, n), direction)
    if direction == RIGHT
        temp1 = n - 1
        temp2 = temp1 ÷ 2
        temp3 = i - temp2
        return CartesianIndices((temp3 : temp3 + temp1, j : j + m - 1))
    elseif direction == UP
        temp1 = n - 1
        temp2 = temp1 ÷ 2
        temp3 = j - temp2
        return CartesianIndices((i - m + 1 : i, temp3 : temp3 + temp1))
    elseif direction == LEFT
        temp1 = n - 1
        temp2 = temp1 ÷ 2
        temp3 = i - temp2
        return CartesianIndices((temp3 : temp3 + temp1, j - m + 1 : j))
    else
        temp1 = n - 1
        temp2 = temp1 ÷ 2
        temp3 = j - temp2
        return CartesianIndices((i : i + m - 1, temp3 : temp3 + temp1))
    end
end

function map_index((i,j), (m, n), direction)
    if direction == RIGHT
        return (j, n-i+1)
    elseif direction == UP
        return (m-i+1, n-j+1)
    elseif direction == LEFT
        return (m-j+1, i)
    else
        return (i,j)
    end
end

function get_sub_tile_map(tile_map, position, window_size)
    num_objects = size(tile_map, 1)
    sub_tile_map = falses(num_objects, window_size...)
    get_sub_tile_map!(sub_tile_map, tile_map, position, window_size)
    return sub_tile_map
end

function get_sub_tile_map(tile_map, position, window_size, direction)
    num_objects = size(tile_map, 1)
    sub_tile_map = falses(num_objects, window_size...)
    get_sub_tile_map!(sub_tile_map, tile_map, position, window_size, direction)
    return sub_tile_map
end

function get_sub_tile_map!(sub_tile_map, tile_map, position, window_size)
    _, height, width = size(tile_map)

    window_region = get_window_region(position.I, window_size)

    valid_region = CartesianIndices((1 : height, 1 : width))

    @views for key in keys(window_region)
        pos = window_region[key]
        if pos in valid_region
            sub_tile_map[:, key] .= tile_map[:, pos]
        end
    end

    return nothing
end

function get_sub_tile_map!(sub_tile_map, tile_map, position, window_size, direction)
    _, height, width = size(tile_map)

    window_region = get_window_region(position.I, window_size, direction)

    valid_region = CartesianIndices((1 : height, 1 : width))

    @views for key in keys(window_region)
        pos = window_region[key]
        if pos in valid_region
            sub_tile_map[:, map_index(key.I, window_size, direction)...] .= tile_map[:, pos]
        end
    end

    return nothing
end

#####
##### sampling tile map positions
#####

function sample_empty_position(rng, tile_map, region, max_tries)
    position = rand(rng, region)

    for i in 1:max_tries
        if any(@view tile_map[:, position])
            position = rand(rng, region)
        else
            return position
        end
    end

    @warn "Could not sample an empty position in max_tries = $(max_tries). Returning non-empty position: $(position)"

    return position
end

function sample_empty_position(rng, tile_map, region)
    max_tries = 1024 * length(region)
    position = sample_empty_position(rng, tile_map, region, max_tries)
    return position
end

function sample_empty_position(rng, tile_map, max_tries::Integer)
    _, height, width = size(tile_map)
    region = CartesianIndices((1 : height, 1 : width))
    position = sample_empty_position(rng, tile_map, region, max_tries)
    return position
end

function sample_empty_position(rng, tile_map)
    _, height, width = size(tile_map)
    region = CartesianIndices((1 : height, 1 : width))
    max_tries = 1024 * height * width
    position = sample_empty_position(rng, tile_map, region, max_tries)
    return position
end

function sample_two_positions_without_replacement(rng, region)
    position1 = rand(rng, region)
    position2 = rand(rng, region)

    while position1 == position2
        position2 = rand(rng, region)
    end

    return position1, position2
end


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

"""
    S = fourierSeries(T, aₖ, kInds)

Create a 'AMFMmodel' by specifiying components drawn from a Fourier Series
  where 'T' is the fundemental period, 'aₖ' is a function specifying the Fourier
  series coefficients, and 'kInds' is a vector containing values of `k` to
  consider in the partial sum.

# Examples
```@example
using ISA, Plots
T = 1.0
aₖ(k) = 1.0
kInds = collect(-25:25)
z = fourierSeries(T, aₖ, kInds)
plot(z; timeaxis=-1.0:0.001:1.0)
```
"""
function fourierSeries(T::Real, aₖ::Function, kInds::Array{Int,1}=Vector(-100:100))
  function getIA(k)
    function a(t::Real)
      ia = abs(aₖ(k))
      return ia::Real
    end
    return a
  end
  function getIF(k)
      function ω(t::Real)
        omega = 2*π*k/T
        return omega::Real
      end
    return ω
  end
  function getPh(k)
    ph = angle(aₖ(k))
    return ph
  end
  S = [ AMFMcomp( getIA(k), getIF(k), getPh(k) ) for k ∈ kInds]
  return compSet(S)
end

# CONSTRUCTIONS
fourierSeries(T::Real, aₖ::Function, kInds::StepRangeLen) = fourierSeries(T, aₖ, collect(kInds))
fourierSeries(T::Real, aₖ::Function, kInds::UnitRange) = fourierSeries(T, aₖ, collect(kInds))


text:
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials, Test
	const PDMP = PiecewiseDeterministicMarkovProcesses

function AnalyticalSampleCHV(xc0, xd0, ti, nj::Int64)
	xch = [xc0[1]]
	xdh = [xd0[1]]
	th  = [ti]
	t = ti
	while length(th)<nj
		xc = xch[end]
		xd = xdh[end]
		S = -log(rand())
		if mod(xd,2) == 0
			t += 1/10*log(1+10*S/xc)
			push!(xch,xc + 10 * S )
		else
			t += 1/(3xc)*(exp(3S)-1)
			push!(xch,xc * exp(-3S) )
		end
		push!(xdh,xd + 1 )
		push!(th,t)

		S = -log(rand())
	end
	return th, xch, xdh
end

function AnalyticalSampleRejection(xc0,xd0,ti,nj::Int64; verbose = false)
	verbose && printstyled(color=:red,"--> Start analytical method\n")
	xch = [xc0[1]]
	xdh = [xd0[1]]
	th  = [ti]
	t = ti
	xc = xc0[1]
	njumps = 1
	rt = zeros(2)
	lambda_star = R!(rt,xc0,xd0,ti,Float64[],true)[2]
	rate		= R!(rt,xc0,xd0,ti,Float64[],true)[1]
	S = -log(rand()) / lambda_star
	while njumps < nj
		xd  = [xdh[end] ,1]

		t += S
		if mod(xd[1],2) == 0
			xc = xc * exp(10 * S)#xc + S
		else
			xc = xc / (3 * S * xc + 1)#xc * exp(-10S)
		end
		verbose && println("--> S = $S, t = $t, xc = $xc, xd = $(xd[1]), λ_* = ", lambda_star)
		#reject?
		lambda_star = R!(rt,[xc],xd,ti,Float64[],true)[2]
		rate		= R!(rt,[xc],xd,ti,Float64[],true)[1]
		@assert rate <= lambda_star "Wrong bound"
		reject = rand() < (1 - rate / lambda_star)
		S = -log(rand()) / lambda_star
		if ~reject
			verbose && println("----> Jump!, ratio = ", rate / lambda_star)
			@assert rate <= lambda_star "Wrong bound"
			push!(th,t)
			push!(xch,xc)
			push!(xdh,xdh[end] + 1)
			njumps += 1
			# dummy call to rand to emulate sampling pfsample
			dum = -log(rand())
		end
	end
	return th, xch, xdh
end

function F!(ẋ, xc, xd, parms, t)
	if mod(xd[1], 2)==0
		ẋ[1] = 10xc[1]
	else
		ẋ[1] = -3xc[1]^2
	end
end

R(x) = x

function R!(rate, xc, xd, parms, t, issum::Bool)
	# rate function
	if issum == false
		rate[1] = R(xc[1])
		rate[2] = parms[1]
		return 0., parms[1] + 50.
	else
		return R(xc[1]) + parms[1], parms[1] + 50.
	end
end

xc0 = [1.0]
	xd0 = [0, 0]

	nu = [1 0;0 -1]
	parms = [.0]
	ti = 0.322156
	tf = 100000.
	nj = 50

errors = Float64[]

Random.seed!(8)
	res_a_chv = AnalyticalSampleCHV(xc0,xd0,ti,nj)
Random.seed!(8)
	res_a_rej = AnalyticalSampleRejection(xc0,xd0,ti,nj)

problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
println("\n\nComparison of solvers")
	for ode in [(:cvode,"cvode"),(:lsoda,"lsoda"),(CVODE_BDF(),"CVODEBDF"),(CVODE_Adams(),"CVODEAdams"),(Tsit5(),"tsit5"),(Rodas4P(autodiff=false),"rodas4P-noAutoDiff"),(Rodas4P(),"rodas4P-AutoDiff"),(Rosenbrock23(),"RS23"),(AutoTsit5(Rosenbrock23(autodiff=true)),"AutoTsit5-RS23")]
	Random.seed!(8)
	res =  PDMP.solve(problem, CHV(ode[1]); n_jumps = nj)
	println("--> norm difference = ", norm(res.time - res_a_chv[1], Inf64), "  - solver = ",ode[2])
	@test norm(res.time - res_a_chv[1],Inf64) < 1e3
end

println("\n\nComparison of solvers - rejection")
	for ode in [(:cvode,"cvode"),(:lsoda,"lsoda"),(CVODE_BDF(),"CVODEBDF"),(CVODE_Adams(),"CVODEAdams"),(Tsit5(),"tsit5"),(Rodas4P(autodiff=false),"rodas4P-noAutoDiff"),(Rodas4P(),"rodas4P-AutoDiff"),(Rosenbrock23(),"RS23"),(AutoTsit5(Rosenbrock23(autodiff=true)),"AutoTsit5-RS23")]
	Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	res =  PDMP.solve(problem, Rejection(ode[1]); n_jumps = 4, verbose = false)
	println("--> norm difference = ", norm(res.time - res_a_rej[1][1:4], Inf64), "  - solver = ",ode[2])
	@test norm(res.time - res_a_rej[1][1:4], Inf64) < 0.0043
end

Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	res1 = PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj)
	Random.seed!(8)
	res2 = PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj)

println("Alternate between calls CHV - Rej")
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	res_chv = PDMP.solve(problem, CHV(Tsit5()); n_jumps = 50)
	Random.seed!(8)
	res_rej = PDMP.solve(problem, Rejection(Tsit5()); n_jumps = 4)
	println("--> norm diff (CHV) = ", norm(res_chv.time - res_a_chv[1], Inf64))
	println("--> norm diff (Rej) = ", norm(res_rej.time - res_a_rej[1][1:4], Inf64))
	@test norm(res_chv.time - res_a_chv[1], Inf64) < 1e-3
	@test norm(res_rej.time - res_a_rej[1][1:4], Inf64) < 0.0043

# here, we write the jump problem with a function
function Delta!(xc, xd, t, parms, ind_reaction::Int64)
	if ind_reaction == 1
		xd[1] += 1
	else
		xd[2] -= 1
	end
	nothing
end

println("\n\nComparison of solvers, with function Delta")
	for ode in [(:cvode,"cvode"),(:lsoda,"lsoda"),(CVODE_BDF(),"CVODEBDF"),(CVODE_Adams(),"CVODEAdams"),(Tsit5(),"tsit5"),(Rodas4P(autodiff=false),"rodas4P-noAutoDiff"),(Rodas4P(),"rodas4P-AutoDiff"),(Rosenbrock23(),"RS23"),(AutoTsit5(Rosenbrock23()),"AutoTsit5RS23")]
	Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!,  Delta!, 2, xc0, xd0, parms, (ti, tf))
	res =  PDMP.solve(problem, CHV(ode[1]); n_jumps = nj)
	println("--> norm difference = ", norm(res.time - res_a_chv[1],Inf64), "  - solver = ", ode[2])
	push!(errors, norm(res.time - res_a_chv[1], Inf64))
end

Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	res =  PDMP.solve(problem, CHV(:lsoda); n_jumps = nj)

# test for allocations, should not depend on the requested number of jumps
Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	alloc1 =  PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, false))
	alloc1 =  @allocated PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, false))
	Random.seed!(8)
	alloc2 =  @allocated PDMP.solve(problem, CHV(Tsit5()); n_jumps = 2nj, save_positions = (false, false))
	println("--> allocations = ", (alloc1, alloc2))

# test for many calls to solve, the trajectories should be the same
problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	res = PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, true))
	restime1 = res.time
	Random.seed!(8)
	res12 = PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, true))


text:
using CodeTracer
using Test

@testset "CodeTracer.jl" begin
    # Write your tests here.
end


text:
module DeprecatedTest

using Test
using MathOptInterface
using ConstraintProgrammingExtensions

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIT = MOI.DeprecatedTest
const CP = ConstraintProgrammingExtensions

# Test sets.
include("alldifferent.jl")
include("antidomain.jl")
include("binpacking.jl")
include("count.jl")
include("countdistinct.jl")
include("differentfrom.jl")
include("domain.jl")
include("element.jl")
include("equivalence.jl")
include("equivalencenot.jl")
include("ifthenelse.jl")
include("implication.jl")
include("indicator.jl")
include("inverse.jl")
include("lexicographically.jl")
include("minimumdistance.jl")
include("reification.jl")
include("strictly.jl")
include("strictly_lexicographically.jl")
include("truefalse.jl")

end


text:
module LispSyntax

include("parser.jl")
export sx, desx, codegen, @lisp_str, assign_reader_dispatch, include_lisp

# Internal types
mutable struct s_expr
  vector
end

sx(x...) = s_expr([x...])
==(a :: s_expr, b :: s_expr) = a.vector == b.vector

function desx(s)
  if typeof(s) == s_expr
    return map(desx, s.vector)
  elseif isa(s, Dict)
    return Dict(desx(x[1])=>desx(x[2]) for x in s)
  elseif isa(s, Set)
    return Set(desx(v) for v in s)
  else
    return s
  end
end

function lispify(s)
  if isa(s, s_expr)
    return "(" * join(map(lispify, s.vector), " ") * ")"
  else
    return "$s"
  end
end

# convert the input tuple to an array
construct_sexpr(items...) = Any[items...]

function assign_reader_dispatch(sym, fn)
  reader_table[sym] = fn
end

function quasiquote(s)
  if isa(s, Array) && length(s) == 2 && s[1] == :splice
    codegen(s[2])
  elseif isa(s, Array) && length(s) == 2 && s[1] == :splice_seq
    Expr(:..., codegen(s[2]))
  elseif isa(s, Array)
    Expr(:call, construct_sexpr, map(quasiquote, s)...)
  elseif isa(s, Symbol)
    Expr(:quote, s)
  else
    s
  end
end

function quote_it(s)
  if isa(s, Array)
    Expr(:call, construct_sexpr, map(s -> quote_it(s), s)...)
  elseif isa(s, Symbol)
   QuoteNode(s)
  else
    s
  end
end

function codegen(s)
  if isa(s, Symbol)
      s
  elseif isa(s, Dict)
    coded_s = [:($(codegen(x[1])) => $(codegen(x[2]))) for x in s]
    Expr(:call, Dict, coded_s...)
  elseif isa(s, Set)
    coded_s = [codegen(x) for x in s]
    Expr(:call, Set, Expr(:vect, coded_s...))
  elseif s isa Expr && s.head == :escape
    # Special case to allow use of `esc` in lisp syntax macros
    esc(codegen(s.args[1]))
  elseif !isa(s, Array) # constant
    s
  elseif length(s) == 0 # empty array
    s
  elseif s[1] == :if
    if length(s) == 3
      :($(codegen(s[2])) && $(codegen(s[3])))
    elseif length(s) == 4
      :($(codegen(s[2])) ? $(codegen(s[3])) : $(codegen(s[4])))
    else
      throw("illegal if statement $s")
    end
  elseif s[1] == :def
    length(s) == 3 || error("Malformed def: Length of list must be == 3")
    :(global $(s[2]) = $(codegen(s[3])))
  elseif s[1] == :let
    bindings = [ :($(s[2][i]) = $(codegen(s[2][i+1]))) for i = 1:2:length(s[2]) ]
    coded_s  = map(codegen, s[3:end])
    Expr(:let, Expr(:block, bindings...), Expr(:block, coded_s...))
  elseif s[1] == :while
    coded_s = map(codegen, s[2:end])
    Expr(:while, coded_s[1], Expr(:block, coded_s[2:end]...))
  elseif s[1] == :for
    bindings = [ :($(s[2][i]) = $(codegen(s[2][i+1]))) for i = 1:2:length(s[2]) ]
    coded_s  = map(codegen, s[3:end])
    Expr(:for, Expr(:block, bindings...), Expr(:block, coded_s...))
  elseif s[1] == :do
    Expr(:block, map(codegen, s[2:end])...)
  elseif s[1] == :global
    Expr(:global, s[2:end]...)
  elseif s[1] == :quote
    quote_it(s[2])
  elseif s[1] == :import
    Expr(:using, [Expr(:., x) for x in s[2:end]]...)
  elseif s[1] == :splice
    throw("missplaced ~ (splice)")
  elseif s[1] == :splice_seq
    throw("missplaced ~@ (splice_seq)")
  elseif s[1] == :quasi
    quasiquote(s[2])
  elseif s[1] == :lambda || s[1] == :fn
    length(s) >= 3 || error("Malformed lambda/fn: list length must be >= 3")
    coded_s = map(codegen, s[3:end])
    Expr(:function, Expr(:tuple, s[2]...), Expr(:block, coded_s...))
  elseif s[1] == :defn
    # NB: This lowering of `defn` makes a julia function which may be a closure
    # if used in local scope. This is a semantic mismatch with clojure where
    # `defn` binds a lambda to a mutable global name. We could do it the
    # closure way, but it would generate much worse code.
    coded_s = map(codegen, s[4:end])
    Expr(:function, Expr(:call, s[2], s[3]...), Expr(:block, coded_s...))
  elseif s[1] == :defmacro
    # NB: Clojure macros are unhygenic by default, the opposite of Julia.
    # Choose Julia semantics and allow the use of `esc` intermixed with lists
    # (see `Expr(:escape)` handling above).
    Expr(:macro, Expr(:call, s[2], s[3]...),
         begin
             sexpr = Expr(:block, map(codegen, s[4:end])...)
             Expr(:block, Expr(:call, codegen, sexpr))
         end)
  elseif s[1] == :defmethod
    # TODO
  else
    coded_s = map(codegen, s)
    if (typeof(coded_s[1]) == Symbol && occursin(r"^@.*$", string(coded_s[1]))) ||
       (typeof(coded_s[1]) == Expr && occursin(r"^@.*$", string(coded_s[1].args[1])))
      Expr(:macrocall, coded_s[1], nothing, coded_s[2:end]...)
    else
      Expr(:call, coded_s[1], coded_s[2:end]...)
    end
  end
end

"This is an internal helper function, do not call outside of package"
function lisp_eval_helper(str :: AbstractString)
  s = desx(LispSyntax.read(str))
  return codegen(s)
end

macro lisp_str(str)
  return esc(lisp_eval_helper(str))
end

"""
    include_lisp(mod, source)

Parse expressions from `source` (which may be an `IO` or file name) and
evaluate them sequentially as top level code in module `mod`.
"""
function include_lisp(mod::Module, filename::AbstractString)
    open(filename) do io
        include_lisp(mod, io)
    end
end

function include_lisp(mod::Module, io::IO)
    content = Base.read(io, String)
    res = nothing
    for sxpr in parse_one(content, top_level)
        ex = codegen(desx(sxpr))
        res = Base.eval(mod, ex)
    end
    res
end

include("repl.jl")

end # module


text:
struct DynamicsData{T}
	m
	ip::InteriorPoint
    z_subset_init::Vector{T}
	θ_params
	h::T
	diff_idx::Int
	x_view::SubArray{T,1,Vector{T},Tuple{Vector{Int}},false}
	δx_view::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}
	δu_view::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}
	δx_cache::Vector{Matrix{T}}
	δu_cache::Vector{Matrix{T}}
	z_cache::Vector{Vector{T}}
end

function dynamics_data(m, h, T,
        r!, rz!, rθ!, rz, rθ;
        idx_ineq = collect(1:0),
        idx_soc = Vector{Int}[],
        z_subset_init = ones(size(rz)[1] - m.dim.q),
		θ_params = [],
		diff_idx = -1,
		opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = false))

    nz = size(rz)[1]
    nθ = size(rθ)[2]
	z = zeros(nz)
	θ = zeros(nθ)

	ip = interior_point(z, θ,
		idx_ineq = idx_ineq,
        idx_soc = idx_soc,
		r! = r!,
		rz! = rz!,
		rθ! = rθ!,
		rz = rz,
		rθ = rθ,
		opts = opts)

	ip.opts.diff_sol = false

	x_view = view(ip.z, collect(1:m.dim.q))
	δx_view = view(ip.δz, collect(1:m.dim.q), collect(1:2 * m.dim.q))
	δu_view = view(ip.δz, collect(1:m.dim.q), collect(2 * m.dim.q .+ (1:m.dim.u)))
	δx_cache = [zeros(m.dim.q, 2 * m.dim.q) for t = 1:T-1]
	δu_cache = [zeros(m.dim.q, m.dim.u) for t = 1:T-1]

	z_cache = [zero(ip.z) for t = 1:T-1]

	DynamicsData(m, ip, z_subset_init, θ_params, h,
		diff_idx, x_view, δx_view, δu_view, δx_cache, δu_cache,
		z_cache)
end

function f!(d::DynamicsData, q0, q1, u1, mode = :dynamics)
	ip = d.ip
	h = d.h

	ip.z .= copy([q1; d.z_subset_init])
	ip.θ .= copy([q0; q1; u1; h; d.θ_params])

	status = interior_point_solve!(ip)

	!status && (@warn "dynamics failure (res norm: $(norm(ip.r, Inf))) \n
		               z = $(ip.z), \n
					   θ = $(ip.θ)")
end

function f(d::DynamicsData, q0, q1, u1, t)
	f!(d, q0, q1, u1, :dynamics)
	differentiate_solution!(d.ip, z = (d.diff_idx == -1 ? d.ip.z : d.ip.z_cache[d.diff_idx]))
	d.δx_cache[t] .= d.δx_view
	d.δu_cache[t] .= d.δu_view
	d.z_cache[t] .= d.ip.z
	return d.x_view
end

function fx1(d::DynamicsData, q0, q1, u1, t)
	return d.δx_cache[t]
end

function fu1(d::DynamicsData, q0, q1, u1, t)
	return d.δu_cache[t]
end

struct ImplicitDynamics{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int
	dynamics::DynamicsData
	cache::Dict
end

function ImplicitDynamics{I,T}(n, m, d, dynamics; cache=Dict()) where {I,T}
	ImplicitDynamics{I,T}(n, m, d, dynamics, cache)
end

function fd(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = view(x, 1:nq)
	q1 = view(x, nq .+ (1:nq))

	q2 = f(model.dynamics, q0, q1, u, t)

	return [q1; q2]
end

function fdx(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = view(x, 1:nq)
	q1 = view(x, nq .+ (1:nq))
	dq2dx1 = fx1(model.dynamics, q0, q1, u, t)

	return [zeros(nq, nq) I; dq2dx1]
end

function fdu(model::ImplicitDynamics{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.m.dim.q
	q0 = view(x, 1:nq)
	q1 = view(x, nq .+ (1:nq))
	dq2du1 = fu1(model.dynamics, q0, q1, u, t)
	return [zeros(nq, model.m); dq2du1]
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
function _pwf2pm_load!(pm_data::Dict, pwf_data::Dict, i::Int)
    bus = pwf_data["DBAR"]["$i"]

    sub_data = Dict{String,Any}()

    sub_data["load_bus"] = bus["NUMBER"]
    sub_data["pd"] = pop!(bus, "ACTIVE CHARGE")
    sub_data["qd"] = pop!(bus, "REACTIVE CHARGE")
    sub_data["status"] = 1

    sub_data["source_id"] = ["load", sub_data["load_bus"], "1 "]
    sub_data["index"] = length(pm_data["load"]) + 1

    idx = string(sub_data["index"])
    pm_data["load"][idx] = sub_data
end

function _pwf2pm_load!(pm_data::Dict, pwf_data::Dict, bus::Dict)
    sub_data = Dict{String,Any}()

    sub_data["load_bus"] = bus["NUMBER"]
    sub_data["pd"] = pop!(bus, "ACTIVE CHARGE")
    sub_data["qd"] = pop!(bus, "REACTIVE CHARGE")
    sub_data["status"] = 1

    sub_data["source_id"] = ["load", sub_data["load_bus"], "1 "]
    sub_data["index"] = length(pm_data["load"]) + 1

    idx = string(sub_data["index"])
    pm_data["load"][idx] = sub_data
end

function _pwf2pm_load!(pm_data::Dict, pwf_data::Dict)

    pm_data["load"] = Dict{String, Any}()
    if haskey(pwf_data, "DBAR")
        for (i,bus) in pwf_data["DBAR"]
            if bus["ACTIVE CHARGE"] > 0.0 || bus["REACTIVE CHARGE"] > 0.0 || bus["TYPE"] == bus_type_raw_to_pwf[bus_type_str_to_num["PQ"]]
                _pwf2pm_load!(pm_data, pwf_data, bus["NUMBER"])
            end
        end
    end
end


text:
module Options_Test

using Compat
using Compat.Test
using EAGOBranchBound

x = BnBSolver()

set_Branch_Scheme!(x,"best")
@test x.Node_Select == EAGOBranchBound.NS_best
@test x.Branch_Sto == EAGOBranchBound.BM_depth_best!

set_Branch_Scheme!(x,"breadth")
@test x.Node_Select == EAGOBranchBound.NS_depth_breadth
@test x.Branch_Sto == EAGOBranchBound.BM_breadth!

set_Branch_Scheme!(x,"depth")
@test x.Node_Select == EAGOBranchBound.NS_depth_breadth
@test x.Branch_Sto == EAGOBranchBound.BM_depth_best!

set_Bisect_Func!(x,"relative midpoint",-1)
@test x.Bisect_Func == EAGOBranchBound.Bisect_Rel
set_Bisect_Func!(x,"absolute midpoint",-1)
@test x.Bisect_Func == EAGOBranchBound.Bisect_Abs

set_Verbosity!(x,"Normal")
@test x.Verbosity == "Normal"

@test_throws ErrorException set_Branch_Scheme!(x,"random")
@test_throws ErrorException  set_Verbosity!(x,"random")
@test_throws ErrorException  set_Bisect_Func!(x,"random",-1)
@test_throws ErrorException  set_Bisect_Func!(x,"random",1)

@test_nowarn set_Bisect_Func!(x,"relative midpoint",3)
@test_nowarn set_Bisect_Func!(x,"absolute midpoint",3)

end


text:
"""
    @polyeval(z, p, N)

Evaluate the polynomial ``\\sum_k c[k] z^{k-1}`` for the coefficients `c[1]`, `c[2]`, ..., `c[N]`;
that is, the coefficients are given in ascending order by power of `z`.  This macro expands
to efficient inline code that uses either Horner's method.

```jldoctest

julia> p = randn(5)
5-element Array{Float64,1}:
  0.106455
  0.0351716
 -0.204764
 -0.87642
  0.260407

julia> @polyeval(0.5, p, 5)
-0.020427265391245605

```
"""
macro polyeval(x, p, N)
    if N==1
        return :($(esc(p))[1])
    elseif N==2
        return :($(esc(p))[1] + $(esc(x))*$(esc(p))[2])
    else
        nm1 = N-1
        ex = :($(esc(p))[$nm1] + $(esc(x))*$(esc(p))[$N])
        for i in N-2:-1:1
            ex = :($(esc(p))[$i] + $(esc(x))*$ex)
        end
    end
    return ex
end

"""
    ConvergenceError([msg, val, niter, err])

Exception that is to be thrown when an iterative procedure fails
to converge.

"""
mutable struct ConvergenceError <: Exception
    msg::String
    val::Float64
    niter::Int
    err::Float64
    ConvergenceError() = new("Calculations failed to converge", NaN, -1, NaN)
    ConvergenceError(msg) = new(msg, NaN, -1, NaN)
    ConvergenceError(msg, val) = new(msg, val, -1, NaN)
    ConvergenceError(msg, val, niter) = new(msg, val, niter, NaN)
    ConvergenceError(msg, val, niter, err) = new(msg, val, niter, err)
end

"""
    calcz(vm, b0, c0, [EPS, [MAXITER, [relax]]])

Calculates the compressibility factor using a virial equation:

`z = 1 + b₀/z + c₀/vₘ²`

The algorithm initially uses as an initial guess only the b₀ virial
coefficient. Then it uses a Newton-Raphson algorithm to compute the
compressibility factor

 * `b0` b₀ parameter of the virial equation
 * `c0` c₀ parameter of the virial equation
 * `EPS` Convergence criterium
 * `MAXITER` Maximim number of iterations
 * `relax` Sub-relaxation parameter if convergence is difficult.
"""
function calcz(b0, c0, EPS=1e-8, MAXITER=200, relax=1.0)

    z = (1.0 + sqrt(1.0 + 4*b0)) / 2.0
    dz = 0.0

    for i = 1:MAXITER
        f = -c0 + z*(-b0 + z*(-1.0 + z))
        df = -b0 + z*(-2.0 + 3.0*z)
        dz = -f / df
        if abs(dz) < EPS
            return z + dz
        end
        z = z + relax*dz
    end
    throw(ConvergenceError("Compressibility factor failed to converge properly!",
                           z, MAXITER, abs(dz)))

end


text:
# Nibble operations
# =================
#
# Copyright Ben J. Ward and Luis Yanes.
# License is MIT.

"""
    enumerate_nibbles(abxor::UInt64)

Count the number of set bits in each nibble (aligned 4 bit segments) of an unsigned 64 bit integer.

E.g. An input of:

0100 0010 0001 0110 1100 1110 1101 1111

Would result in:

0001 0001 0001 0010 0010 0011 0011 0100

This is used to identify different occurances of certain bit patterns.
"""
@inline function enumerate_nibbles(x::UInt64)
    x = x - ((x >> 1) & 0x5555555555555555)
    return (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
end

"""
    count_nonzero_nibbles(x::UInt64)

Count the number of nibbles (aligned 4 bit segments) in an unsigned integer which have at least one
bit set.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 15.
"""
@inline function count_nonzero_nibbles(x::UInt64)
    out = UInt64(0)
    out |= x & 0x1111111111111111
    out |= (x & 0x2222222222222222) >> 1
    out |= (x & 0x4444444444444444) >> 2
    out |= (x & 0x8888888888888888) >> 3
    return count_ones(out)
end

"""
    count_zero_nibbles(x::UInt64)

Counts the number of nibbles (aligned 4 bit segments) in a UInt64 `x` that have all their bits unset i.e.
nibbles of 0000.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 1.
"""
@inline function count_zero_nibbles(x::UInt64)
    return 16 - count_nonzero_nibbles(x)
end

"""
    count_one_nibbles(x::UInt64)

Counts the number of nibbles (aligned 4 bit segments) in a UInt64 `x` that have all their bits set i.e.
all nibbles of 1111.

E.g. An input of:

0x0F11F111F11111F1

Would give the answer: 4.
"""
@inline function count_one_nibbles(x::UInt64)
    out = x & 0x1111111111111111
    out &= (x & 0x2222222222222222) >> 1
    out &= (x & 0x4444444444444444) >> 2
    out &= (x & 0x8888888888888888) >> 3
    return count_ones(out)
end

"""
    nibble_mask(x::UInt64, value::UInt64)

Create a mask for the nibbles (aligned 4 bit segments) in a 64 bit integer `x`
that match a given value dictated by the pattern in `value`.
"""
@inline function nibble_mask(value::UInt64, x::UInt64)
    # XOR with the desired values. So matching nibbles will be 0000.
    x $= value
    # Horizontally OR the nibbles.
    x |= (x >> 1)
    x |= (x >> 2)
    # AND removes junk, we then widen x by multiplication and return
    # the inverse.
    x &= 0x1111111111111111
    x *= 15
    return ~x
end


text:
include("BPX1D.jl")
include("BPX2D.jl")
include("containers.jl")
include("Losses.jl")
using SparseArrays

function transform_T(T_prec::T_preconditioner)
    prec = preconditioner(T_prec.A, T_prec.Ind, T_prec.BCs, T_prec.L, T_prec.name, T_prec.symmetry)
end

function assemble_matrix(parameters::AbstractArray{<:Real, 1}, prec::preconditioner)
    B = nothing
    if prec.name == "full BPX"
        B = I + sum([parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]*BPX_layer_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs) for l=1:(prec.L-1)])
    elseif prec.name == "symmetric BPX"
        B = I + sum([parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]][1]*BPX_layer_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs) for l=1:(prec.L-1)])
    elseif prec.name == "reduced BPX"
        B = I + sum([parameters[2^prec.L+l]*BPX_layer_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs) for l=1:(prec.L-1)])
    elseif prec.name == "reduced+symmetric BPX"
        B = I + sum([parameters[2^(prec.L-1)+l]*BPX_layer_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs) for l=1:(prec.L-1)])
    elseif prec.name == "full Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        P = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[1][1][1]:prec.Ind[1][1][2]], parameters[prec.Ind[1][2][1]:prec.Ind[1][2][2]], prec.BCs)
        A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
        B += parameters[prec.Ind[1][3][1]:prec.Ind[1][3][2]][1]*P*A_P*transpose(P)
        for l=2:(prec.L-1)
            P = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs)
            A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            B += parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]*P*A_P*transpose(P)
        end
    elseif prec.name == "symmetric Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        P = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[1][1][1]:prec.Ind[1][1][2]], vcat(reverse(parameters[prec.Ind[1][1][1]:prec.Ind[1][1][2]])[2:end], 0), prec.BCs)
        A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
        B += parameters[prec.Ind[1][2][1]:prec.Ind[1][2][2]][1]*P*A_P*transpose(P)
        for l=2:(prec.L-1)
            P = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs)
            A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            B += parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]][1]*P*A_P*transpose(P)
        end
    elseif prec.name == "reduced Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        P = BPX_layer_half_matrix(prec.L, parameters[1:(2^(prec.L-1))], parameters[(2^(prec.L-1)+1):(2^prec.L)], prec.BCs)
        A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
        B += parameters[2^prec.L+1]*P*A_P*transpose(P)
        for l=2:(prec.L-1)
            P = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs)
            A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            B += parameters[2^prec.L+l]*P*A_P*transpose(P)
        end
    elseif prec.name == "reduced+symmetric Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        P = BPX_layer_half_matrix(prec.L, parameters[1:(2^(prec.L-1))], vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0), prec.BCs)
        A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
        B += parameters[2^(prec.L-1)+1]*P*A_P*transpose(P)
        for l=2:(prec.L-1)
            P = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs)
            A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            B += parameters[2^(prec.L-1)+1]*P*A_P*transpose(P)
        end
    elseif prec.name == "2D full BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]], parameters[prec.Ind[l][4][1]:prec.Ind[l][4][2]], prec.BCs[3:4])
            P = kron(P_y, P_x)
            B += (P*transpose(P))*parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]][1]
        end
    elseif prec.name == "2D symmetric BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], vcat(reverse(parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]])[2:end], 0), prec.BCs[3:4])
            P = kron(P_y, P_x)
            B += (P*transpose(P))*parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]
        end
    elseif prec.name == "2D reduced BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L)+1):(3*2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(3*2^(prec.L-1)+1):(2^(prec.L+1))])[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            B += (P*transpose(P))*parameters[2^(prec.L+1)+l]
        end
    elseif prec.name == "2D reduced+symmetric BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            B += (P*transpose(P))*parameters[2^(prec.L)+l]
        end
    elseif prec.name == "2D full Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]], parameters[prec.Ind[l][4][1]:prec.Ind[l][4][2]], prec.BCs[3:4])
            P = kron(P_y, P_x)
            if l == 1
                A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
            else
                A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            end
            B += P*A_P*transpose(P)*parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]][1]
        end
    elseif prec.name == "2D symmetric Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], vcat(reverse(parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]])[2:end], 0), prec.BCs[3:4])
            P = kron(P_y, P_x)
            if l == 1
                A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
            else
                A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            end
            B += P*A_P*transpose(P)*parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]
        end
    elseif prec.name == "2D reduced Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L)+1):(3*2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(3*2^(prec.L-1)+1):(2^(prec.L+1))])[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            if l == 1
                A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
            else
                A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            end
            B += P*A_P*transpose(P)*parameters[2^(prec.L+1)+l]
        end
    elseif prec.name == "2D reduced+symmetric Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            if l == 1
                A_P = sparse(inv(Array(transpose(P)*prec.A*P) + 1e-6*I))
            else
                A_P = spdiagm(0 => diag(transpose(P)*prec.A*P).^-1)
            end
            B += P*A_P*transpose(P)*parameters[2^(prec.L)+l]
        end
    elseif prec.name == "2D tridiagonal BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            N_x, N_y = size(P_x)[2], size(P_y)[2]
            T_x = spdiagm(0 => parameters[2^(prec.L)+1+6*(l-1)]*ones(N_x), 1 => parameters[2^(prec.L)+2+6*(l-1)]*ones(N_x-1), -1 => parameters[2^(prec.L)+3+6*(l-1)]*ones(N_x-1))
            T_y = spdiagm(0 => parameters[2^(prec.L)+4+6*(l-1)]*ones(N_y), 1 => parameters[2^(prec.L)+5+6*(l-1)]*ones(N_y-1), -1 => parameters[2^(prec.L)+6+6*(l-1)]*ones(N_y-1))
            P = kron(P_y, P_x)*kron(T_y, T_x)
            # P = kron(P_y, P_x)*(kron(spdiagm(0 => ones(N_y)), T_x) + kron(T_y, spdiagm(0 => ones(N_x))))
            B += parameters[6*(prec.L-1)+2^(prec.L)+l]*P*transpose(P)
        end
    elseif prec.name == "2D double Schwarz"
        B = spdiagm(0 => diag(prec.A).^-1)
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P_x_ = BPX_layer_half_matrix(prec.L, (parameters[(2^prec.L+1):(2^prec.L+2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^prec.L+1):(2^prec.L+2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y_ = BPX_layer_half_matrix(prec.L, (parameters[(2^prec.L+2^(prec.L-1)+1):(2^prec.L+2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^prec.L+2^(prec.L-1)+1):(2^prec.L+2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            P_ = kron(P_y_, P_x_)
            if l == 1
                A_P = sparse(inv(Array(transpose(P_)*prec.A*P_) + 1e-6*I))
            else
                A_P = spdiagm(0 => diag(transpose(P_)*prec.A*P_).^-1)
            end
            B += P*A_P*transpose(P)*parameters[2^(prec.L+1)+l]
        end
    elseif prec.name == "2D reduced+symmetric+rescaled BPX"
        B = spdiagm(0 => diag(prec.A).^(-1/2))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P = kron(P_y, P_x)
            ψ, χ = BPX_weights(prec.L, l)
            P_x_ = BPX_layer_half_matrix(prec.L, ψ, χ, prec.BCs[1:2])
            P_y_ = BPX_layer_half_matrix(prec.L, ψ, χ, prec.BCs[3:4])
            P_ = kron(P_y_, P_x_)
            A_P = spdiagm(0 => diag(transpose(P_)*prec.A*P_).^(-1/2))
            B += P*A_P*transpose(P)*parameters[2^(prec.L)+l]
        end
    elseif prec.name == "2D reduced+symmetric+semi BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        shift = 1
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            l_ = min(l + shift, (prec.L - 1))
            if l_ == (prec.L - 1)
                P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l_-1):2^(l_-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l_-1):2^(l_-1):end], prec.BCs[3:4])
                P_y = spdiagm(0 => ones(size(P_y)[1]))
            else
                P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l_-1):2^(l_-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l_-1):2^(l_-1):end], prec.BCs[3:4])
            end
            P = kron(P_y, P_x)
            B += (P*transpose(P))*parameters[2^(prec.L)+l]
        end
    end
    return B
end

function assemble_matrix(parameters::AbstractArray{<:Real, 1}, T_prec::T_preconditioner)
    assemble_matrix(parameters, transform_T(T_prec))
end

function apply_T_preconditioner(X::AbstractArray{<:Real, 2}, parameters::AbstractArray{<:Real, 1}, prec::T_preconditioner)
    Y = X .+ 0.0
    if prec.name == "2D full BPX"
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]], parameters[prec.Ind[l][4][1]:prec.Ind[l][4][2]], prec.BCs[3:4])
            Y += parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]][1]*(((P_y*(transpose(P_y)*X))*P_x)*transpose(P_x))
        end
    elseif prec.name == "2D symmetric BPX"
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], vcat(reverse(parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]])[2:end], 0), prec.BCs[3:4])
            Y += parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]*((P_y*((transpose(P_y)*X)*P_x))*transpose(P_x))
        end
    elseif prec.name == "2D reduced BPX"
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L)+1):(3*2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(3*2^(prec.L-1)+1):(2^(prec.L+1))])[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            Y += parameters[2^(prec.L+1)+l]*((P_y*((transpose(P_y)*X)*P_x))*transpose(P_x))
        end
    elseif prec.name == "2D reduced+symmetric BPX"
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            Y += parameters[2^(prec.L)+l]*((P_y*((transpose(P_y)*X)*P_x))*transpose(P_x))
        end
    elseif prec.name == "2D full Schwarz"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]], parameters[prec.Ind[l][4][1]:prec.Ind[l][4][2]], prec.BCs[3:4])
            if l == 1
                A_inv = inv(Array(sum([kron((transpose(P_y)*c[1])*P_y, (transpose(P_x)*c[2])*P_x) for c = prec.C]) + 1e-6*I))
                y = ((transpose(P_y)*X)*P_x)
                Y += ((P_y*transpose(reshape(A_inv*reshape(transpose(y), (size(A_inv)[1], 1)), (size(y)[2], size(y)[1]))))*transpose(P_x))*parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]][1]
            else
                D_inv = 1 ./ sum([(diag(transpose(P_y)*c[1]*P_y))*transpose(diag(transpose(P_x)*c[2]*P_x)) for c = prec.C])
                Y += ((P_y*(((transpose(P_y)*X)*P_x) .* D_inv))*transpose(P_x))*parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]][1]
            end
        end
    elseif prec.name == "2D symmetric Schwarz"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]], vcat(reverse(parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]])[2:end], 0), prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]], vcat(reverse(parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]])[2:end], 0), prec.BCs[3:4])
            if l == 1
                A_inv = inv(Array(sum([kron((transpose(P_y)*c[1])*P_y, (transpose(P_x)*c[2])*P_x) for c = prec.C]) + 1e-6*I))
                y = ((transpose(P_y)*X)*P_x)
                Y += ((P_y*transpose(reshape(A_inv*reshape(transpose(y), (size(A_inv)[1], 1)), (size(y)[2], size(y)[1]))))*transpose(P_x))*parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]
            else
                D_inv = 1 ./ sum([(diag(transpose(P_y)*c[1]*P_y))*transpose(diag(transpose(P_x)*c[2]*P_x)) for c = prec.C])
                Y += ((P_y*(((transpose(P_y)*X)*P_x) .* D_inv))*transpose(P_x))*parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]][1]
            end
        end
    elseif prec.name == "2D reduced Schwarz"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L)+1):(3*2^(prec.L-1))])[2^(l-1):2^(l-1):end], (parameters[(3*2^(prec.L-1)+1):(2^(prec.L+1))])[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            if l == 1
                A_inv = inv(Array(sum([kron((transpose(P_y)*c[1])*P_y, (transpose(P_x)*c[2])*P_x) for c = prec.C]) + 1e-6*I))
                y = ((transpose(P_y)*X)*P_x)
                Y += ((P_y*transpose(reshape(A_inv*reshape(transpose(y), (size(A_inv)[1], 1)), (size(y)[2], size(y)[1]))))*transpose(P_x))*parameters[2^(prec.L+1)+l]
            else
                D_inv = 1 ./ sum([(diag(transpose(P_y)*c[1]*P_y))*transpose(diag(transpose(P_x)*c[2]*P_x)) for c = prec.C])
                Y += ((P_y*(((transpose(P_y)*X)*P_x) .* D_inv))*transpose(P_x))*parameters[2^(prec.L+1)+l]
            end
        end
    elseif prec.name == "2D reduced+symmetric Schwarz"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            if l == 1
                A_inv = inv(Array(sum([kron((transpose(P_y)*c[1])*P_y, (transpose(P_x)*c[2])*P_x) for c = prec.C]) + 1e-6*I))
                y = ((transpose(P_y)*X)*P_x)
                Y += ((P_y*transpose(reshape(A_inv*reshape(transpose(y), (size(A_inv)[1], 1)), (size(y)[2], size(y)[1]))))*transpose(P_x))*parameters[2^(prec.L)+l]
            else
                D_inv = 1 ./ sum([(diag(transpose(P_y)*c[1]*P_y))*transpose(diag(transpose(P_x)*c[2]*P_x)) for c = prec.C])
                Y += ((P_y*(((transpose(P_y)*X)*P_x) .* D_inv))*transpose(P_x))*parameters[2^(prec.L)+l]
            end
        end
    elseif prec.name == "2D double Schwarz"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            P_x_ = BPX_layer_half_matrix(prec.L, (parameters[(2^prec.L+1):(2^prec.L+2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^prec.L+1):(2^prec.L+2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y_ = BPX_layer_half_matrix(prec.L, (parameters[(2^prec.L+2^(prec.L-1)+1):(2^prec.L+2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^prec.L+2^(prec.L-1)+1):(2^prec.L+2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            if l == 1
                A_inv = inv(Array(sum([kron((transpose(P_y_)*c[1])*P_y_, (transpose(P_x_)*c[2])*P_x_) for c = prec.C]) + 1e-6*I))
                y = ((transpose(P_y)*X)*P_x)
                Y += ((P_y*transpose(reshape(A_inv*reshape(transpose(y), (size(A_inv)[1], 1)), (size(y)[2], size(y)[1]))))*transpose(P_x))*parameters[2^(prec.L)+l]
            else
                D_inv = 1 ./ sum([(diag(transpose(P_y_)*c[1]*P_y_))*transpose(diag(transpose(P_x_)*c[2]*P_x_)) for c = prec.C])
                Y += ((P_y*(((transpose(P_y)*X)*P_x) .* D_inv))*transpose(P_x))*parameters[2^prec.L+2^(prec.L)+l]
            end
        end
    elseif prec.name == "2D tridiagonal BPX"
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            N_x, N_y = size(P_x)[2], size(P_y)[2]
            T_x = spdiagm(0 => parameters[2^(prec.L)+1+6*(l-1)]*ones(N_x), 1 => parameters[2^(prec.L)+2+6*(l-1)]*ones(N_x-1), -1 => parameters[2^(prec.L)+3+6*(l-1)]*ones(N_x-1))
            T_y = spdiagm(0 => parameters[2^(prec.L)+4+6*(l-1)]*ones(N_y), 1 => parameters[2^(prec.L)+5+6*(l-1)]*ones(N_y-1), -1 => parameters[2^(prec.L)+6+6*(l-1)]*ones(N_y-1))
            P = kron(P_y, P_x)*(kron(spdiagm(0 => ones(N_y)), T_x) + kron(T_y, spdiagm(0 => ones(N_x))))
            Y += parameters[6*(prec.L-1)+2^(prec.L)+l]*(P_y*(T_y*(transpose(T_y)*(((transpose(P_y)*(X*P_x))*T_x)*transpose(T_x)))))*transpose(P_x)
            # z = transpose(P_y)*(X*P_x)
            # Y += parameters[6*(prec.L-1)+2^(prec.L)+l]*((P_y*((z*T_x)*transpose(T_x) + T_y*(transpose(T_y)*z) + (transpose(T_y)*z)*transpose(T_x) + (T_y*z)*T_x))*transpose(P_x))
        end
    elseif prec.name == "2D reduced+symmetric+rescaled BPX"
        Y ./= transpose(reshape(diag(prec.A), (size(X)[2], size(X)[1])))
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[3:4])
            ψ, χ = BPX_weights(prec.L, l)
            P_x_ = BPX_layer_half_matrix(prec.L, ψ, χ, prec.BCs[1:2])
            P_y_ = BPX_layer_half_matrix(prec.L, ψ, χ, prec.BCs[3:4])
            D_inv = 1 ./ sum([(diag(transpose(P_y_)*c[1]*P_y_))*transpose(diag(transpose(P_x_)*c[2]*P_x_)) for c = prec.C])
            Y += ((P_y*(((transpose(P_y)*X)*P_x) .* sqrt.(D_inv)))*transpose(P_x))*parameters[2^(prec.L)+l]
        end
    elseif prec.name == "2D reduced+symmetric+semi BPX"
        B = spdiagm(0 => ones(size(prec.A)[1]))
        shift = 1
        for l=1:(prec.L-1)
            P_x = BPX_layer_half_matrix(prec.L, (parameters[1:(2^(prec.L-1))])[2^(l-1):2^(l-1):end], (vcat(reverse(parameters[1:(2^(prec.L-1))])[2:end], 0))[2^(l-1):2^(l-1):end], prec.BCs[1:2])
            l_ = min(l + shift, (prec.L - 1))
            if l_ == (prec.L - 1)
                P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l_-1):2^(l_-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l_-1):2^(l_-1):end], prec.BCs[3:4])
                P_y = spdiagm(0 => ones(size(P_y)[1]))
            else
                P_y = BPX_layer_half_matrix(prec.L, (parameters[(2^(prec.L-1)+1):(2^prec.L)])[2^(l_-1):2^(l_-1):end], (vcat(reverse(parameters[(2^(prec.L-1)+1):(2^prec.L)])[2:end], 0))[2^(l_-1):2^(l_-1):end], prec.BCs[3:4])
            end
            Y += parameters[2^(prec.L)+l]*((P_y*((transpose(P_y)*X)*P_x))*transpose(P_x))
        end
    end
    return Y
end

function indices(BCs::Array{String, 1}, L::Integer, name::String)
    ind = []
    c = 0
    if name[1:4] == "full"
        for l=1:(L-1)
            vec_size = 2^(L-l)
            push!(ind, [[c+1, c+vec_size], [c+vec_size+1, c+2*vec_size], [c+2*vec_size+1, c+2*vec_size+1]])
            c += 2*vec_size+1
        end
    elseif name[1:3] == "sym"
        for l=1:(L-1)
            vec_size = 2^(L-l)
            push!(ind, [[c+1, c+vec_size], [c+vec_size+1, c+vec_size+1]])
            c += vec_size+1
        end
    elseif name[1:7] == "2D full"
        for l=1:(L-1)
            vec_size = 2^(L-l)
            push!(ind, [[c+1, c+vec_size], [c+vec_size+1, c+2*vec_size], [c+2*vec_size+1, c+3*vec_size], [c+3*vec_size+1, c+4*vec_size], [c+4*vec_size+1, c+4*vec_size+1]])
            c += 4*vec_size+1
        end
    elseif name[1:6] == "2D sym"
        for l=1:(L-1)
            vec_size = 2^(L-l)
            push!(ind, [[c+1, c+vec_size], [c+vec_size+1, c+2*vec_size], [c+2*vec_size+1, c+2*vec_size+1]])
            c += 2*vec_size+1
        end
    end
    return ind
end

function default_parameters(prec::preconditioner, ω::Real)
    if prec.name[1:4] == "full"
        parameters = Array{Float64, 1}(undef, 2^(prec.L+1)+prec.L-5+1)
        for l=1:(prec.L-1)
            N_ = 2^l+1
            h_ = 1/(N_-1)
            ψ, χ = BPX_weights(prec.L, l, prec.BCs)
            parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]] .= ψ
            parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]] .= χ
            parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]] .= 1
        end
    elseif prec.name[1:3] == "sym"
        parameters = Array{Float64, 1}(undef, 2^prec.L+prec.L-3+1)
        for l=1:(prec.L-1)
            ψ, _ = BPX_weights(prec.L, l, prec.BCs)
            parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]] .= ψ
            parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]] .= 1
        end
    elseif prec.name[1:8] == "reduced "
        parameters = Array{Float64, 1}(undef, 2*2^(prec.L-1)+prec.L-1+1); parameters .= 1
        ψ, χ  = BPX_weights(prec.L, 1, prec.BCs)
        parameters[1:2^(prec.L-1)] .= ψ
        parameters[(2^(prec.L-1)+1):2^(prec.L)] .= χ
    elseif prec.name[1:min(17, length(prec.name))] == "reduced+symmetric"
        parameters = Array{Float64, 1}(undef, 2^(prec.L-1)+prec.L-1+1); parameters .= 1
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs)
        parameters[1:2^(prec.L-1)] .= ψ
    elseif prec.name[1:7] == "2D full"
        parameters = Array{Float64, 1}(undef, 2^(prec.L+2)+prec.L-8)
        for l=1:(prec.L-1)
            ψ, χ = BPX_weights(prec.L, l, prec.BCs[1:2])
            parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]] .= ψ
            parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]] .= χ
            ψ, χ = BPX_weights(prec.L, l, prec.BCs[3:4])
            parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]] .= ψ
            parameters[prec.Ind[l][4][1]:prec.Ind[l][4][2]] .= χ
            if (prec.name[(end-2):end] == "BPX")
                parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]] .= 1/2^prec.L
            else
                parameters[prec.Ind[l][5][1]:prec.Ind[l][5][2]] .= 1
            end
        end
    elseif prec.name[1:6] == "2D sym"
        parameters = Array{Float64, 1}(undef, 2^(prec.L+1)+prec.L-4)
        for l=1:(prec.L-1)
            ψ, _ = BPX_weights(prec.L, l, prec.BCs[1:2])
            parameters[prec.Ind[l][1][1]:prec.Ind[l][1][2]] .= ψ
            ψ, _ = BPX_weights(prec.L, l, prec.BCs[3:4])
            parameters[prec.Ind[l][2][1]:prec.Ind[l][2][2]] .= ψ
            if (prec.name[(end-2):end] == "BPX")
                parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]] .= 1/2^prec.L
            else
                parameters[prec.Ind[l][3][1]:prec.Ind[l][3][2]] .= 1
            end
        end
    elseif prec.name[1:11] == "2D reduced "
        parameters = Array{Float64, 1}(undef, 2^(prec.L+1)+prec.L-1+1)
        if (prec.name[(end-2):end] == "BPX")
            parameters .=  1/2^prec.L
        else
            parameters .=  1
        end
        ψ, χ  = BPX_weights(prec.L, 1, prec.BCs[1:2])
        parameters[1:2^(prec.L-1)] .= ψ
        parameters[(2^(prec.L-1)+1):2^(prec.L)] .= χ
        ψ, χ  = BPX_weights(prec.L, 1, prec.BCs[3:4])
        parameters[2^(prec.L)+1:(3*2^(prec.L-1))] .= ψ
        parameters[(3*2^(prec.L-1)+1):2^(prec.L+1)] .= χ
    elseif prec.name[1:min(20, length(prec.name))] == "2D reduced+symmetric"
        parameters = Array{Float64, 1}(undef, 2^prec.L+prec.L-1+1)
        if (prec.name[(end-2):end] == "BPX")
            for l = 1:(prec.L-1)
                N = 2^prec.L+1
                h = 1/(N-1)
                N_ = 2^(prec.L-l)+1
                h_ = 1/(N_-1)
                parameters[end-l] = (h/h_)
            end
        else
            parameters .=  1
        end
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[1:2])
        parameters[1:2^(prec.L-1)] .= ψ
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[3:4])
        parameters[(2^(prec.L-1)+1):2^(prec.L)] .= ψ
    elseif prec.name[1:14] == "2D tridiagonal"
        parameters = Array{Float64, 1}(undef, 6*(prec.L-1)+2^(prec.L)+prec.L)
        parameters .=  1/2^prec.L
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[1:2])
        parameters[1:2^(prec.L-1)] .= ψ
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[3:4])
        parameters[(2^(prec.L-1)+1):2^(prec.L)] .= ψ
        parameters[(2^(prec.L)+1):6:(6*(prec.L-1)+2^(prec.L))] .= 1
        parameters[(2^(prec.L)+2):6:(6*(prec.L-1)+2^(prec.L))] .= 1e-4
        parameters[(2^(prec.L)+3):6:(6*(prec.L-1)+2^(prec.L))] .= 1e-4
        parameters[(2^(prec.L)+4):6:(6*(prec.L-1)+2^(prec.L))] .= 1
        parameters[(2^(prec.L)+5):6:(6*(prec.L-1)+2^(prec.L))] .= 1e-4
        parameters[(2^(prec.L)+6):6:(6*(prec.L-1)+2^(prec.L))] .= 1e-4
    elseif prec.name == "2D double Schwarz"
        parameters = Array{Float64, 1}(undef, 2^(prec.L+1)+prec.L)
        parameters .=  1
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[1:2])
        parameters[1:2^(prec.L-1)] .= ψ
        parameters[(1+2^(prec.L)):(2^(prec.L)+2^(prec.L-1))] .= ψ
        ψ, _  = BPX_weights(prec.L, 1, prec.BCs[3:4])
        parameters[(2^(prec.L-1)+1):2^(prec.L)] .= ψ
        parameters[(2^(prec.L)+2^(prec.L-1)+1):2^(prec.L+1)] .= ψ
    end
    parameters[end] = ω
    return parameters
end

function default_parameters(T_prec::T_preconditioner, ω::Real)
    default_parameters(transform_T(T_prec), ω)
end

function matching_preconditioner(prec::preconditioner)
    if (prec.name[1:2] != "2D")
        if (prec.name[(end-2):end] == "BPX")
            B = BPX_matrix(prec.L, prec.BCs)
        elseif prec.name[(end-6):end] == "Schwarz"
            B = Schwarz_matrix(prec.A, prec.L, prec.BCs)
        end
    else
        if (prec.name[(end-2):end] == "BPX")
            B = BPX_matrix_2D(prec.L, prec.BCs)
        elseif prec.name[(end-6):end] == "Schwarz"
            B = Schwarz_matrix_2D(prec.A, prec.L, prec.BCs)
        end
    end
    return B
end

function matching_preconditioner(T_prec::T_preconditioner)
    matching_preconditioner(transform_T(T_prec))
end

function initialise_preconditioner(A::SparseMatrixCSC{<:Real, <:Integer}, L::Integer, BCs::Array{String, 1}, name::String, symmetry::String)
    Ind = indices(BCs, L, name)
    prec = preconditioner(A, Ind, BCs, L, name, symmetry)
    parameters = default_parameters(prec, 1.0)
    B = assemble_matrix(parameters, prec)
    A_prec(x) = B*A*B*x
    if symmetry == "left"
        A_prec(x) = B*A*x
    end
    ω = 1/ρ_error(size(B)[1], 20, A_prec, 10)
    parameters[end] = ω
    return prec, B, parameters
end

function initialise_preconditioner(C::Array{Array{SparseMatrixCSC{Float64,Int64},1},1}, L::Integer, BCs::Array{String, 1}, name::String, symmetry::String)
    Ind = indices(BCs, L, name)
    A = sum(kron(c[1], c[2]) for c = C)
    prec = T_preconditioner(A, C, Ind, BCs, L, name, symmetry)
    B = matching_preconditioner(prec)
    A_prec(x) = B*A*B*x
    if symmetry == "left"
        A_prec(x) = B*A*x
    end
    ω = 1/ρ_error(size(B)[1], 20, A_prec, 10)
    parameters = default_parameters(prec, ω)
    return prec, B, parameters
end

function power_method(f::Function, D::Integer, N_sweeps::Integer)
    x = randn(D)
    α = 1
    for i = 1:N_sweeps
        x = f(x)
        α = maximum(abs.(x))
        x ./= α
    end
    return α
end

function measure_condition_number(parameters::AbstractArray{<:Real, 1}, prec::preconditioner)
    M = assemble_matrix(parameters::AbstractArray{<:Real, 1}, prec::preconditioner)
    if prec.symmetry == "symmetric"
        T = M*prec.A*M
    elseif prec.symmetry == "left"
        T = M*prec.A
    end
    σ = sort(abs.(eigvals(Array(T))))
    return σ[end]/σ[1]
end

function measure_condition_number(parameters::AbstractArray{<:Real, 1}, T_prec::T_preconditioner)
    measure_condition_number(parameters, transform_T(T_prec))
end

function measure_prec_condition_number(parameters::AbstractArray{<:Real, 1}, prec::preconditioner)
    M = matching_preconditioner(prec::preconditioner)
    if prec.symmetry == "symmetric"
        T = M*prec.A*M
    elseif prec.symmetry == "left"
        T = M*prec.A
    end
    σ = sort(abs.(eigvals(Array(T))))
    return σ[end]/σ[1]
end

function measure_prec_condition_number(parameters::AbstractArray{<:Real, 1}, T_prec::T_preconditioner)
    measure_prec_condition_number(parameters, transform_T(T_prec))
end

function measure_s_condition_number(parameters::AbstractArray{<:Real, 1}, prec::preconditioner, params::optimization_parameters)
    ρ = spectral_radius(parameters, prec, params)
    M = assemble_matrix(parameters, prec)
    α = scalar_loss_function([1/ρ, ], M, prec, params)
    return 1/(1-α)
end

function measure_s_condition_number(parameters::AbstractArray{<:Real, 1}, T_prec::T_preconditioner, params::optimization_parameters)
    measure_s_condition_number(parameters, transform_T(T_prec), params)
end

function test_matching_prec(equation::Function, parameters::Array{<:Real, 1}, BCs::Array{String, 1}, L::Integer, name::String, symmetry::String, tol::Real=1e-11)
    A = equation(L, BCs, parameters)
    prec = preconditioner(A, [[[1, 2], [3, 4], [5, 6]]], BCs, L, name, symmetry)
    M = matching_preconditioner(prec)
    if prec.symmetry == "symmetric"
        vals = svdvals(Array(M*prec.A*M))
    elseif prec.symmetry == "left"
        vals = svdvals(Array(M*prec.A))
    end
    prec_condition_number = vals[1]/vals[end-(vals[end]<=tol)]
    vals = svdvals(Array(prec.A))
    condition_number = vals[1]/vals[end-(vals[end]<=tol)]
    return condition_number, prec_condition_number
end

include("equations.jl")
include("tensor_equations.jl")

function test_T_precs(equation::Function, eq_params::Array{<:Real, 1},  BCs::Array{String, 1}, L::Integer, name::String, symmetry::String)
    C = equation(L, BCs, eq_params)
    T_prec, B, parameters = initialise_T_preconditioner(C, L, BCs, name, symmetry)
    N = 2^L+1
    DOF_x, DOF_y = N - 2 + sum(BCs[1:2] .== "N"), N - 2 + sum(BCs[3:4] .== "N")
    x, y = randn(DOF_x*DOF_y, 1), zeros(DOF_x*DOF_y, 1)
    z = transpose(reshape(x, (DOF_x, DOF_y)))
    y = B*x
    z = apply_T_preconditioner(z, parameters, T_prec)
    z = reshape(transpose(z), (DOF_x*DOF_y, 1))
    println(norm(z-y))
end

function test_precs(equation::Function, eq_params::Array{<:Real, 1},  BCs::Array{String, 1}, L::Integer, name::String, symmetry::String)
    A = equation(L, BCs, eq_params)
    prec, B, parameters = initialise_preconditioner(A, L, BCs, name, symmetry)
    B_prec = assemble_matrix(parameters, prec)
    N = 2^L+1
    DOF_x, DOF_y = N - 2 + sum(BCs[1:2] .== "N"), N - 2 + sum(BCs[3:4] .== "N")
    x, y, z = randn(DOF_x*DOF_y, 1), zeros(DOF_x*DOF_y, 1), zeros(DOF_x*DOF_y, 1)
    z = B_prec*x
    y = B*x
    println(norm(z-y))
end

function test()
    # "2D full BPX"
    # "2D reduced BPX"
    # "2D symmetric BPX"
    # "2D reduced+symmetric BPX"
    for BCs = Iterators.product(["D", "N"], ["D", "N"], ["D", "N"], ["D", "N"])
        println(BCs)
        test_T_precs(T_Laplace_2D, [10.0, 10.0, ], [BCs...], 5, "2D full Schwarz", "left")
        #test_precs(Anisotropic_Laplace_2D, [1.0, 1.0, ], [BCs...], 4, "2D full BPX", "left")
        println("=========")
    end
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
using Documenter, DualMatrixTools, LinearAlgebra

makedocs(
    sitename="DualMatrixTools Documentation",
    # options
    modules = [DualMatrixTools]
)

deploydocs(
    repo = "github.com/briochemc/DualMatrixTools.jl.git",
)


text:
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
using SemanticModels.CategoryTheory
import SemanticModels.CategoryTheory: undecorate, ⊔
using SemanticModels.PetriModels
using SemanticModels.PetriCospans
import SemanticModels: model

"""    fluxes(m::Petri.Model)

a PetriNet interpreter that computes the mass action kinetics of a petri net.

Usage: pass this to ODEProblem to set up an ODE for a given model.
"""
function fluxes(m::Petri.Model)
    S = m.S
    T = m.Δ
    nS = length(S)
    nT = length(T)
    ϕ = zeros(Float64, nT)
    f(du, u, p, t) = begin
        for (i, t) in enumerate(T)
            ins = t[1]
            # TODO: accomodate multiplicites here
            ϕ[i] = p[i]*prod(u[ins])
        end
        for i in S
            du[i] = 0
        end
        for (i, t) in enumerate(T)
            ins = t[1]
            out = t[2]
            for s in ins
                # TODO: accomodate multiplicites here
                du[s] -= ϕ[i]
            end
            for s in out
                # TODO: accomodate multiplicites here
                du[s] += ϕ[i]
            end
        end
        return du
    end
    return f
end

u₀(m::Petri.Model) = begin
    zeros(Float64, length(m.S))
end

u₀(m::Petri.Model, initialS, initialI=1) = begin
    u0=zeros(Float64, length(m.S))
    u0[1] = initialS
    u0[2] = initialI
    return u0
end

function savedata(sol::ODESolution, colnames, fname::String)
    open(fname, "w") do fp
        println(fp, "time, $colnames")
        map(tuples(sol)) do (u,t)
            print(fp, "$t")
            map(u) do x
                print(fp, ",$x")
            end
            println(fp, "")
        end
    end
end

model(c::PetriCospan) = left(c.f).d[1].model

X = FinSet(1)
# Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
# m_seir = left(Fseir.f).d[1]

Psir = PetriModel(
          Petri.Model(1:3,[
            ([1,2],[2,2]), # exposure
            ([2],[3]),     # recovery
            ], missing, missing))

m = Psir.model

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[2]  = 1

β = [10/sum(u0), 1/5]

tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())

savedata(sol, "S,I,R", "sirdata.csv")

Psird = PetriModel(
          Petri.Model(1:3,[
            ([1,2],[2,2]), # exposure
            ([2],[3]),     # recovery
            ([2],[4]),     # death
            ], missing, missing))

m = Psir.model

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[2]  = 1

β = [10/sum(u0), 1/5, 1/10]

tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())

savedata(sol, "S,I,R,D", "sirddata.csv")

Pseir = PetriModel(
          Petri.Model(1:5,[
            ([1,2],[3,2]), # exposure
            ([3],[2]),     # onset
            ([2],[4]),     # recovery
            ], missing, missing))

m = Pseir.model

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[2]  = 1

β = [10/sum(u0), 1/2, 1/5]

tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())

savedata(sol, "S,I,E,R", "seirdata.csv")
Pseird = PetriModel(
          Petri.Model(1:5,[
            ([1,2],[3,2]), # exposure
            ([3],[2]),     # onset
            ([2],[4]),     # recovery
            ([2],[5]),     # death
            ], missing, missing))

m = Pseird.model

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[2]  = 1

seirdparams() = begin
    βseird = [10/sum(u0), 1/2, 1/5, 1/16]
end

tspan = (0,100.0)
prob = ODEProblem(fluxes(m), u0, tspan, seirdparams())
sol = OrdinaryDiffEq.solve(prob, alg=Tsit5())

savedata(sol, "S,I,E,R,D", "seirddata.csv")


text:
using RecipesBase
@recipe function f(X::Trajectory{<:AbstractVector,<:AbstractVector})
  X.t, X.x
end


text:
@testset "Delete DB" begin
  println(testname("Delete DB tests"))
  print("  [  ] Delete test database: $database ")
  result = deletedb(cl, database)
  @test result["ok"] == true
  println("\r  [OK] Delete test database: $database")
end

max_stars_repo_name: xpqz/couchzilla

max_issues_repo_name: xpqz/couchzilla

max_forks_repo_name: xpqz/couchzilla


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

using Markdown
using InteractiveUtils

# ╔═╡ 786a41c2-63ea-11ec-3ec3-473864f1051d
using Pkg

# ╔═╡ 6f9c2a09-371e-4c5a-a9d3-acb53ea222ac
Pkg.activate(".")

# ╔═╡ 2106f13d-f759-444a-8c2e-3bfd2bcb3478
Pkg.add("CSV")

# ╔═╡ f58a6aa0-30f7-4362-bf67-64fbd3376ed7
using HTTP, CSV, DataFrames

# ╔═╡ 152d990a-7de0-4ecc-b223-0934d579281d
using Images

# ╔═╡ 520aacd9-5a6d-40bb-996d-026bf9499962
csv_reader = CSV.File("./runs/dataset/test.csv")

# ╔═╡ 66fee52c-5f97-4840-8609-d0ec4d7f7fdc
dataall = csv_reader |> Tables.matrix

# ╔═╡ baf17f10-8033-49a7-a80a-8f0e5f36e6de
dataall[1,:]

# ╔═╡ d56cbd95-5ef7-4457-ad81-333174f63fe9
images = map(i-> Gray.(reshape(dataall[i,:], (28,28))./255), 1:100)

# ╔═╡ fc98a3f0-6275-4078-b568-d74868df8bb2
r = HTTP.post("http://localhost:8001/digitreg", [], string(dataall[5,:]))

# ╔═╡ 895b23e7-2319-4920-8a9f-ec31735682d2
 a = String(r.body)

# ╔═╡ 5bd3dffc-f152-411d-9567-107a64327868
parse(Int64, a)

# ╔═╡ Cell order:
# ╠═786a41c2-63ea-11ec-3ec3-473864f1051d
# ╠═6f9c2a09-371e-4c5a-a9d3-acb53ea222ac
# ╠═2106f13d-f759-444a-8c2e-3bfd2bcb3478
# ╠═f58a6aa0-30f7-4362-bf67-64fbd3376ed7
# ╠═520aacd9-5a6d-40bb-996d-026bf9499962
# ╠═66fee52c-5f97-4840-8609-d0ec4d7f7fdc
# ╠═baf17f10-8033-49a7-a80a-8f0e5f36e6de
# ╠═152d990a-7de0-4ecc-b223-0934d579281d
# ╠═d56cbd95-5ef7-4457-ad81-333174f63fe9
# ╠═fc98a3f0-6275-4078-b568-d74868df8bb2
# ╠═895b23e7-2319-4920-8a9f-ec31735682d2
# ╠═5bd3dffc-f152-411d-9567-107a64327868

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
function tcm_predictMNL(N, u, assortments)
	M = size(assortments,1)

	v = zeros(Float64,size(assortments))

    for m = 1 : M
        v[m,N] = 1.0
        for p = 1 : N-1
            v[m,p] = exp(u[p]) * assortments[m,p]
        end
        v[m,:] = v[m,:] / sum(v[m,:])
    end

    return v
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
module Fields

using Gridap.Arrays
import Gridap.Arrays: testvalue
import Gridap.Arrays: inverse_map
import Gridap.Arrays: get_children
import Gridap.Arrays: testitem

using Gridap.Helpers: @abstractmethod, @notimplemented
using Gridap.Helpers: @notimplementedif, @unreachable, @check

using Gridap.Algebra: mul!

using Gridap.TensorValues
using Gridap.Algebra

using LinearAlgebra: mul!, Transpose, diag

using ForwardDiff
using FillArrays
using NLsolve
using Test
using StaticArrays
using LinearAlgebra

import LinearAlgebra: det, inv, transpose, tr, cross
import LinearAlgebra: ⋅, dot

import Base: +, -, *, /
import Gridap.TensorValues: ⊗, ⊙, symmetric_part, outer, meas

import Gridap.Arrays: IndexStyle
import Gridap.Arrays: return_cache
import Gridap.Arrays: return_type
import Gridap.Arrays: testargs
import Gridap.Arrays: return_value
import Gridap.Arrays: evaluate!
import Gridap.Arrays: lazy_map
import Gridap.Arrays: array_cache

export evaluate
export evaluate!
export return_type
export return_cache

export Field
export GenericField
export ConstantField
export constant_field
export FieldGradient
export FieldGradientArray
export ZeroField
export MockField
export MockFieldArray
export Point
export inverse_map

export AffineMap
export affine_map

export gradient
export ∇
export ∇∇
export gradient_type
export push_∇
export push_∇∇
export pinvJt

export curl
export grad2curl
export laplacian
export divergence
export DIV
export Δ
export ε
export symmetric_gradient

export test_field
export test_field_array

export linear_combination
export integrate
export IntegrationMap

export ArrayBlock
export VectorBlock
export MatrixBlock
export BlockMap

export VoidField
export VoidFieldMap
export VoidBasis
export VoidBasisMap

export DensifyInnerMostBlockLevelMap

include("FieldsInterfaces.jl")

include("FieldArrays.jl")

include("MockFields.jl")

include("AffineMaps.jl")

include("ApplyOptimizations.jl")

include("DiffOperators.jl")

include("AutoDiff.jl")

include("ArrayBlocks.jl")

include("InverseFields.jl")

include("DensifyInnerMostBlockLevelMaps.jl")

end


text:
using Documenter, PWS

makedocs(sitename="PWS.jl")

deploydocs(
    repo = "github.com/manuel-rhdt/PWS.jl.git",
)


text:
module PhotonTest

@info "Running Photon Unit Tests"

using Photon
using Test

@testset "Photon" begin

@info "Testing Basics"
include("./basic.jl")

@info "Testing Dense Layers"
include("./dense.jl")

@info "Testing Convolutional Layers"
include("./conv.jl")

@info "Testing Recurrent Layers"
include("./recurrent.jl")

@info "Testing Large models"
include("./modelzoo.jl")

@info "Testing Training"
include("./training.jl")
include("./optimise.jl")
include("./loss.jl")

@info "Testing Metrics"
include("./metrics.jl")

@info "Testing Callbacks"
include("./callbacks.jl")

@info "Testing Complex Scenarios"
include("./complex.jl")

@info "Testing dataset and dataloader"
include("./data.jl")

end

end