text:
# deps.jl
using Luxor, Colors, FileIO, ImageView
using ColorVectorSpace

type primitive
    center::Point
    majlen::Float64
    minlen::Float64
    angle::Float64
    name::String
end

import Base.+
function Base.:+(a::primitive, b::primitive)
    center = a.center + b.center
    maj = a.majlen + b.majlen
    min = a.minlen + b.minlen
    ang = a.angle + b.angle

    primitive(center, maj, min, ang, a.name)
end

import Base.copy
function Base.:copy(a::primitive)
    b = primitive(a.center, a.majlen, a.minlen, a.angle, a.name)
    return b
end
# Define helper functions

# Load float image
function loadFloatImage(path)
    # return float(channelview(load(path)))
    return load(path)
end

# Load bit (binary) image
function loadBitImage(path)
    return (Images.Gray.(load(path)) .> 0.5)
end

function genRandShape(shape::String, bounds::Vector)
    cvsdiag = norm(bounds)
    center = randompoint(1, 1, bounds[1], bounds[2]) # needs to be an Int
    majlen = rand(1:.1:cvsdiag) # granularity of .1 pixel
    minlen = clamp(rand(1:.1:cvsdiag), 10, majlen)
    ang = rand(1:.05:2*π) # granularity of .05 radians

    return primitive(center, majlen, minlen, ang, shape)
end

function makeShape(shape::primitive, canvas_size::Vector)
    # Take in a shape and a tuple of canvas size
    # Return a canvas sized mask with the full shape masked out
    if shape.name == "Triangle"
        return triangleRasterize(shape, canvas_size)
    elseif shape.name == "Rectangle"
        return rectangleRasterize(shape, canvas_size)
    elseif shape.name == "Ellipse"
        return ellipseRasterize(shape, canvas_size)
    elseif shape.name == "Curve"
        return curveRasterize(shape, canvas_size)
    end
end

function applyMask(canvas, shape::primitive, src_img)
    # Take in a canvas, a shape, and the source image.
    m, n = size(src_img)
    # Figure out what pixels on the canvas are under the shape given parameters
    mask = makeShape(shape, [m,n])
    # calculate fill color as a weighted sum of the colors of pixels in
    # the source image that lie within the shape.
    color = mean(src_img[Bool.(mask)])
    # Set the color of `canvas` under the mask to the value `color`
    # extract the mean color under the mask
    offset = mean(canvas[Bool.(mask)])
    c_fill = color - offset
    masked_canvas = copy(canvas)
    return masked_canvas += c_fill .* Float64.(mask)
end

function MSE(image, canvas)
    # Take in an image and a canvas
    # Return the Variance between the canvas and the image
    h = (image - canvas)[:]
    return sum([red.(h)'*red.(h), green.(h)'*green.(h), blue.(h)'*blue.(h)])
end

function mutateShape(shape::primitive)
    # Take in shape, randomly mutate, and return
    # granularity
    g = 16

    c_offset = Point(rand(-g:g), rand(-g:g))
    maj_offset = rand(-g:g)
    min_offset = rand(-g:g)
    ang_offset = 2*π*rand(0:2*g)/360

    offset = primitive(c_offset, maj_offset, min_offset, ang_offset, shape.name)
    return shape + offset
end

# Trinagle Specific Functions
function triangleRasterize(tri::primitive, bounds::Vector)
    Drawing(bounds[2], bounds[1], "mask.png")
    background("black")
    sethue("white")

    v0, v1, v2 = ngon(tri.center, tri.minlen, 3, tri.angle, vertices=true)
    poly([v0 + tri.majlen - tri.minlen, v1, v2], :fill)
    finish()
    return loadBitImage("mask.png")
end

# Rectangle Specific Functions
function rectangleRasterize(rec::primitive, bounds::Vector)
    Drawing(bounds[2], bounds[1], "mask.png")
    background("black")
    sethue("white")

    rotate(rec.angle)
    translate(rec.center)
    rect(O, rec.majlen, rec.minlen, :fill)
    finish()
    return loadBitImage("mask.png")
end

# Ellipse Specific Functions
function ellipseRasterize(ellip::primitive, bounds::Vector)
    Drawing(bounds[2], bounds[1], "mask.png")
    background("black")
    sethue("white")

    rotate(ellip.angle)
    translate(ellip.center)
    ellipse(O, ellip.majlen, ellip.minlen, :fill)
    finish()
    return loadBitImage("mask.png")
end

# Curve Specific Functions
function curveRasterize(cve::primitive, bounds::Vector)
    Drawing(bounds[2], bounds[1], "mask.png")
    background("black")

    v0, v1, v2 = ngon(cve.center, cve.majlen, 3, cve.angle, vertices=true)
    sethue("white")
    setline(3)
    arc2r(v0, v1, v2, :stroke)
    finish()
    return loadBitImage("mask.png")
end

# Hill Climbing
function hillClimb(canvas, shape::primitive, image, max_age::Int)
    best_canvas = copy(canvas)
    best_error = MSE(image, best_canvas)
    step = 0

    for age = 0:max_age
        new_shape = mutateShape(shape)
        new_canvas = applyMask(canvas, new_shape, image)
        new_error = MSE(image, new_canvas)
        if new_error < best_error
            best_error = copy(new_error)
            best_canvas = copy(new_canvas)
            age = -1
            shape = copy(new_shape)
        end
        step += 1
    end
    return best_canvas
end

max_stars_repo_name: mcbaron/primitive

max_issues_repo_name: mcbaron/primitive

max_forks_repo_name: mcbaron/primitive


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

# Classification of MNIST dataset
## with the convolutional neural network known as LeNet5.
## This script also combines various
## packages from the Julia ecosystem with Flux.
using Flux, Statistics
using Flux.Data: DataLoader
using Flux: onehotbatch, onecold, @epochs
using Flux.Losses: logitcrossentropy
using Base: @kwdef
using CUDA
using MLDatasets
using OperatorFlux

function getdata(args, device)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    # Loading Dataset
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    # Reshape Data in order to flatten each image into a linear array
    xtrain = reshape(xtrain, 1, 28, 28, :)
    xtest = reshape(xtest, 1, 28, 28, :)

    # One-hot-encode the labels
    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    # Create DataLoaders (mini-batch iterators)
    train_loader =
        DataLoader((xtrain, ytrain), batchsize = args.batchsize, shuffle = true)
    test_loader = DataLoader((xtest, ytest), batchsize = args.batchsize)

    return train_loader, test_loader
end

function build_model(; nclasses = 10)
    trafo = FourierTransform(modes = (8, 8))
    return Chain(
        SpectralKernelOperator(trafo, 1 => 1, gelu),
        SpectralKernelOperator(trafo, 1 => 1, gelu),
        flatten,
        Dense(784, nclasses),
    )
end

function loss_and_accuracy(data_loader, model, device)
    acc = 0
    ls = 0.0f0
    num = 0
    for (x, y) in data_loader
        x, y = device(x), device(y)
        ŷ = model(x)
        ls += logitcrossentropy(ŷ, y, agg = sum)
        acc += sum(onecold(ŷ) .== onecold(y))
        num += size(x)[end]
    end
    return ls / num, acc / num
end

@kwdef mutable struct Args
    η::Float64 = 3e-4       # learning rate
    batchsize::Int = 256    # batch size
    epochs::Int = 10        # number of epochs
    use_cuda::Bool = true   # use gpu (if cuda available)
end

function train(; kws...)
    args = Args(; kws...) # collect options in a struct for convenience

    if CUDA.functional() && args.use_cuda
        @info "Training on CUDA GPU"
        CUDA.allowscalar(false)
        device = gpu
    else
        @info "Training on CPU"
        device = cpu
    end

    # Create test and train dataloaders
    train_loader, test_loader = getdata(args, device)

    # Construct model
    model = build_model() |> device
    ps = Flux.params(model) # model's trainable parameters

    ## Optimizer
    opt = ADAM(args.η)

    ## Training
    for epoch in 1:(args.epochs)
        for (x, y) in train_loader
            x, y = device(x), device(y) # transfer data to device
            gs = gradient(() -> logitcrossentropy(model(x), y), ps) # compute gradient
            Flux.Optimise.update!(opt, ps, gs) # update parameters
        end

        # Report on train and test
        train_loss, train_acc = loss_and_accuracy(train_loader, model, device)
        test_loss, test_acc = loss_and_accuracy(test_loader, model, device)
        println("Epoch=$epoch")
        println("  train_loss = $train_loss, train_accuracy = $train_acc")
        println("  test_loss = $test_loss, test_accuracy = $test_acc")
    end
end

### Run training
if abspath(PROGRAM_FILE) == @__FILE__
    train()
end
# train(η=0.01) # can change hyperparameters

ext: jl


text:
module YAXArrays
global const YAXDefaults = (
    workdir = Ref("./"),
    recal = Ref(false),
    chunksize = Ref{Any}(:input),
    max_cache = Ref(1e8),
    cubedir = Ref(""),
    subsetextensions = [],
)
global const workdir = YAXDefaults.workdir
global const recal = YAXDefaults.recal
function __init__()
    YAXDefaults.workdir[] = get(ENV, "YAXARRAY_WORKDIR", "./")
    YAXDefaults.max_cache[] = parse(Float64, get(ENV, "YAXARRAY_MAX_CACHE", "100")) * 1e6
    YAXDefaults.cubedir[] = get(ENV, "YAXARRAY_CUBEDIR", "")
end
YAXdir(x::String) = YAXDefaults.workdir[] = x
recalculate(x::Bool) = YAXDefaults.recal[] = x
recalculate() = YAXDefaults.recal[]
YAXdir() = YAXDefaults.workdir[]
export YAXdir

include("YAXTools.jl")
include("Cubes/Cubes.jl")
include("DatasetAPI/Datasets.jl")
include("DAT/DAT.jl")

using Reexport: @reexport
using YAXArrayBase: getattributes

@reexport using Dates: Date, DateTime
@reexport using IntervalSets: (..)
@reexport using .Cubes
@reexport using .Cubes.Axes

@reexport using .DAT
@reexport using .Datasets
@reexport using .YAXTools: @loadOrGenerate # from YAXTools

@deprecate saveCube(data, filename) savecube(data, filename)
@deprecate loadCube(filename) loadcube(filename)
@deprecate rmCube(filename) rmcube(filename)
@deprecate exportcube(data, filename; kwargs...) savecube(
    data,
    filename;
    backend = :netcdf,
    kwargs...,
)
@deprecate cubeproperties(x) getattributes(x)
@deprecate concantenateCubes(args...) concatenatecubes(args...)

# include("precompile.jl")
# _precompile_()

end # module


text:
# Randomized Dynamic Mode Decomposition

function rdmd(X)
end


text:
module GatekeepHelperCustomColorHeartDoor

using ..Ahorn, Maple

@mapdef Entity "GatekeepHelper/CustomColorHeartDoor" CustomColorHeartDoor(x::Integer, y::Integer, width::Integer=40, height::Integer=Maple.defaultBlockHeight,
    requires::Integer=0, startHidden::Bool=false, color::String="18668F")

const placements = Ahorn.PlacementDict(
    "Custom Color Heart Door (GatekeepHelper)" => Ahorn.EntityPlacement(
        CustomColorHeartDoor,
        "rectangle"
    )
)

Ahorn.editingOptions(entity::CustomColorHeartDoor) = Dict{String, Any}(
)

Ahorn.minimumSize(entity::CustomColorHeartDoor) = 40, 8
Ahorn.resizable(entity::CustomColorHeartDoor) = true, false

Ahorn.nodeLimits(entity::CustomColorHeartDoor) = 0, 1

const heartPadding = 4

function Ahorn.selection(entity::CustomColorHeartDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)

    roomWidth, roomHeight = room.size
    width = get(entity.data, "width", 40)

    nodes = get(entity.data, "nodes", ())

    if isempty(nodes)
        return Ahorn.Rectangle(x, 0, width, roomHeight)

    else
        nx, ny = Int.(nodes[1])

        return [Ahorn.Rectangle(x, 0, width, roomHeight), Ahorn.Rectangle(nx - 8, ny, width + 16, 8)]
    end
end

function heartsWidth(heartSprite::Ahorn.Sprite, hearts::Integer)
    return hearts * (heartSprite.width + heartPadding) - heartPadding
end

function heartsPossible(width::Integer, edgeSprite::Ahorn.Sprite, heartSprite::Ahorn.Sprite, required::Integer)
    rowWidth = width - 2 * edgeSprite.width

    for i in 0:required
        if heartsWidth(heartSprite, i) > rowWidth
            return i - 1
        end
    end

    return required
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::CustomColorHeartDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])
        dy = ny - y

        Ahorn.drawRectangle(ctx, x, ny, width, 1, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
        Ahorn.drawRectangle(ctx, x, 2 * y - ny, width, 1, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))

        Ahorn.drawRectangle(ctx, nx - 8, ny, width + 16, 8, (1.0, 0.0, 0.0, 1.0), (0.0, 0.0, 0.0, 0.0))
    end
end

# Not completely accurate on heart positions, but good enough
function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CustomColorHeartDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 40))
    roomWidth, roomHeight = room.size
    hearts = Int(get(entity.data, "requires", 0))

    edgeSprite = Ahorn.getSprite("objects/heartdoor/edge", "Gameplay")
    heartSprite = Ahorn.getSprite("objects/heartdoor/icon00", "Gameplay")

    color = Ahorn.argb32ToRGBATuple(parse(Int, get(entity.data, "color", "18668F"), base=16))[1:3] ./ 255
    Ahorn.drawRectangle(ctx, x, 0, width, roomHeight, color, (0.0, 0.0, 0.0, 0.0))

    for i in 0:edgeSprite.height:roomHeight
        Ahorn.Cairo.save(ctx)

        Ahorn.drawImage(ctx, edgeSprite, x + width - edgeSprite.width, i)
        Ahorn.scale(ctx, -1, 1)
        Ahorn.drawImage(ctx, edgeSprite, -x - edgeSprite.width, i)

        Ahorn.Cairo.restore(ctx)
    end

    if hearts > 0
        fits = heartsPossible(width, edgeSprite, heartSprite, hearts)
        rows = ceil(Int, hearts / fits)

        for row in 1:rows
            fits = heartsPossible(width, edgeSprite, heartSprite, hearts)
            drawWidth = heartsWidth(heartSprite, fits)

            startX = x + round(Int, (width - drawWidth) / 2) + edgeSprite.width - 1
            startY = y - round(Int, rows / 2 * (heartSprite.height + heartPadding)) - heartPadding - 1

            for col in 1:fits
                drawX = (col - 1) * (heartSprite.width + heartPadding) - heartPadding
                drawY = row * (heartSprite.height + heartPadding) - heartPadding

                Ahorn.drawImage(ctx, heartSprite, startX + drawX, startY + drawY)
            end

            hearts -= fits
        end
    end
end

end


text:
using Test
using HaltonSequences

@testset "values" begin
    @test haltonvalue(1,2) ≈ 0.5
    @test haltonvalue(111111,2) ≈ 0.876182556152344
    @test haltonvalue(789,19) ≈ 0.534917626476163
    @test Halton{Rational}(17)[24] == 120//289
end

@testset "integrals" begin
    @test sum( Halton(5)[100:20100] )/20000 ≈ 0.5 atol=0.002
    p = HaltonPoint(4,length=32000)
    @test sum(p)/32000 ≈ [0.5,0.5,0.5,0.5] atol=0.001
    h = Halton{Float32}(29,start=1234,length=22000)
    @test sum( h.^2 )/22000 ≈ 1/3 atol=0.001
end


text:
#   Debug.Eval:
# ===============
# eval in debug scope

module Eval
using Debug.AST, Debug.Runtime, Debug.Graft
export debug_eval, Scope

debug_eval(scope::ModuleScope, ex) = scope.eval(ex)
function debug_eval(scope::LocalScope, ex)
    e = child(Expr(:let, ex), NoEnv()) # todo: actually wrap ex in a let?
    grafted = graft(e, scope, ex)

    assigned = setdiff(e.assigned, scope.env.assigned)
    if !isempty(e.defined)
        error("debug_eval: cannot define $(tuple(e.defined...)) in top scope")
    elseif !isempty(assigned)
        error("debug_eval: cannot assign $(tuple(assigned...)) in top scope")
    end

    eval = get_eval(scope)
    eval(Expr(:let, grafted))
end

end # module


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: greengrassv2
using AWS.Compat
using AWS.UUIDs

"""
    cancel_deployment(deployment_id)
    cancel_deployment(deployment_id, params::Dict{String,<:Any})

Cancels a deployment. This operation cancels the deployment for devices that haven't yet
received it. If a device already received the deployment, this operation doesn't change
anything for that device.

# Arguments
- `deployment_id`: The ID of the deployment.

"""
cancel_deployment(deploymentId; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/deployments/$(deploymentId)/cancel"; aws_config=aws_config)
cancel_deployment(deploymentId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/deployments/$(deploymentId)/cancel", params; aws_config=aws_config)

"""
    create_component_version()
    create_component_version(params::Dict{String,<:Any})

Creates a component. Components are software that run on AWS IoT Greengrass core devices.
After you develop and test a component on your core device, you can use this operation to
upload your component to AWS IoT Greengrass. Then, you can deploy the component to other
core devices. You can use this operation to do the following:    Create components from
recipes  Create a component from a recipe, which is a file that defines the component's
metadata, parameters, dependencies, lifecycle, artifacts, and platform capability. For more
information, see AWS IoT Greengrass component recipe reference in the AWS IoT Greengrass V2
Developer Guide. To create a component from a recipe, specify inlineRecipe when you call
this operation.    Create components from Lambda functions  Create a component from an AWS
Lambda function that runs on AWS IoT Greengrass. This creates a recipe and artifacts from
the Lambda function's deployment package. You can use this operation to migrate Lambda
functions from AWS IoT Greengrass V1 to AWS IoT Greengrass V2. This function only accepts
Lambda functions that use the following runtimes:   Python 2.7 – python2.7    Python 3.7
– python3.7    Python 3.8 – python3.8    Java 8 – java8    Node.js 10 – nodejs10.x
  Node.js 12 – nodejs12.x    To create a component from a Lambda function, specify
lambdaFunction when you call this operation.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"inlineRecipe"`: The recipe to use to create the component. The recipe defines the
  component's metadata, parameters, dependencies, lifecycle, artifacts, and platform
  compatibility. You must specify either inlineRecipe or lambdaFunction.
- `"lambdaFunction"`: The parameters to create a component from a Lambda function. You must
  specify either inlineRecipe or lambdaFunction.
- `"tags"`: A list of key-value pairs that contain metadata for the resource. For more
  information, see Tag your resources in the AWS IoT Greengrass V2 Developer Guide.
"""
create_component_version(; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/createComponentVersion"; aws_config=aws_config)
create_component_version(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/createComponentVersion", params; aws_config=aws_config)

"""
    create_deployment(target_arn)
    create_deployment(target_arn, params::Dict{String,<:Any})

Creates a continuous deployment for a target, which is a AWS IoT Greengrass core device or
group of core devices. When you add a new core device to a group of core devices that has a
deployment, AWS IoT Greengrass deploys that group's deployment to the new device. You can
define one deployment for each target. When you create a new deployment for a target that
has an existing deployment, you replace the previous deployment. AWS IoT Greengrass applies
the new deployment to the target devices. Every deployment has a revision number that
indicates how many deployment revisions you define for a target. Use this operation to
create a new revision of an existing deployment. This operation returns the revision number
of the new deployment when you create it. For more information, see the Create deployments
in the AWS IoT Greengrass V2 Developer Guide.

# Arguments
- `target_arn`: The ARN of the target AWS IoT thing or thing group.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"components"`: The components to deploy. This is a dictionary, where each key is the
  name of a component, and each key's value is the version and configuration to deploy for
  that component.
- `"deploymentName"`: The name of the deployment. You can create deployments without names.
  If you create a deployment without a name, the AWS IoT Greengrass V2 console shows the
  deployment name as <targetType>:<targetName>, where targetType and targetName
  are the type and name of the deployment target.
- `"deploymentPolicies"`: The deployment policies for the deployment. These policies define
  how the deployment updates components and handles failure.
- `"iotJobConfiguration"`: The job configuration for the deployment configuration. The job
  configuration specifies the rollout, timeout, and stop configurations for the deployment
  configuration.
- `"tags"`: A list of key-value pairs that contain metadata for the resource. For more
  information, see Tag your resources in the AWS IoT Greengrass V2 Developer Guide.
"""
create_deployment(targetArn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/deployments", Dict{String, Any}("targetArn"=>targetArn); aws_config=aws_config)
create_deployment(targetArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/deployments", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("targetArn"=>targetArn), params)); aws_config=aws_config)

"""
    delete_component(arn)
    delete_component(arn, params::Dict{String,<:Any})

Deletes a version of a component from AWS IoT Greengrass.  This operation deletes the
component's recipe and artifacts. As a result, deployments that refer to this component
version will fail. If you have deployments that use this component version, you can remove
the component from the deployment or update the deployment to use a valid version.

# Arguments
- `arn`: The ARN of the component version.

"""
delete_component(arn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/greengrass/v2/components/$(arn)"; aws_config=aws_config)
delete_component(arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/greengrass/v2/components/$(arn)", params; aws_config=aws_config)

"""
    delete_core_device(core_device_thing_name)
    delete_core_device(core_device_thing_name, params::Dict{String,<:Any})

Deletes a AWS IoT Greengrass core device, which is an AWS IoT thing. This operation removes
the core device from the list of core devices. This operation doesn't delete the AWS IoT
thing. For more information about how to delete the AWS IoT thing, see DeleteThing in the
AWS IoT API Reference.

# Arguments
- `core_device_thing_name`: The name of the core device. This is also the name of the AWS
  IoT thing.

"""
delete_core_device(coreDeviceThingName; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/greengrass/v2/coreDevices/$(coreDeviceThingName)"; aws_config=aws_config)
delete_core_device(coreDeviceThingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/greengrass/v2/coreDevices/$(coreDeviceThingName)", params; aws_config=aws_config)

"""
    describe_component(arn)
    describe_component(arn, params::Dict{String,<:Any})

Retrieves metadata for a version of a component.

# Arguments
- `arn`: The ARN of the component version.

"""
describe_component(arn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/metadata"; aws_config=aws_config)
describe_component(arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/metadata", params; aws_config=aws_config)

"""
    get_component(arn)
    get_component(arn, params::Dict{String,<:Any})

Gets the recipe for a version of a component. Core devices can call this operation to
identify the artifacts and requirements to install a component.

# Arguments
- `arn`: The ARN of the component version.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"recipeOutputFormat"`: The format of the recipe.
"""
get_component(arn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)"; aws_config=aws_config)
get_component(arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)", params; aws_config=aws_config)

"""
    get_component_version_artifact(arn, artifact_name)
    get_component_version_artifact(arn, artifact_name, params::Dict{String,<:Any})

Gets the pre-signed URL to download a public component artifact. Core devices call this
operation to identify the URL that they can use to download an artifact to install.

# Arguments
- `arn`: The ARN of the component version. Specify the ARN of a public component version.
- `artifact_name`: The name of the artifact. You can use the GetComponent operation to
  download the component recipe, which includes the URI of the artifact. The artifact name is
  the section of the URI after the scheme. For example, in the artifact URI
  greengrass:SomeArtifact.zip, the artifact name is SomeArtifact.zip.

"""
get_component_version_artifact(arn, artifactName; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/artifacts/$(artifactName)"; aws_config=aws_config)
get_component_version_artifact(arn, artifactName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/artifacts/$(artifactName)", params; aws_config=aws_config)

"""
    get_core_device(core_device_thing_name)
    get_core_device(core_device_thing_name, params::Dict{String,<:Any})

Retrieves metadata for a AWS IoT Greengrass core device.

# Arguments
- `core_device_thing_name`: The name of the core device. This is also the name of the AWS
  IoT thing.

"""
get_core_device(coreDeviceThingName; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)"; aws_config=aws_config)
get_core_device(coreDeviceThingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)", params; aws_config=aws_config)

"""
    get_deployment(deployment_id)
    get_deployment(deployment_id, params::Dict{String,<:Any})

Gets a deployment. Deployments define the components that run on AWS IoT Greengrass core
devices.

# Arguments
- `deployment_id`: The ID of the deployment.

"""
get_deployment(deploymentId; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/deployments/$(deploymentId)"; aws_config=aws_config)
get_deployment(deploymentId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/deployments/$(deploymentId)", params; aws_config=aws_config)

"""
    list_component_versions(arn)
    list_component_versions(arn, params::Dict{String,<:Any})

Retrieves a paginated list of all versions for a component.

# Arguments
- `arn`: The ARN of the component version.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
"""
list_component_versions(arn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/versions"; aws_config=aws_config)
list_component_versions(arn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components/$(arn)/versions", params; aws_config=aws_config)

"""
    list_components()
    list_components(params::Dict{String,<:Any})

Retrieves a paginated list of component summaries. This list includes components that you
have permission to view.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
- `"scope"`: The scope of the components to list. Default: PRIVATE
"""
list_components(; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components"; aws_config=aws_config)
list_components(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/components", params; aws_config=aws_config)

"""
    list_core_devices()
    list_core_devices(params::Dict{String,<:Any})

Retrieves a paginated list of AWS IoT Greengrass core devices.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
- `"status"`: The core device status by which to filter. If you specify this parameter, the
  list includes only core devices that have this status. Choose one of the following options:
     HEALTHY – The AWS IoT Greengrass Core software and all components run on the core
  device without issue.    UNHEALTHY – The AWS IoT Greengrass Core software or a component
  is in a failed state on the core device.
- `"thingGroupArn"`: The ARN of the AWS IoT thing group by which to filter. If you specify
  this parameter, the list includes only core devices that are members of this thing group.
"""
list_core_devices(; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices"; aws_config=aws_config)
list_core_devices(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices", params; aws_config=aws_config)

"""
    list_deployments()
    list_deployments(params::Dict{String,<:Any})

Retrieves a paginated list of deployments.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"historyFilter"`: The filter for the list of deployments. Choose one of the following
  options:    ALL – The list includes all deployments.    LATEST_ONLY – The list includes
  only the latest revision of each deployment.   Default: LATEST_ONLY
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
- `"targetArn"`: The ARN of the target AWS IoT thing or thing group.
"""
list_deployments(; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/deployments"; aws_config=aws_config)
list_deployments(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/deployments", params; aws_config=aws_config)

"""
    list_effective_deployments(core_device_thing_name)
    list_effective_deployments(core_device_thing_name, params::Dict{String,<:Any})

Retrieves a paginated list of deployment jobs that AWS IoT Greengrass sends to AWS IoT
Greengrass core devices.

# Arguments
- `core_device_thing_name`: The name of the core device. This is also the name of the AWS
  IoT thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
"""
list_effective_deployments(coreDeviceThingName; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)/effectiveDeployments"; aws_config=aws_config)
list_effective_deployments(coreDeviceThingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)/effectiveDeployments", params; aws_config=aws_config)

"""
    list_installed_components(core_device_thing_name)
    list_installed_components(core_device_thing_name, params::Dict{String,<:Any})

Retrieves a paginated list of the components that a AWS IoT Greengrass core device runs.

# Arguments
- `core_device_thing_name`: The name of the core device. This is also the name of the AWS
  IoT thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per paginated request.
- `"nextToken"`: The token to be used for the next set of paginated results.
"""
list_installed_components(coreDeviceThingName; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)/installedComponents"; aws_config=aws_config)
list_installed_components(coreDeviceThingName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/greengrass/v2/coreDevices/$(coreDeviceThingName)/installedComponents", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Retrieves the list of tags for an AWS IoT Greengrass resource.

# Arguments
- `resource_arn`: The ARN of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    resolve_component_candidates(component_candidates, platform)
    resolve_component_candidates(component_candidates, platform, params::Dict{String,<:Any})

Retrieves a list of components that meet the component, version, and platform requirements
of a deployment. AWS IoT Greengrass core devices call this operation when they receive a
deployment to identify the components to install. This operation identifies components that
meet all dependency requirements for a deployment. If the requirements conflict, then this
operation returns an error and the deployment fails. For example, this occurs if component
A requires version >2.0.0 and component B requires version <2.0.0 of a component
dependency. When you specify the component candidates to resolve, AWS IoT Greengrass
compares each component's digest from the core device with the component's digest in the
AWS Cloud. If the digests don't match, then AWS IoT Greengrass specifies to use the version
from the AWS Cloud.  To use this operation, you must use the data plane API endpoint and
authenticate with an AWS IoT device certificate. For more information, see AWS IoT
Greengrass endpoints and quotas.

# Arguments
- `component_candidates`: The list of components to resolve.
- `platform`: The platform to use to resolve compatible components.

"""
resolve_component_candidates(componentCandidates, platform; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/resolveComponentCandidates", Dict{String, Any}("componentCandidates"=>componentCandidates, "platform"=>platform); aws_config=aws_config)
resolve_component_candidates(componentCandidates, platform, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/greengrass/v2/resolveComponentCandidates", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("componentCandidates"=>componentCandidates, "platform"=>platform), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds tags to an AWS IoT Greengrass resource. If a tag already exists for the resource, this
operation updates the tag's value.

# Arguments
- `resource_arn`: The ARN of the resource to tag.
- `tags`: A list of key-value pairs that contain metadata for the resource. For more
  information, see Tag your resources in the AWS IoT Greengrass V2 Developer Guide.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes a tag from an AWS IoT Greengrass resource.

# Arguments
- `resource_arn`: The ARN of the resource to untag.
- `tag_keys`: A list of keys for tags to remove from the resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = greengrassv2("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)


text:
function bifurcate_test_functional()
    employees = Dict() # no employees :-)
    ok, shifts = bifurcate_unassigned_scheduler(test_delivery_unassigned_env)
    @test week_sum_coverage(test_delivery_unassigned_env) * TEST_THRESHOLD > week_hours_scheduled(shifts)
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
# # Advanced LRP usage
# One of the design goals of ExplainableAI.jl is to combine ease of use and
# extensibility for the purpose of research.
#
#
# This example will show you how to implement custom LRP rules and register custom layers
# and activation functions.
# For this purpose, we will quickly load our model from the previous section:
using ExplainableAI
using Flux
using MLDatasets
using ImageCore
using BSON

model = BSON.load("../model.bson", @__MODULE__)[:model]

index = 10
x, y = MNIST.testdata(Float32, index)
input = reshape(x, 28, 28, 1, :);

# ## Custom LRP composites
# Instead of creating an LRP-analyzer from a single rule (e.g. `LRP(model, GammaRule())`),
# we can also assign rules to each layer individually.
# For this purpose, we create an array of rules that matches the length of the Flux chain:
rules = [
    ZBoxRule(),
    GammaRule(),
    GammaRule(),
    EpsilonRule(),
    EpsilonRule(),
    EpsilonRule(),
    ZeroRule(),
    ZeroRule(),
]

analyzer = LRP(model, rules)
heatmap(input, analyzer)

# Since some Flux Chains contain other Flux Chains, ExplainableAI provides
# a utility function called [`flatten_model`](@ref).
#
#md # !!! warning "Flattening models"
#md #     Not all models can be flattened, e.g. those using
#md #     `Parallel` and `SkipConnection` layers.

# ## Custom LRP rules
# Let's define a rule that modifies the weights and biases of our layer on the forward pass.
# The rule has to be of type `AbstractLRPRule`.
struct MyGammaRule <: AbstractLRPRule end

# It is then possible to dispatch on the utility functions  [`modify_params`](@ref) and [`modify_denominator`](@ref)
# with the rule type `MyCustomLRPRule` to define custom rules without writing any boilerplate code.
# To extend internal functions, import them explicitly:
import ExplainableAI: modify_params

function modify_params(::MyGammaRule, W, b)
    ρW = W + 0.25 * relu.(W)
    ρb = b + 0.25 * relu.(b)
    return ρW, ρb
end

# We can directly use this rule to make an analyzer!
analyzer = LRP(model, MyGammaRule())
heatmap(input, analyzer)

# We just implemented our own version of the ``γ``-rule in 7 lines of code!
# The outputs match perfectly:
analyzer = LRP(model, GammaRule())
heatmap(input, analyzer)

# If the layer doesn't use weights and biases `W` and `b`, ExplainableAI provides a
# lower-level variant of [`modify_params`](@ref) called [`modify_layer`](@ref).
# This function is expected to take a layer and return a new, modified layer.

#md # !!! warning "Using modify_layer"
#md #
#md #     Use of the function `modify_layer` will overwrite functionality of `modify_params`
#md #     for the implemented combination of rule and layer types.
#md #     This is due to the fact that internally, `modify_params` is called by the default
#md #     implementation of `modify_layer`.
#md #
#md #     Therefore it is recommended to only extend `modify_layer` for a specific rule
#md #     and a specific layer type.

# ## Custom layers and activation functions
# ### Model checks for humans
# Good model checks and presets should allow novice users to apply XAI methods
# in a "plug & play" manner according to best practices.
#
# Let's say we define a layer that doubles its input:
struct MyDoublingLayer end
(::MyDoublingLayer)(x) = 2 * x

mylayer = MyDoublingLayer()
mylayer([1, 2, 3])

# Let's append this layer to our model:
model = Chain(model..., MyDoublingLayer())

# Creating an LRP analyzer, e.g. `LRPZero(model)`, will throw an `ArgumentError`
# and print a summary of the model check in the REPL:
# ```julia-repl
# ┌───┬───────────────────────┬─────────────────┬────────────┬────────────────┐
# │   │ Layer                 │ Layer supported │ Activation │ Act. supported │
# ├───┼───────────────────────┼─────────────────┼────────────┼────────────────┤
# │ 1 │ flatten               │            true │     —      │           true │
# │ 2 │ Dense(784, 100, relu) │            true │    relu    │           true │
# │ 3 │ Dense(100, 10)        │            true │  identity  │           true │
# │ 4 │ MyDoublingLayer()     │           false │     —      │           true │
# └───┴───────────────────────┴─────────────────┴────────────┴────────────────┘
#   Layers failed model check
#   ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
#
#   Found unknown layers MyDoublingLayer() that are not supported by ExplainableAI's LRP implementation yet.
#
#   If you think the missing layer should be supported by default, please submit an issue (https://github.com/adrhill/ExplainableAI.jl/issues).
#
#   These model checks can be skipped at your own risk by setting the LRP-analyzer keyword argument skip_checks=true.
#
#   [...]
# ```

# LRP should only be used on deep rectifier networks and ExplainableAI doesn't
# recognize `MyDoublingLayer` as a compatible layer.
# By default, it will therefore return an error and a model check summary
# instead of returning an incorrect explanation.
#
# However, if we know `MyDoublingLayer` is compatible with deep rectifier networks,
# we can register it to tell ExplainableAI that it is ok to use.
# This will be shown in the following section.

#md # !!! warning "Skipping model checks"
#md #
#md #     All model checks can be skipped at the user's own risk by setting the LRP-analyzer
#md #     keyword argument `skip_checks=true`.

# ### Registering custom layers
# The error in the model check will stop after registering our custom layer type
# `MyDoublingLayer` as "supported" by ExplainableAI.
#
# This is done using the function [`LRP_CONFIG.supports_layer`](@ref),
# which should be set to return `true` for the type `MyDoublingLayer`:
LRP_CONFIG.supports_layer(::MyDoublingLayer) = true

# Now we can create and run an analyzer without getting an error:
analyzer = LRPZero(model)
heatmap(input, analyzer)

#md # !!! note "Registering functions"
#md #
#md #     Flux's `Chains` can also contain functions, e.g. `flatten`.
#md #     This kind of layer can be registered as
#md #     ```julia
#md #     LRP_CONFIG.supports_layer(::typeof(mylayer)) = true
#md #     ```

# ### Registering activation functions
# The mechanism for registering custom activation functions is analogous to that of custom layers:
myrelu(x) = max.(0, x)
model = Chain(Flux.flatten, Dense(784, 100, myrelu), Dense(100, 10))

# Once again, creating an LRP analyzer for this model will throw an `ArgumentError`
# and display the following model check summary:
# ```julia-repl
# julia> analyzer = LRPZero(model3)
# ┌───┬─────────────────────────┬─────────────────┬────────────┬────────────────┐
# │   │ Layer                   │ Layer supported │ Activation │ Act. supported │
# ├───┼─────────────────────────┼─────────────────┼────────────┼────────────────┤
# │ 1 │ flatten                 │            true │     —      │           true │
# │ 2 │ Dense(784, 100, myrelu) │            true │   myrelu   │          false │
# │ 3 │ Dense(100, 10)          │            true │  identity  │           true │
# └───┴─────────────────────────┴─────────────────┴────────────┴────────────────┘
#   Activations failed model check
#   ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
#
#   Found layers with unknown or unsupported activation functions myrelu. LRP assumes that the model is a "deep rectifier network" that only contains ReLU-like activation functions.
#
#   If you think the missing activation function should be supported by default, please submit an issue (https://github.com/adrhill/ExplainableAI.jl/issues).
#
#   These model checks can be skipped at your own risk by setting the LRP-analyzer keyword argument skip_checks=true.
#
#   [...]
# ```

# Registation works by defining the function [`LRP_CONFIG.supports_activation`](@ref) as `true`:
LRP_CONFIG.supports_activation(::typeof(myrelu)) = true

# now the analyzer can be created without error:
analyzer = LRPZero(model)

# ## How it works internally
# Internally, ExplainableAI dispatches to low level functions
# ```julia
# function lrp!(rule, layer, Rₖ, aₖ, Rₖ₊₁)
#     Rₖ .= ...
# end
# ```
# These functions use the arguments `rule` and `layer` to dispatch
# `modify_params` and `modify_denominator` on the rule and layer type.
# They in-place modify a pre-allocated array of the input relevance `Rₖ`
# based on the input activation `aₖ` and output relevance `Rₖ₊₁`.
#
# Calling `analyze` then applies a foward-pass of the model, keeping track of
# the activations `aₖ` for each layer `k`.
# The relevance `Rₖ₊₁` is then set to the output neuron activation and the rules are applied
# in a backward-pass over the model layers and previous activations.

# ### Generic rule implementation using automatic differentiation
# The generic LRP rule–of which the ``0``-, ``\epsilon``- and ``\gamma``-rules are special cases–reads[^1][^2]:
# ```math
# R_{j}=\sum_{k} \frac{a_{j} \cdot \rho\left(w_{j k}\right)}{\epsilon+\sum_{0, j} a_{j} \cdot \rho\left(w_{j k}\right)} R_{k}
# ```
#
# where ``\rho`` is a function that modifies parameters – what we have so far called `modify_params`.
#
# The computation of this propagation rule can be decomposed into four steps:
# ```math
# \begin{array}{lr}
# \forall_{k}: z_{k}=\epsilon+\sum_{0, j} a_{j} \cdot \rho\left(w_{j k}\right) & \text { (forward pass) } \\
# \forall_{k}: s_{k}=R_{k} / z_{k} & \text { (element-wise division) } \\
# \forall_{j}: c_{j}=\sum_{k} \rho\left(w_{j k}\right) \cdot s_{k} & \text { (backward pass) } \\
# \forall_{j}: R_{j}=a_{j} c_{j} & \text { (element-wise product) }
# \end{array}
# ```
#
# For deep rectifier networks, the third step can also be written as the gradient computation
# ```math
# c_{j}=\left[\nabla\left(\sum_{k} z_{k}(\boldsymbol{a}) \cdot s_{k}\right)\right]_{j}
# ```
#
# and can be implemented via automatic differentiation (AD).
#
# This equation is implemented in ExplainableAI as the default method
# for all layer types that don't have a specialized implementation.
# We will refer to it as the "AD fallback".
#
# [^1]: G. Montavon et al., [Layer-Wise Relevance Propagation: An Overview](https://link.springer.com/chapter/10.1007/978-3-030-28954-6_10)
# [^2]: W. Samek et al., [Explaining Deep Neural Networks and Beyond: A Review of Methods and Applications](https://ieeexplore.ieee.org/document/9369420)

# ### AD fallback
# The default LRP fallback for unknown layers uses AD via [Zygote](https://github.com/FluxML/Zygote.jl).
# For `lrp!`, we end up with something that looks very similar to the previous four step computation:
# ```julia
# function lrp!(rule, layer, Rₖ, aₖ, Rₖ₊₁)
#     layerᵨ = modify_layer(rule, layer)
#     c = gradient(aₖ) do a
#             z = layerᵨ(a)
#             s = Zygote.@ignore Rₖ₊₁ ./ modify_denominator(rule, z)
#             z ⋅ s
# 	      end |> only
#     Rₖ .= aₖ .* c
# end
# ```
#
# You can see how `modify_layer` and `modify_denominator` dispatch on the rule and layer type.
# This is how we implemented our own `MyGammaRule`.
# Unknown layers that are registered in the `LRP_CONFIG` use this exact function.

# ### Specialized implementations
# We can also implement specialized versions of `lrp!` based on the type of `layer`,
# e.g. reshaping layers.
#
# Reshaping layers don't affect attributions. We can therefore avoid the computational
# overhead of AD by writing a specialized implementation that simply reshapes back:
# ```julia
# function lrp!(::AbstractLRPRule, ::ReshapingLayer, Rₖ, aₖ, Rₖ₊₁)
#     Rₖ .= reshape(Rₖ₊₁, size(aₖ))
# end
# ```
#
# Since the rule type didn't matter in this case, we didn't specify it.
#
# We can even implement the generic rule as a specialized implementation for `Dense` layers:
# ```julia
# function lrp!(rule::AbstractLRPRule, layer::Dense, Rₖ, aₖ, Rₖ₊₁)
#     ρW, ρb = modify_params(rule, get_params(layer)...)
#     ãₖ₊₁ = modify_denominator(rule, ρW * aₖ + ρb)
#     @tullio Rₖ[j] = aₖ[j] * ρW[k, j] / ãₖ₊₁[k] * Rₖ₊₁[k] # Tullio ≈ fast einsum
# end
# ```
#
# For maximum low-level control beyond `modify_layer`, `modify_params` and `modify_denominator`,
# you can also implement your own `lrp!` function and dispatch
# on individual rule types `MyRule` and layer types `MyLayer`:
# ```julia
# function lrp!(rule::MyRule, layer::MyLayer, Rₖ, aₖ, Rₖ₊₁)
#     Rₖ .= ...
# end
# ```


text:
using WriteVTK
using Test

function main()
   # Suppose that the mesh is made of 5 points:
  cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
           MeshCell(VTKCellTypes.VTK_QUAD, [2, 4, 3, 5])]
  x=[.95,.16,.07,.21,.60]
  y=[.32,.55,.87,.12,.85]
  vtufile = "simulation/simulation_1.vtu"
  rm(vtufile,force=true)
  @time pvtk = pvtk_grid("simulation", x, y, cells; part=1,nparts=1) # 2D
  pvtk["Pressure"] = x
  pvtk["Processor"] = [1,2]
  pvtk["Temperature",VTKPointData()] = y
  pvtk["Id",VTKCellData()] = [2,1]
  @time outfiles = vtk_save(pvtk)
  @test isfile(vtufile)
  println("Saved:  ", join(outfiles, "  "))
  outfiles
end
main()


text:
module Schedulers

include("../../docstrings.jl")

import ..Rule
import ..EGraph
using ..Util

include("./abstractscheduler.jl")
include("./backoffscheduler.jl")
include("./simplescheduler.jl")
include("./scoredscheduler.jl")

export AbstractScheduler
export SimpleScheduler
export BackoffScheduler
export ScoredScheduler

export cansaturate
export shouldskip
export readstep!
export writestep!

end


text:
using Nettle
const hmacstate = HMACState[]
function hmac(s1,s2,s3,s4)
    if isempty(hmacstate)
        return ""
    else
        update!(hmacstate[1], s1)
        update!(hmacstate[1], s2)
        update!(hmacstate[1], s3)
        update!(hmacstate[1], s4)
        return hexdigest!(hmacstate[1])
    end
end


text:
function curandGetVersion()
    ver = Ref{Cint}()
    curandGetVersion(ver)
    return ver[]
end

function curandGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end


text:
abstract type abstract_reduct_factors end
abstract type abstract_pars end

mutable struct Class_kkt_solver_options <: abstract_pars
    ItRefine_Num::Int64
    ItRefine_BigFloat::Bool
    saddle_err_tol::Float64
    kkt_solver_type::Symbol
    kkt_system_rescale::Symbol
    linear_solver_type::Symbol
    linear_solver_safe_mode::Bool
    linear_solver_recycle::Bool
    ma97_u::Float64
    stable_reduct_factors::abstract_reduct_factors
    aggressive_reduct_factors::abstract_reduct_factors

    function Class_kkt_solver_options()
      this = new()

      this.ItRefine_Num = 3
      this.ItRefine_BigFloat = false
      this.saddle_err_tol = Inf

      this.kkt_system_rescale = :none
      this.ma97_u = 1e-8
      #this.kkt_system_rescale = :u_only
      #this.kkt_system_rescale = :u_and_x
      #this.kkt_system_rescale = :u_and_x
      if true
        this.kkt_solver_type = :schur #_direct
        this.linear_solver_type = :julia
        #this.linear_solver_type = :mumps
      else
        this.kkt_solver_type = :symmetric
        this.linear_solver_type = :mumps
      end
      this.linear_solver_safe_mode = false #true
      this.linear_solver_recycle = false

      # Don't change these parameters except for experimentation
      this.stable_reduct_factors = Reduct_stable()
      this.aggressive_reduct_factors = Reduct_affine()

      return this
    end
end

mutable struct Class_line_search_parameters <: abstract_pars
  kkt_reduction_factor::Float64
  kkt_include_comp::Bool
  filter_type::Symbol

  predict_reduction_factor::Float64

  fraction_to_boundary_predict::Float64
  fraction_to_boundary::Float64
  fraction_to_boundary_predict_exp::Float64
  fraction_to_boundary_linear::Float64

  backtracking_factor::Float64
  num_backtracks::Int64

  agg_gamma::Symbol

  move_primal_seperate_to_dual::Bool
  dual_ls::Int64

  # to add
  comp_feas::Float64
  comp_feas_agg::Float64
  min_step_size_stable::Float64
  min_step_size_agg_ratio::Float64

  function Class_line_search_parameters()
    this = new()
    # LINE SEARCH
    this.kkt_reduction_factor = 0.2
    this.kkt_include_comp = true
    this.filter_type = :test2

    this.predict_reduction_factor = 0.1

    this.fraction_to_boundary = 0.1 #0.1
    this.fraction_to_boundary_predict = 0.2 #0.25
    this.fraction_to_boundary_linear = 0.1
    this.fraction_to_boundary_predict_exp = 0.5

    this.backtracking_factor = 0.5
    this.num_backtracks = 60;

    this.agg_gamma = :mehrotra_stb
    #this.agg_gamma = :affine

    this.move_primal_seperate_to_dual = true
    this.dual_ls = 1

    this.comp_feas = 1/100.0
    this.comp_feas_agg = 1/50.0
    this.min_step_size_stable = 0.5^5.0
    this.min_step_size_agg_ratio = 1e-4

    return this
  end
end

mutable struct Class_IPM_parameters
  # fill in
end

mutable struct Class_termination_parameters <: abstract_pars
    max_it::Int64
    max_time::Float64
    tol_opt::Float64
    tol_unbounded::Float64
    grad_max::Float64
    tol_inf_1::Float64
    tol_inf_2::Float64
    dual_scale_threshold::Float64
    dual_scale_mode::Symbol

    function Class_termination_parameters()
        this = new()

        this.max_it = 3000
        this.max_time = 60.0^2
        this.tol_opt = 1e-6
        this.tol_unbounded = 1e-12
        this.grad_max = Inf
        this.tol_inf_1 = 1e-3
        this.tol_inf_2 = 1e-6
        this.dual_scale_threshold = 100.0;
        this.dual_scale_mode = :max_dual

        return this
    end
end

mutable struct Class_delta_parameters <: abstract_pars
    max::Float64
    #max_it::Int64
    start::Float64
    dec::Float64
    inc::Float64
    zero::Float64 #get_mu(iter) / ((1e2 + LinearAlgebra.norm(get_x(iter),Inf)) * 1e2)
    min::Float64

    function Class_delta_parameters()
        this = new()
        this.max = 10.0^(50.0)
        #this.max_it = 200
        this.start = 1e-6
        this.zero = 0.0
        this.min = 1e-12
        this.inc = 8.0
        this.dec = 1.0 / pi

        return this
    end
end

mutable struct Class_init_parameters <: abstract_pars
    mu_scale::Float64
    mehotra_scaling::Bool
    init_style::Symbol
    start_satisfying_bounds::Bool
    dual_threshold::Float64
    linear_scale::Float64
    nl_ineq_scale::Float64
    nl_eq_scale::Float64
    dual_max::Float64
    dual_min::Float64

    function Class_init_parameters()
        this = new()
        #mode = :experimentation
        mode = :standard
        if mode == :standard
          this.mu_scale = 1.0
          this.mehotra_scaling = true
          this.init_style = :gertz # #
          this.dual_threshold = 1.0
          this.start_satisfying_bounds = true
          this.linear_scale = 1.0
          this.nl_ineq_scale = 1.0
          this.nl_eq_scale = 1.0
          this.dual_max = 1e3
          this.dual_min = 1e-2
        else
          this.mu_scale = 1.0
          this.mehotra_scaling = false
          this.init_style = :mehrotra
          this.dual_threshold = 1.0
          this.start_satisfying_bounds = true
          this.linear_scale = 1e-2
          this.nl_ineq_scale = 1e-1
          this.nl_eq_scale = 1.0
        end

        return this
    end
end

mutable struct Class_testing <: abstract_pars
    response_to_failure::Symbol
    function Class_testing()
        this = new()
        #this.response_to_failure = :default
        this.response_to_failure = :lag_delta_inc

        return this
    end
end

mutable struct Class_parameters <: abstract_pars
    term::Class_termination_parameters
    init::Class_init_parameters
    delta::Class_delta_parameters
    ls::Class_line_search_parameters
    kkt::Class_kkt_solver_options
    test::Class_testing

    # debugging
    output_level::Int64
    debug_mode::Int64
    throw_error_nans::Bool

    # IPM GENERAL
    aggressive_dual_threshold::Float64
    max_it_corrections::Int64

    superlinear_theory_mode::Bool
    agg_protection_factor::Float64

    kkt_include_comp::Bool
    primal_bounds_dual_feas::Bool
    a_norm_penalty::Float64

    eps_mach::Float64

    function Class_parameters()
        this = new()

        ############################
        ### groups of parameters ###
        ############################
        this.delta = Class_delta_parameters()
        this.term = Class_termination_parameters()
        this.init = Class_init_parameters()
        this.ls = Class_line_search_parameters()
        this.kkt = Class_kkt_solver_options()
        this.test = Class_testing()

        ########################
        ## general parameters ##
        ########################

        # debugging
        this.output_level = 2
        this.debug_mode = 0
        this.throw_error_nans = false

        # switching condition
        this.aggressive_dual_threshold = 1.0 # kappa_{1}
        this.primal_bounds_dual_feas = false

        # algorithm parameters
        this.max_it_corrections = 2
        this.superlinear_theory_mode = false
        this.agg_protection_factor = 0.9

        # penalty parameter size, changes log barrier to stop iterates growing too big when feasible region is unbounded
        # e.g. min x s.t. x, y >= 0, without any changes y -> \infty
        this.a_norm_penalty = 1e-4

        # machine precision
        this.eps_mach = 1e-16

        return this
    end
end

#::IOStream
function write_pars(stream, par::Class_parameters)
    write(stream, "PAR \t\t\t\t\t\t VALUE \n")

    for fieldname in fieldnames(par)
        fieldval = getfield(par, fieldname)
        if isa(fieldval,abstract_pars)
          write(stream, "$(pd(fieldname,40)) \n")
          for subfieldname in fieldnames(fieldval)
            subfieldval = getfield(fieldval, subfieldname)
            write(stream, "\t $(pd(subfieldname,40)) \t $(subfieldval) \n")
          end
        else
          write(stream, "$(pd(fieldname,40)) \t $(fieldval) \n")
        end
    end
end


text:
# Invertible conditional HINT network from Kruse et. al (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: February 2020

export NetworkMultiScaleConditionalHINT

"""
    CH = NetworkMultiScaleConditionalHINT(nx, ny, n_in, batchsize, n_hidden,  L, K; split_scales=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create a conditional HINT network for data-driven generative modeling based
 on the change of variables formula.

 *Input*:

 - `nx`, `ny`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensors `X` and `Y`

 - `n_hidden`: number of hidden units in residual blocks

 - `L`: number of scales (outer loop)

 - `K`: number of flow steps per scale (inner loop)

 - `split_scales`: if true, split output in half along channel dimension after each scale. Feed one half through the next layers,
    while saving the remaining channels for the output.

 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)

 - `p1`, `p2`: respective padding sizes for residual block layers

 - `s1`, `s2`: respective strides for residual block layers

 *Output*:

 - `CH`: conditional HINT network

 *Usage:*

 - Forward mode: `Zx, Zy, logdet = CH.forward(X, Y)`

 - Inverse mode: `X, Y = CH.inverse(Zx, Zy)`

 - Backward mode: `ΔX, X = CH.backward(ΔZx, ΔZy, Zx, Zy)`

 *Trainable parameters:*

 - None in `CH` itself

 - Trainable parameters in activation normalizations `CH.AN_X[i]` and `CH.AN_Y[i]`,
 and in coupling layers `CH.CL[i]`, where `i` ranges from `1` to `depth`.

 See also: [`ActNorm`](@ref), [`ConditionalLayerHINT!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct NetworkMultiScaleConditionalHINT <: InvertibleNetwork
    AN_X::AbstractArray{ActNorm, 2}
    AN_Y::AbstractArray{ActNorm, 2}
    CL::AbstractArray{ConditionalLayerHINT, 2}
    XY_dims::Union{Array{Tuple, 1}, Nothing}
    L::Int64
    K::Int64
    split_scales::Bool
end

@Flux.functor NetworkMultiScaleConditionalHINT

# Constructor
function NetworkMultiScaleConditionalHINT(nx::Int64, ny::Int64, n_in::Int64, batchsize::Int64,
                                          n_hidden::Int64, L::Int64, K::Int64;
                                          split_scales=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

    AN_X = Array{ActNorm}(undef, L, K)
    AN_Y = Array{ActNorm}(undef, L, K)
    CL = Array{ConditionalLayerHINT}(undef, L, K)
    if split_scales
        XY_dims = Array{Tuple}(undef, L-1)
        channel_factor = 2
    else
        XY_dims = nothing
        channel_factor = 4
    end

    # Create layers
    for i=1:L
        for j=1:K
            AN_X[i, j] = ActNorm(n_in*4; logdet=true)
            AN_Y[i, j] = ActNorm(n_in*4; logdet=true)
            CL[i, j] = ConditionalLayerHINT(Int(nx/2^i), Int(ny/2^i), n_in*4, n_hidden, batchsize;
                                            permute=true, k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
        end
        n_in *= channel_factor
    end

    return NetworkMultiScaleConditionalHINT(AN_X, AN_Y, CL, XY_dims, L, K, split_scales)
end

# Concatenate states Zi and final output
function cat_states(XY_save::AbstractArray{Array, 2}, X::AbstractArray{Float32, 4}, Y::AbstractArray{Float32, 4})
    X_full = []
    Y_full = []
    for j=1:size(XY_save, 1)
        X_full = cat(X_full, vec(XY_save[j, 1]); dims=1)
        Y_full = cat(Y_full, vec(XY_save[j, 2]); dims=1)
    end
    X_full = cat(X_full, vec(X); dims=1)
    Y_full = cat(Y_full, vec(Y); dims=1)
    return Float32.(X_full), Float32.(Y_full)  # convert to Array{Float32, 1}
end

function cat_states(Y_save::AbstractArray{Array, 1}, Y::AbstractArray{Float32, 4})
    Y_full = []
    for j=1:size(Y_save, 1)
        Y_full = cat(Y_full, vec(Y_save[j, 1]); dims=1)
    end
    Y_full = cat(Y_full, vec(Y); dims=1)
    return Float32.(Y_full)
end

# Split 1D vector in latent space back to states Zi
function split_states(XY_dims::AbstractArray{Tuple, 1}, X_full::AbstractArray{Float32, 1}, Y_full::AbstractArray{Float32, 1})
    L = length(XY_dims) + 1
    XY_save = Array{Array}(undef, L-1, 2)
    count = 1
    for j=1:L-1
        XY_save[j, 1] = reshape(X_full[count: count + prod(XY_dims[j])-1], XY_dims[j])
        XY_save[j, 2] = reshape(Y_full[count: count + prod(XY_dims[j])-1], XY_dims[j])
        count += prod(XY_dims[j])
    end
    X = reshape(X_full[count: count + prod(XY_dims[end])-1], Int.(XY_dims[end].*(.5, .5, 4, 1)))
    Y = reshape(Y_full[count: count + prod(XY_dims[end])-1], Int.(XY_dims[end].*(.5, .5, 4, 1)))
    return XY_save, X, Y
end

# Split 1D vector in latent space back to states Zi
function split_states(XY_dims::AbstractArray{Tuple, 1}, Y_full::AbstractArray{Float32, 1})
    L = length(XY_dims) + 1
    Y_save = Array{Array}(undef, L-1)
    count = 1
    for j=1:L-1
        Y_save[j] = reshape(Y_full[count: count + prod(XY_dims[j])-1], XY_dims[j])
        count += prod(XY_dims[j])
    end
    Y = reshape(Y_full[count: count + prod(XY_dims[end])-1], Int.(XY_dims[end].*(.5, .5, 4, 1)))
    return Y_save, Y
end

# Forward pass and compute logdet
function forward(X, Y, CH::NetworkMultiScaleConditionalHINT)
    CH.split_scales && (XY_save = Array{Array}(undef, CH.L-1, 2))
    logdet = 0f0
    for i=1:CH.L
        X = wavelet_squeeze(X)
        Y = wavelet_squeeze(Y)
        for j=1:CH.K
            X_, logdet1 = CH.AN_X[i, j].forward(X)
            Y_, logdet2 = CH.AN_Y[i, j].forward(Y)
            X, Y, logdet3 = CH.CL[i, j].forward(X_, Y_)
            logdet += (logdet1 + logdet2 + logdet3)
        end
        if CH.split_scales && i < CH.L    # don't split after last iteration
            X, Zx = tensor_split(X)
            Y, Zy = tensor_split(Y)
            XY_save[i, :] = [Zx, Zy]
            CH.XY_dims[i] = size(Zx)
        end
    end
    CH.split_scales && ((X, Y) = cat_states(XY_save, X, Y))
    return X, Y, logdet
end

# Inverse pass and compute gradients
function inverse(Zx, Zy, CH::NetworkMultiScaleConditionalHINT)
    CH.split_scales && ((XY_save, Zx, Zy) = split_states(CH.XY_dims, Zx, Zy))
    for i=CH.L:-1:1
        if CH.split_scales && i < CH.L
            Zx = tensor_cat(Zx, XY_save[i, 1])
            Zy = tensor_cat(Zy, XY_save[i, 2])
        end
        for j=CH.K:-1:1
            Zx_, Zy_ = CH.CL[i, j].inverse(Zx, Zy)
            Zy = CH.AN_Y[i, j].inverse(Zy_; logdet=false)
            Zx = CH.AN_X[i, j].inverse(Zx_; logdet=false)
        end
        Zx = wavelet_unsqueeze(Zx)
        Zy = wavelet_unsqueeze(Zy)
    end
    return Zx, Zy
end

# Backward pass and compute gradients
function backward(ΔZx, ΔZy, Zx, Zy, CH::NetworkMultiScaleConditionalHINT)

    # Split data and gradients
    if CH.split_scales
        ΔXY_save, ΔZx, ΔZy = split_states(CH.XY_dims, ΔZx, ΔZy)
        XY_save, Zx, Zy = split_states(CH.XY_dims, Zx, Zy)
    end

    for i=CH.L:-1:1
        if CH.split_scales && i < CH.L
            ΔZx = tensor_cat(ΔZx, ΔXY_save[i, 1])
            ΔZy = tensor_cat(ΔZy, ΔXY_save[i, 2])
            Zx = tensor_cat(Zx, XY_save[i, 1])
            Zy = tensor_cat(Zy, XY_save[i, 2])
        end
        for j=CH.K:-1:1
            ΔZx_, ΔZy_, Zx_, Zy_ = CH.CL[i, j].backward(ΔZx, ΔZy, Zx, Zy)
            ΔZx, Zx = CH.AN_X[i, j].backward(ΔZx_, Zx_)
            ΔZy, Zy = CH.AN_Y[i, j].backward(ΔZy_, Zy_)
        end
        ΔZx = wavelet_unsqueeze(ΔZx)
        ΔZy = wavelet_unsqueeze(ΔZy)
        Zx = wavelet_unsqueeze(Zx)
        Zy = wavelet_unsqueeze(Zy)
    end
    return ΔZx, ΔZy, Zx, Zy
end

# Forward pass and compute logdet
function forward_Y(Y, CH::NetworkMultiScaleConditionalHINT)
    CH.split_scales && (Y_save = Array{Array}(undef, CH.L-1))
    for i=1:CH.L
        Y = wavelet_squeeze(Y)
        for j=1:CH.K
            Y_, logdet2 = CH.AN_Y[i, j].forward(Y)
            Y = CH.CL[i, j].forward_Y(Y_)
        end
        if CH.split_scales && i < CH.L    # don't split after last iteration
            Y, Zy = tensor_split(Y)
            Y_save[i] = Zy
            CH.XY_dims[i] = size(Zy)
        end
    end
    CH.split_scales && (Y = cat_states(Y_save, Y))
    return Y

end

# Inverse pass and compute gradients
function inverse_Y(Zy, CH::NetworkMultiScaleConditionalHINT)
    CH.split_scales && ((Y_save, Zy) = split_states(CH.XY_dims, Zy))
    for i=CH.L:-1:1
        if CH.split_scales && i < CH.L
            Zy = tensor_cat(Zy, Y_save[i])
        end
        for j=CH.K:-1:1
            Zy_ = CH.CL[i, j].inverse_Y(Zy)
            Zy = CH.AN_Y[i, j].inverse(Zy_; logdet=false)
        end
        Zy = wavelet_unsqueeze(Zy)
    end
    return Zy
end

# Clear gradients
function clear_grad!(CH::NetworkMultiScaleConditionalHINT)
    depth = length(CH.CL)
    for j=1:depth
        clear_grad!(CH.AN_X[j])
        clear_grad!(CH.AN_Y[j])
        clear_grad!(CH.CL[j])
    end
end

# Get parameters
function get_params(CH::NetworkMultiScaleConditionalHINT)
    depth = length(CH.CL)
    p = []
    for j=1:depth
        p = cat(p, get_params(CH.AN_X[j]); dims=1)
        p = cat(p, get_params(CH.AN_Y[j]); dims=1)
        p = cat(p, get_params(CH.CL[j]); dims=1)
    end
    return p
end

# Put parameters
function put_params!(CH::NetworkMultiScaleConditionalHINT, Params::Array{Any,1})
    depth = length(CH.CL)
    idx_s = 1
    counter = 0

    for j = 1:depth
        idx_s += counter
        counter = 2
        put_params!(CH.AN_X[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 2
        put_params!(CH.AN_Y[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 0
        nlayers_x = length(CH.CL[j].CL_X.CL)
        isnothing(CH.CL[j].CL_X.C) ? counter += 5*nlayers_x : counter += 5*nlayers_x + 3
        nlayers_y = length(CH.CL[j].CL_Y.CL)
        isnothing(CH.CL[j].CL_Y.C) ? counter += 5*nlayers_y : counter += 5*nlayers_y + 3
        counter += 5
        ~isnothing(CH.CL[j].C_X) && (counter += 3)
        ~isnothing(CH.CL[j].C_Y) && (counter += 3)

        put_params!(CH.CL[j], Params[idx_s:idx_s+counter-1])
    end
end


text:
using SparseArrays

@inline function speye(dtype, n)
    return spdiagm(0=>ones(dtype, n))
end

# Assuming x-direction is periodic
struct BasicMatrixOperators

    U_pts
    V_pts
    W_pts
    T_pts
    UV_pts

    # Nomenclature:
    #
    # [new-grid][direction][old-grid]
    #
    # U_W_T : sending variable westward from T grid to U grid

    U_I_U
    V_I_V
    T_I_T
    UV_I_UV

    U_W_T
    U_E_T
    U_W_U
    U_E_U
    V_W_V
    V_E_V
    V_W_UV
    V_E_UV

    V_S_T
    V_N_T
    U_S_U
    U_N_U
    V_S_V
    V_N_V
    U_S_UV
    U_N_UV

    T_S_V
    T_N_V
    T_W_U
    T_E_U

    UV_E_V
    UV_W_V
    UV_N_U
    UV_S_U

    U_SW_V
    U_SE_V
    U_NW_V
    U_NE_V
    V_SW_U
    V_SE_U
    V_NW_U
    V_NE_U

    T_SW_UV
    T_SE_UV
    T_NW_UV
    T_NE_UV

    UV_SW_T
    UV_SE_T
    UV_NW_T
    UV_NE_T

    T_N_T
    T_S_T
    T_E_T
    T_W_T

    T_UP_T
    T_DN_T
    U_UP_U
    U_DN_U
    V_UP_V
    V_DN_V

    T_UP_W
    T_DN_W

    W_UP_T
    W_DN_T

    function BasicMatrixOperators(;
        Nx             :: Int64,
        Ny             :: Int64,
        Nz             :: Int64,
    )

       #println("Making MatrixOperators")
        # Making operator
        U_dim = (Nz, Nx, Ny)
        V_dim = (Nz, Nx, Ny+1)
        W_dim = (Nz+1, Nx, Ny)
        T_dim = (Nz, Nx, Ny)
        UV_dim = (Nz, Nx, Ny+1)

        U_pts = reduce(*, U_dim)
        V_pts = reduce(*, V_dim)
        W_pts = reduce(*, W_dim)
        T_pts = reduce(*, T_dim)
        UV_pts = reduce(*, UV_dim)

        U_I_U = speye(Float64, U_pts)
        V_I_V = speye(Float64, V_pts)
        W_I_W = speye(Float64, W_pts)
        T_I_T = speye(Float64, T_pts)
        UV_I_UV = speye(Float64, UV_pts)

        U_I_U_expand = vcat(U_I_U, zeros(Float64, 1, U_pts))
        V_I_V_expand = vcat(V_I_V, zeros(Float64, 1, V_pts))
        W_I_W_expand = vcat(W_I_W, zeros(Float64, 1, W_pts))
        T_I_T_expand = vcat(T_I_T, zeros(Float64, 1, T_pts))
        UV_I_UV_expand = vcat(UV_I_UV, zeros(Float64, 1, UV_pts))

        num_U = zeros(Int64, U_dim...)
        num_V = zeros(Int64, V_dim...)
        num_W = zeros(Int64, W_dim...)
        num_T = zeros(Int64, T_dim...)
        num_UV = zeros(Int64, UV_dim...)

        num_U[:] = 1:length(num_U)
        num_V[:] = 1:length(num_V)
        num_W[:] = 1:length(num_W)
        num_T[:] = 1:length(num_T)
        num_UV[:] = 1:length(num_UV)

        U = num_U * 0
        V = num_V * 0
        W = num_W * 0
        T = num_T * 0
        F = num_T * 0

        #smb = SparseMatrixBuilder(Nx*(Ny+1)*(Nz+1)*4)
        function build!(id_mtx, idx; wipe=:none)
           #println("Build!")
            local result
            rows = size(id_mtx)[1]
            if wipe == :n
                idx[:, :, end] .= rows
            elseif wipe == :s
                idx[:, :,   1] .= rows
            elseif wipe == :t
                idx[1, :, :] .= rows
            elseif wipe == :b
                idx[end, :, :] .= rows
            elseif wipe != :none
                throw(ErrorException("Wrong keyword"))
            end

            # using transpose speeds up by 100 times
            tp = transpose(id_mtx) |> sparse
            result = transpose(tp[:, view(idx, :)]) |> sparse
            #result = id_mtx[view(idx, :), :]
            #dropzeros!(result)

            idx .= 0 # clean so that debug is easir when some girds are not assigned
            return result
        end

       #println("Making shifting operators")
        # east, west passing mtx
        U[:, :, :] = num_T;                             U_W_T = build!(T_I_T, U)
        U[:, :, :] = circshift(num_T, (0, 1, 0));       U_E_T = build!(T_I_T, U)

        U[:] = circshift(num_U, (0, -1, 0));            U_W_U = build!(U_I_U, U)
        U[:] = circshift(num_U, (0,  1, 0));            U_E_U = build!(U_I_U, U)

        V[:] = circshift(num_V, (0, -1, 0));            V_W_V = build!(V_I_V, V)
        V[:] = circshift(num_V, (0,  1, 0));            V_E_V = build!(V_I_V, V)

        V[:, :, :] = circshift(num_UV, (0, -1, 0));      V_W_UV = build!(UV_I_UV, V)
        V[:, :, :] = num_UV;                             V_E_UV = build!(UV_I_UV, V)

        # north, south passing mtx
        V[:, :, 1:Ny, ] = num_T;                        V_S_T = build!(T_I_T_expand, V; wipe=:n)
        V[:, :, 2:Ny+1] = num_T;                        V_N_T = build!(T_I_T_expand, V; wipe=:s)

        U[:, :, 1:Ny-1] = view(num_U, :, :, 2:Ny, );    U_S_U = build!(U_I_U_expand, U; wipe=:n)
        U[:, :, 2:Ny  ] = view(num_U, :, :, 1:Ny-1);    U_N_U = build!(U_I_U_expand, U; wipe=:s)

        V[:, :, 1:Ny, ] = view(num_V, :, :, 2:Ny+1);    V_S_V = build!(V_I_V_expand, V; wipe=:n)
        V[:, :, 2:Ny+1] = view(num_V, :, :, 1:Ny  );    V_N_V = build!(V_I_V_expand, V; wipe=:s)

        U[:, :, :]      = view(num_UV, :, :, 2:Ny+1);    U_S_UV = build!(UV_I_UV, U)
        U[:, :, :]      = view(num_UV, :, :, 1:Ny  );    U_N_UV = build!(UV_I_UV, U)

        # inverse directions
        T_S_V = V_N_T' |> sparse
        T_N_V = V_S_T' |> sparse
        T_W_U = U_E_T' |> sparse
        T_E_U = U_W_T' |> sparse

        UV_E_V = V_W_UV' |> sparse
        UV_W_V = V_E_UV' |> sparse
        UV_N_U = U_S_UV' |> sparse
        UV_S_U = U_N_UV' |> sparse

        # diagonal passing mtx
        U_SW_V = U_W_T * T_S_V
        U_SE_V = U_E_T * T_S_V
        U_NW_V = U_W_T * T_N_V
        U_NE_V = U_E_T * T_N_V

        V_SW_U = V_S_T * T_W_U
        V_SE_U = V_S_T * T_E_U
        V_NW_U = V_N_T * T_W_U
        V_NE_U = V_N_T * T_E_U

        T_SW_UV = T_S_V * V_W_UV
        T_SE_UV = T_S_V * V_E_UV
        T_NW_UV = T_N_V * V_W_UV
        T_NE_UV = T_N_V * V_E_UV

        UV_SW_T = T_NE_UV' |> sparse
        UV_SE_T = T_NW_UV' |> sparse
        UV_NW_T = T_SE_UV' |> sparse
        UV_NE_T = T_SW_UV' |> sparse

        # T to T operators
        T_N_T = T_N_V * V_N_T
        T_S_T = T_S_V * V_S_T
        T_E_T = T_E_U * U_E_T
        T_W_T = T_W_U * U_W_T

        # upward, downward passing mtx
        T[1:Nz-1, :, :] = view(num_T, 2:Nz,   :, :);    T_UP_T = build!(T_I_T_expand, T; wipe=:b)
        T[2:Nz,   :, :] = view(num_T, 1:Nz-1, :, :);    T_DN_T = build!(T_I_T_expand, T; wipe=:t)

        U[1:Nz-1, :, :] = view(num_U, 2:Nz  , :, :);    U_UP_U = build!(U_I_U_expand, U; wipe=:b)
        U[2:Nz,   :, :] = view(num_U, 1:Nz-1, :, :);    U_DN_U = build!(U_I_U_expand, U; wipe=:t)

        V[1:Nz-1, :, :] = view(num_V, 2:Nz,   :, :);    V_UP_V = build!(V_I_V_expand, V; wipe=:b)
        V[2:Nz  , :, :] = view(num_V, 1:Nz-1, :, :);    V_DN_V = build!(V_I_V_expand, V; wipe=:t)

        T[:, :, :] = view(num_W, 2:Nz+1, :, :);         T_UP_W = build!(W_I_W_expand, T)
        T[:, :, :] = view(num_W, 1:Nz  , :, :);         T_DN_W = build!(W_I_W_expand, T)

        # inverse directions
        W_DN_T = T_UP_W' |> sparse
        W_UP_T = T_DN_W' |> sparse

        return new(
            U_pts, V_pts, W_pts, T_pts, UV_pts,
            U_I_U, V_I_V, T_I_T, UV_I_UV,

            U_W_T, U_E_T,
            U_W_U, U_E_U,
            V_W_V, V_E_V,
            V_W_UV, V_E_UV,

            V_S_T, V_N_T,
            U_S_U, U_N_U,
            V_S_V, V_N_V,
            U_S_UV, U_N_UV,

            T_S_V, T_N_V,
            T_W_U, T_E_U,

            UV_E_V, UV_W_V,
            UV_N_U, UV_S_U,

            U_SW_V, U_SE_V,
            U_NW_V, U_NE_V,
            V_SW_U, V_SE_U,
            V_NW_U, V_NE_U,

            T_SW_UV, T_SE_UV,
            T_NW_UV, T_NE_UV,

            UV_SW_T, UV_SE_T,
            UV_NW_T, UV_NE_T,

            T_N_T, T_S_T,
            T_E_T, T_W_T,

            T_UP_T, T_DN_T,
            U_UP_U, U_DN_U,
            V_UP_V, V_DN_V,
            T_UP_W, T_DN_W,
            W_UP_T, W_DN_T,
        )

    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using LinearAlgebra

#Givens from Problem Statement
#Elastic Modulus(psi),Area(in^2),length(in),Force(lb)
E = 1.9*10^6
A = 8
l = 3*12
l2 = 3*sqrt(2)*12
F = 500

# Torsional Stiffness Function
k(A,E,l) = (A*E)/l

# k_local(Elastic Modulus,Area,length,theta,number of nodes,
# row/column one,row/column two, ==1 left fixed/ ==2 right fixed)
function k_local(E,A,l,theta,i,n1,n2,fixed=0)
  #k_local = zeros(4,4)
  p1 = (n1*2)-1
  p2 = n1*2
  p3 = (n2*2)-1
  p4 = n2*2

  k_global = zeros(2*i,2*i)
  kin = k(A,E,l)
  s = (sind(theta))^2
	c = (cosd(theta))^2
	sc = sind(theta)*cosd(theta)
	k_quad = [c sc;
            sc s]
  k_global[p1:p2,p1:p2] = k*k_quad
  k_global[p1:p2,p3:p4] = -k*k_quad
  k_global[p3:p4,p1:p2] = -k*k_quad
  k_global[p3:p4,p3:p4] = k*k_quad
  # Need to remove rows and columns; check: push!, pop!
  if fixed == 1
      adj = zeros(1,i)
      #adj[n1] = 1
      k_global[p1,:] = adj
      k_global[:,p1] = adj
      k_global[p2,:] = adj
      k_global[:,p2] = adj
    elseif fixed == 2
      adj = zeros(1,i)
      #adj[n2] = 1
      k_global[p3,:] = adj
      k_global[:,p3] = adj
      k_global[p4,:] = adj
      k_global[:,p4] = adj
    end
    return k_global
end

# Unadjusted Global
k_1 = k_local(E,A,l,0,5,1,2,1)
k_2 = k_local(E,A,l,135,5,2,3)
k_3 = k_local(E,A,l,0,5,3,4,1)
k_4 = k_local(E,A,l,90,5,2,4)
k_5 = k_local(E,A,l,45,5,2,5)
k_6 = k_local(E,A,l,0,5,4,5)
kG = k_1+k_2+k_3+k_4+k_5+k_6

#Adjusted Global
kG_adj = kG

# Load Vector
p = Float64[0 0 0 -F 0 -F]
P = transpose(p)

# Angular Displacement Vector
u = kG_adj\P
println("Displacement Vector (in)")
display(u)

# Moment Reaction Vector
P_abs = kG*u_adj
println("Load Reaction Vector (lb):")
display(P_abs)


text:
# test ModelSetup/Medium/initwf.jl
include("typedefine.jl")
include("mesh.jl")
include("initwf.jl")
include("BD.jl")

# 2d elastic #
vp = [5000.0,5300,5800,6000]
vs = [3500.0,3200,4000,4500]
rho = 2.8
depth = [700,200,400,600]
HX = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,vs,rho,depth,HX,pkf,T,ext,iflag,dx,dz,dt,nT)

elwf = initnelwf(Model.nDZ,Model.nHX,ext,iflag)

pml = BD(elwf,Model)

# 2d acoustic #
vp = [5000.0,5300,5800,6000]
rho = 2.8
depth = [700,200,400,600]
HX = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,rho,depth,HX,pkf,T,ext,iflag,dx,dz,dt,nT)

acwf = initnacwf(Model.nDZ,Model.nHX,ext,iflag)

pml = BD(acwf,Model)

# 3d elastic #
vp = [5000.0,5300,5800,6000]
vs = [3500.0,3200,4000,4500]
rho = 2.8
depth = [700,200,400,600]
HX = 1000
HY = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dy = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,vs,rho,depth,HX,HY,pkf,T,ext,iflag,dx,dy,dz,dt,nT)

elwf = initnelwf(Model.nDZ,Model.nHX,Model.nHY,ext,iflag)

pml = BD(elwf,Model)

# 3d acoustic #
vp = [5000.0,5300,5800,6000]
rho = 2.8
depth = [700,200,400,600]
HX = 1000
HY = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dy = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,rho,depth,HX,HY,pkf,T,ext,iflag,dx,dy,dz,dt,nT)

elwf = initnacwf(Model.nDZ,Model.nHX,Model.nHY,ext,iflag)

pml = BD(elwf,Model)


text:
# Load the training data for the specified alpha and seed type

    using MAT

function LoadFlowSeedTrain(epsilon,seedtype)
    if seedtype == 1
        seedtag = "1"
        seednodes = " 100 Target Nodes Given"
    elseif seedtype == 2
        seedtag = "2"
        seednodes = " 1% of Target Given"
    elseif seedtype == 3
        seedtag = "2per"
        seednodes = " 2% of Target Given"
    elseif seedtype == 4
        seedtag = "3per"
        seednodes = " 3% of Target Given"
    else
        println("Give a valid seedtype, 1 through 4")
        return
    end

    loadstring = "../FlowSeed_Train_Output/Seed_"*seedtag*"_epsilon_"*string(epsilon)*".mat"

    mat = matread(loadstring)

    stats1 = mat["stats1"]
    stats2 = mat["stats2"]
    stats3 = mat["stats3"]
    stats4 = mat["stats4"]

    # Load the Region sizes
    mat = matread("../Brain_Regions/BrainRegions_Train23.mat")
    Sizes = vec(round.(Int64,mat["RegionSizeTrain"]))

    # Explanation of what's stored in the "stats"
    # stats1[:,regi] = [Ssize;timer;Sstats[4];Rcond;pr;re;f1]

    # Save the Precision, Recall, and F1 scores
    numregions = size(stats1,2)
    prs = zeros(4,numregions)
    res = zeros(4,numregions)
    f1s = zeros(4,numregions)

    # Precision for 4 different versions of PageRank
    prs[1,:] = stats1[5,:]
    prs[2,:] = stats2[5,:]
    prs[3,:] = stats3[5,:]
    prs[4,:] = stats4[5,:]

    res[1,:] = stats1[6,:]
    res[2,:] = stats2[6,:]
    res[3,:] = stats3[6,:]
    res[4,:] = stats4[6,:]

    f1s[1,:] = stats1[7,:]
    f1s[2,:] = stats2[7,:]
    f1s[3,:] = stats3[7,:]
    f1s[4,:] = stats4[7,:]

    return prs, res, f1s, Sizes, seednodes

end

function LoadPPRTrain(alpha,seedtype)

    if seedtype == 1
        seedtag = "1"
    elseif seedtype == 2
        seedtag = "2"
    elseif seedtype == 3
        seedtag = "2per"
    elseif seedtype == 4
        seedtag = "3per"
    else
        println("Give a valid seedtype, 1 through 4")
        return
    end

    ppr_string = "../PPR_Train_Output/Seed_"*seedtag*"_alpha_"*string(round(Int64,alpha*100))*"sets_stats.mat"
    mat = matread(ppr_string)

    stats7 = mat["stats_tol_7"]
    sets7 = mat["sets_tol_7"]
    stats8 = mat["stats_tol_8"]
    sets8 = mat["sets_tol_8"]
    stats9 = mat["stats_tol_9"]
    sets9 = mat["sets_tol_9"]
    stats10 = mat["stats_tol_10"]
    sets10 = mat["sets_tol_10"]

    # Load the Region sizes
    mat = matread("../Brain_Regions/BrainRegions_Train23.mat")

    Sizes = vec(round.(Int64,mat["RegionSizeTrain"]))

    ## Explanation of what's stored in the "stats"
    # % [tol;numel(PR);tPR;cond;pPR;rPR;fPR];
    # % Row 1: epsilon tolerance (ppr locality parameter)
    # % Row 2: size of output
    # % Row 3: runtime
    # % Row 4: conductance of output
    # % Row 5: precision
    # % Row 6: recall
    # % Row 7: f1 score

    ## Save the Precision, Recall, and F1 scores
    numregions = size(stats7,2)
    prs = zeros(4,numregions)
    res = zeros(4,numregions)
    f1s = zeros(4,numregions)

    # Precision for 4 different versions of PageRank
    prs[1,:] = stats7[5,:]
    prs[2,:] = stats8[5,:]
    prs[3,:] = stats9[5,:]
    prs[4,:] = stats10[5,:]

    res[1,:] = stats7[6,:]
    res[2,:] = stats8[6,:]
    res[3,:] = stats9[6,:]
    res[4,:] = stats10[6,:]

    f1s[1,:] = stats7[7,:]
    f1s[2,:] = stats8[7,:]
    f1s[3,:] = stats9[7,:]
    f1s[4,:] = stats10[7,:]

    return prs, res, f1s, Sizes

end

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
"""
_get_number_of_updated_edges(locations_of_change).\n
get the number of edges that are updated as [Nx,Ny,Nz]\n
sum([Nx,Ny,Nz]) == size of the matrix inverse.\n
"""
function _get_number_of_updated_edges(locations_of_change)
    return count(locations_of_change[1]), count(locations_of_change[2]), count(locations_of_change[3])
end

"""
dielectric = _set_dielectric(background_dielectric, constants, number_of_cells).\n
given the background dielectric make the dielectric used for the VIE method.\n
it has the same size as the Green's function operator and is zero padded.\n
for now the background medium is vacuum, but this could be changed to some homogenous dielectric.\n
"""
function _set_dielectric(background_dielectric, constants, number_of_cells)
    T = ComplexF32
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))
    p, q, r = number_of_cells
    #dielectric vacuum
    vac = one(T)
    # dielectric background
    tmpx = (background_dielectric[:σˣ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵˣ] .-1)[:,2:end-1,2:end-1]
    tmpy = (background_dielectric[:σʸ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵʸ] .-1)[2:end-1,:,2:end-1]
    tmpz = (background_dielectric[:σᶻ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵᶻ] .-1)[2:end-1,2:end-1,:]

    patient = zeros(T, dx, dy, dz, 3)
    patient[1:p+2,1:q+1,1:r+1,1] = tmpx
    patient[1:p+1,1:q+2,1:r+1,2] = tmpy
    patient[1:p+1,1:q+1,1:r+2,3] = tmpz

    return (;vac, patient)
end

"""
dielectric = _set_dielectric_gpu(background_dielectric, constants, number_of_cells).\n
given the background dielectric make the dielectric used for the VIE method on the gpu.\n
it has the same size as the Green's function operator and is zero padded.\n
for now the background medium is vacuum, but this could be changed to some homogenous dielectric.\n
"""
function _set_dielectric_gpu(background_dielectric, constants, number_of_cells)
    T = ComplexF32
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))
    p, q, r = number_of_cells
    #dielectric vacuum
    vac = one(T)
    # dielectric
    tmpx = (background_dielectric[:σˣ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵˣ] .-1)[:,2:end-1,2:end-1]
    tmpy = (background_dielectric[:σʸ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵʸ] .-1)[2:end-1,:,2:end-1]
    tmpz = (background_dielectric[:σᶻ]./(1im*constants.ω*constants.ϵ₀) + background_dielectric[:ϵᶻ] .-1)[2:end-1,2:end-1,:]

    pat = zeros(T, dx, dy, dz, 3)
    pat[1:p+2,1:q+1,1:r+1,1] = tmpx
    pat[1:p+1,1:q+2,1:r+1,2] = tmpy
    pat[1:p+1,1:q+1,1:r+2,3] = tmpz
    patient = pat |> cu

    return (;vac, patient)
end

"""
_volume_integral_equation_wrapper_explicit_return(sol, sources, vacuum_electric_field, C, vector_to_set_current_density, S, tmp, fft_tmp, G, A, w, dielectric, a, AToE, Ig, scaling, planfft, planifft, x, p, r, rt, u, v, q, uq).\n
function that is used to make a linear map of f(x) = x-CS^TZx.\n
Where Zx is calculated using the VIE method.\n
"""
function _volume_integral_equation_wrapper_explicit_return(
    sol,
    sources,
    vacuum_electric_field,
    C,
    vector_to_set_current_density,
    S,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    scaling,
    planfft,
    planifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq
)
  _volume_integral_equation_wrapper!(sol, sources, vacuum_electric_field, C, vector_to_set_current_density, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, planfft, planifft, x, p, r, rt, u, v, q, uq)

  return sol
end

function _volume_integral_equation_wrapper_explicit_return_gpu(
    sol,
    sources,
    vacuum_electric_field,
    C,
    vector_to_set_current_density,
    S,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    scaling,
    planfft,
    planifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq
)
    input_vector = vector_to_set_current_density |> cu
    _volume_integral_equation_wrapper!(sol, sources, vacuum_electric_field, C, input_vector, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, planfft, planifft, x, p, r, rt, u, v, q, uq)

return collect(sol)
end

"""
b, norm(b) = get_b_vector(CVector, SVector, Einc).\n
Calculates the vector b for the minimaztion of ||b-Ax||_2^2.\n
Einc is the electric field on the yee grid ordered in a vector by Einc_x, Einc_y, Einc_z.\n
Returns the normalized vector b and a float to rescale the output of the minimization.\n
"""
function get_b_vector(C, S, incident_electric_field)
    b = C.*(Transpose(S)*incident_electric_field);
    nb = norm(b);
    bnormalized = b./nb;

    return bnormalized, nb
end

"""
sources, vacuum_electric_field, a, A, fft_tmp, planfft, planifft = _allocate_memory_for_volume_integral_equation(number_of_cells, number_of_updates)\n
allocates arrays that are used for the volume integral equation method.\n
"""
function _allocate_memory_for_volume_integral_equation(
    number_of_cells,
    number_of_updates
)
    T = ComplexF32
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    xl = (p+2)*(q+1)*(r+1)
    yl = (p+1)*(q+2)*(r+1)
    zl = (p+1)*(q+1)*(r+2)
    tl = xl + yl + zl

    #source values
    sources = zeros(T, sum(number_of_updates))
    #incident electric fields and vector potential
    vacuum_electric_field, vector_potential_1d = (zeros(T, tl) for _ in 1:2)
    #fft temporaries
    fft_tmp, vector_potential_3d = (zeros(T, dx*dy*dz*3) for _ in 1:2)
    vector_potential_3d = reshape(vector_potential_3d, dx, dy, dz, 3)
    #planning for fft and ifft
    planfft  = FFTW.plan_fft!(vector_potential_3d, (1, 2, 3))
    planifft = FFTW.plan_ifft!(vector_potential_3d, (1, 2, 3))

    return sources, vacuum_electric_field, vector_potential_1d, vector_potential_3d, fft_tmp, planfft, planifft
end

"""
sources, vacuum_electric_field, a, A, fft_tmp, planfft, planifft = _allocate_memory_for_volume_integral_equation_gpu(number_of_cells, number_of_updates)\n
allocates arrays that are used for the volume integral equation method on the gpu.\n
"""
function _allocate_memory_for_volume_integral_equation_gpu(
    number_of_cells,
    number_of_updates
)
    T = ComplexF32
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    xl = (p+2)*(q+1)*(r+1)
    yl = (p+1)*(q+2)*(r+1)
    zl = (p+1)*(q+1)*(r+2)
    tl = xl + yl + zl

    #source values
    sources = CUDA.zeros(T, sum(number_of_updates))
    #incident electric fields and vector potential
    vacuum_electric_field, vector_potential_1d = (CUDA.zeros(T, tl) for _ in 1:2)
    #fft temporaries
    fft_tmp, vector_potential_3d = (CUDA.zeros(T, dx*dy*dz*3) for _ in 1:2)
    vector_potential_3d = reshape(vector_potential_3d, dx, dy, dz, 3)
    #planning for fft and ifft
    planfft  = FFTW.plan_fft!(vector_potential_3d, (1, 2, 3))
    planifft = FFTW.plan_ifft!(vector_potential_3d, (1, 2, 3))

    return sources, vacuum_electric_field, vector_potential_1d, vector_potential_3d, fft_tmp, planfft, planifft
end

"""
x, p, r, rt, u, v, q, uq = _allocate_memory_cgs(number_of_cells)\n
allocates the memory for the conjugate gradient squared method to solve the VIE.\n
"""
function _allocate_memory_cgs(number_of_cells)
    T = ComplexF32
    #vectors
    np, nq, nr = number_of_cells
    xl = (np+2)*(nq+1)*(nr+1)
    yl = (np+1)*(nq+2)*(nr+1)
    zl = (np+1)*(nq+1)*(nr+2)
    tl = xl + yl + zl

    x, p, r, rt, u, v, q, uq = (zeros(T, tl) for _ in 1:8)

    return x, p, r, rt, u, v, q, uq
end

"""
x, p, r, rt, u, v, q, uq = _allocate_memory_cgs_gpu(number_of_cells)\n
allocates the memory on the gpu for the conjugate gradient squared method to solve the VIE.\n
"""
function _allocate_memory_cgs_gpu(number_of_cells)
    T = ComplexF32
    #vectors
    np, nq, nr = number_of_cells
    xl = (np+2)*(nq+1)*(nr+1)
    yl = (np+1)*(nq+2)*(nr+1)
    zl = (np+1)*(nq+1)*(nr+2)
    tl = xl + yl + zl

    x, p, r, rt, u, v, q, uq = (CUDA.zeros(T, tl) for _ in 1:8)

    return x, p, r, rt, u, v, q, uq
end

#=
some code not directly required for the method that this package was developed for.
=#
"""
magnetic_soucres, fft_tmp, magnetic_vector_potential_1d, magnetic_vector_potential_3d, planfft, planifft = _allocate_memory_for_volume_integral_equation_magnetic_fields(number_of_cells)\n
allocates arrays that are used for the volume integral equation method.\n
"""
function _allocate_memory_for_volume_integral_equation_magnetic_fields(number_of_cells)
    T = ComplexF32

    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))
    # for H fields
    xl = (p+1)*(q+2)*(r+2)
    yl = (p+2)*(q+1)*(r+2)
    zl = (p+2)*(q+2)*(r+1)

    tl = xl + yl + zl

    #incident electric fields
    magnetic_soucres = zeros(T, tl)

    #vector potentials
    magnetic_vector_potential_1d = zeros(T, tl)

    fft_tmp = zeros(T, dx*dy*dz*3)
    magnetic_vector_potential_3d = zeros(T, dx, dy, dz, 3)
    #planning for fft and ifft
    planfft = plan_fft!(magnetic_vector_potential_3d, (1, 2, 3))
    planifft = plan_ifft!(magnetic_vector_potential_3d, (1, 2, 3))

    return magnetic_soucres, fft_tmp, magnetic_vector_potential_1d, magnetic_vector_potential_3d, planfft, planifft
end

"""
magnetic_soucres, fft_tmp, magnetic_vector_potential_1d, magnetic_vector_potential_3d, planfft, planifft = _allocate_memory_for_volume_integral_equation_magnetic_fields_gpu(number_of_cells)\n
allocates arrays that are used for the volume integral equation method.\n
"""
function _allocate_memory_for_volume_integral_equation_magnetic_fields_gpu(number_of_cells)
    T = ComplexF32

    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))
    # for H fields
    xl = (p+1)*(q+2)*(r+2)
    yl = (p+2)*(q+1)*(r+2)
    zl = (p+2)*(q+2)*(r+1)

    tl = xl + yl + zl

    #incident electric fields
    magnetic_soucres = CUDA.zeros(T, tl)

    #vector potentials
    magnetic_vector_potential_1d = CUDA.zeros(T, tl)

    #fft temporaries
    fft_tmp = CUDA.zeros(T, dx*dy*dz*3)
    magnetic_vector_potential_3d = CUDA.zeros(T, dx, dy, dz, 3)

    #planning for fft and ifft
    planfft = plan_fft!(magnetic_vector_potential_3d, (1, 2, 3))
    planifft = plan_ifft!(magnetic_vector_potential_3d, (1, 2, 3))

    return magnetic_soucres, fft_tmp, magnetic_vector_potential_1d, magnetic_vector_potential_3d, planfft, planifft
end


text:
"""
    VortexRing(r1, r2, r3, r4, r_c, n̂, ε)

A vortex ring consisting of four points ``r_i, i = 1,…,4``, a collocation point ``r_c``, a normal vector ``n̂``, and a core size ``ε``.
"""
struct VortexRing{T <: Real} <: AbstractVortex
    r1                :: SVector{3,T}
    r2                :: SVector{3,T}
    r3                :: SVector{3,T}
    r4                :: SVector{3,T}
    collocation_point :: SVector{3,T}
    normal            :: SVector{3,T}
    core              :: T
end

function vortex_lines(p1, p2, p3, p4)
    v1 = quarter_point(p1, p2)
    v4 = quarter_point(p4, p3)
    v2 = v1 + p2 - p1
    v3 = v4 + p3 - p4

    v1, v2, v3, v4
end

"""
Constructor for vortex rings on a Panel3D using Lines. The following convention is adopted:

```
    p1 —front leg→ p4
    |               |
left leg       right leg
    ↓               ↓
    p2 —back leg-→ p3
```
"""
function VortexRing(panel :: Panel3D{T}, normal = panel_normal(panel)) where T <: Real
    r_c = collocation_point(panel)
    ε   = 0. # (norm ∘ average_chord)(panel)
    VortexRing{T}(panel.p1, panel.p2, panel.p3, panel.p4, r_c, normal, ε)
end

Base.length(::VortexRing) = 1

"""
Computes the induced velocities at a point ``r`` of a VortexRing with constant strength ``Γ``.
"""
velocity(r, vortex_ring :: VortexRing, Γ) = sum_vortices(r, structtolist(vortex_ring), Γ)