text:
# promote_type/promote_typejoin should work with Identity and Const (aka Either, but also for Option)

# promote_type
# ============

# promote_type only has to deal with concrete types, including Unions in our case

# promote_type Const & Identity
# -----------------------------

Base.promote_rule(::Type{C}, ::Type{I}) where {C<:Const, I<:Identity} = Union{C, I}

# promote_type Const & Either
# -----------------------------

Base.promote_rule(::Type{Const}, ::Type{Either{L, R}}) where {L, R} = Either{<:Any, R}
Base.promote_rule(::Type{Const{L}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}
Base.promote_rule(::Type{Const{L1}}, ::Type{Either{L2, R}}) where {L1, R, L2 <: L1} = Either{L1, R}
Base.promote_rule(::Type{Const{L2}}, ::Type{Either{L1, R}}) where {L1, R, L2 <: L1} = Either{L1, R}
Base.promote_rule(::Type{Const{L2}}, ::Type{Either{L1, R}}) where {L1, R, L2} = Either{<:Any, R}

Base.promote_rule(::Type{Const}, ::Type{Either{L, <:Any}}) where {L} = Either
Base.promote_rule(::Type{Const{L}}, ::Type{Either{L, <:Any}}) where {L} = Either{L, <:Any}
Base.promote_rule(::Type{Const{L1}}, ::Type{Either{L2, <:Any}}) where {L1, L2 <: L1} = Either{L1, <:Any}
Base.promote_rule(::Type{Const{L2}}, ::Type{Either{L1, <:Any}}) where {L1, L2 <: L1} = Either{L1, <:Any}
Base.promote_rule(::Type{Const{L1}}, ::Type{Either{L2, <:Any}}) where {L1, L2} = Either

Base.promote_rule(::Type{Const}, ::Type{Either}) = Either

# promote_type Identity & Either
# -----------------------------

Base.promote_rule(::Type{Identity}, ::Type{Either{L, R}}) where {L, R} = Either{L, <:Any}
Base.promote_rule(::Type{Identity{R}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}
Base.promote_rule(::Type{Identity{R1}}, ::Type{Either{L, R2}}) where {L, R1, R2 <: R1} = Either{L, R1}
Base.promote_rule(::Type{Identity{R2}}, ::Type{Either{L, R1}}) where {L, R1, R2 <: R1} = Either{L, R1}
Base.promote_rule(::Type{Identity{R2}}, ::Type{Either{L, R1}}) where {L, R1, R2} = Either{L, <:Any}

Base.promote_rule(::Type{Identity}, ::Type{Either{<:Any, R}}) where {R} = Either
Base.promote_rule(::Type{Identity{R}}, ::Type{Either{<:Any, R}}) where {R} = Either{<:Any, R}
Base.promote_rule(::Type{Identity{R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2 <: R1} = Either{<:Any, R1}
Base.promote_rule(::Type{Identity{R2}}, ::Type{Either{<:Any, R1}}) where {R1, R2 <: R1} = Either{<:Any, R1}
Base.promote_rule(::Type{Identity{R2}}, ::Type{Either{<:Any, R1}}) where {R1, R2} = Either

Base.promote_rule(::Type{Identity}, ::Type{Either}) = Either

# promote_type Either & Either
# -----------------------------

# It turns out for this symmetric case we need to define both versions of promote_rule, including the flipped one
# This is needed, because `promote_type` assumes promote_rule to fallback to the general fallback returning Union{}
# for its feature "you only need to define one" to work out correctly.

# It turns out further, that there is a major bug in Julia, not allowing for these kinds of dispatch as of now
# see https://github.com/JuliaLang/julia/issues/36804 the issue and possible future updates

# As a workaround, we overload `promote_type` to deal with our symmetric case
function Base.promote_type(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2, R2}
    T = Either{L1, R1}
    S = Either{L2, R2}
    either_promote_type_fix(T, S)
end

function Base.promote_type(::Type{Either{<:Any, R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2}
    T = Either{<:Any, R1}
    S = Either{<:Any, R2}
    either_promote_type_fix(T, S)
end

function Base.promote_type(::Type{Either{L1, <:Any}}, ::Type{Either{L2, <:Any}}) where {L1, L2}
    T = Either{L1, <:Any}
    S = Either{L2, <:Any}
    either_promote_type_fix(T, S)
end

function either_promote_type_fix(T, S)
    # In addition we need to check for MethodAmbiguities, as we cannot resolve them on one side, we just ignore them
    a = @TryCatch MethodError promote_rule(T, S)
    b = @TryCatch MethodError promote_rule(S, T)

    if isfailure(a) && isfailure(b)
        error("Could not `promote_type{T=$T, S=$S}`, as both `promote_rule(T, S)` and `promote_rule(S, T)` result in
               MethodErrors: $a, $b.")
    elseif issuccess(a) && issuccess(b)
        # We use typeintersect, as we get an Either for the non-fitting site,
        # thanks to
        # * Base.promote_rule(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2, R2} = Either
        # * Base.promote_rule(::Type{Either{<:Any, R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2} = Either
        # * Base.promote_rule(::Type{Either{L1, <:Any}}, ::Type{Either{L2, <:Any}}) where {L1, L2} = Either
        typeintersect(a[], b[])
    elseif issuccess(a)
        a[]
    else
        b[]
    end
end

# we need to be cautious here, as we cannot dispatch on Type{<:Either{<:Any, R}} or similar, because R might not be defined
Base.promote_rule(::Type{Either{L, R}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}
Base.promote_rule(::Type{Either{L1, R}}, ::Type{Either{L2, R}}) where {L1, R, L2 <: L1} = Either{L1, R}
Base.promote_rule(::Type{Either{L2, R}}, ::Type{Either{L1, R}}) where {L1, R, L2 <: L1} = Either{L1, R}  # we needed to add this because rule A and B are ambiguous otherwise
Base.promote_rule(::Type{Either{L, R1}}, ::Type{Either{L, R2}}) where {L, R1, R2 <: R1} = Either{L, R1}
Base.promote_rule(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2 <: L1, R2} = Either{L1, <:Any}
Base.promote_rule(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2, R2 <: R1} = Either{<:Any, R1}
Base.promote_rule(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2 <: L1, R2 <: R1} = Either{L1, R1}
Base.promote_rule(::Type{Either{L2, R1}}, ::Type{Either{L1, R2}}) where {L1, R1, L2 <: L1, R2 <: R1} = Either{L1, R1}  # rule A

Base.promote_rule(::Type{Either{L1, R}}, ::Type{Either{L2, R}}) where {L1, R, L2} = Either{<:Any, R}  # rule B
Base.promote_rule(::Type{Either{L, R1}}, ::Type{Either{L, R2}}) where {L, R1, R2} = Either{L, <:Any}
Base.promote_rule(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2, R2} = Either

Base.promote_rule(::Type{Either{L, <:Any}}, ::Type{Either{L, <:Any}}) where {L} = Either{L, <:Any}
Base.promote_rule(::Type{Either{L1, <:Any}}, ::Type{Either{L2, <:Any}}) where {L1, L2 <: L1} = Either{L1, <:Any}
Base.promote_rule(::Type{Either{L1, <:Any}}, ::Type{Either{L2, <:Any}}) where {L1, L2} = Either

Base.promote_rule(::Type{Either{<:Any, R}}, ::Type{Either{<:Any, R}}) where {R} = Either{<:Any, R}
Base.promote_rule(::Type{Either{<:Any, R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2 <: R1} = Either{<:Any, R1}
Base.promote_rule(::Type{Either{<:Any, R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2} = Either

# NOTE we cannot use `Base.promote_rule(::Type{<:Either}, ::Type{<:Either}) = Either` as apparently this interferes
# with more concrete implementations of promote_rule
# promote_type seem to assume, that you really only define promote_rule for concrete types
Base.promote_rule(::Type{Either}, ::Type{Either}) = Either

# promote_typejoin
# ================

# promote_typejoin never converts types, but always returns the next abstract type in the hierarchy, including Unions
# further, promote_typejoin(A, B) and promote_typejoin(B, A) need both to be defined (unlike promote_rule)

# promote_typejoin Const & Identity
# -----------------------------

Base.promote_typejoin(::Type{C}, ::Type{I}) where {C<:Const, I<:Identity} = Union{C, I}
Base.promote_typejoin(::Type{I}, ::Type{C}) where {C<:Const, I<:Identity} = Union{C, I}

# promote_typejoin Const & Either
# -----------------------------

Base.promote_typejoin(::Type{Const}, ::Type{Either{L, R}}) where {L, R} = Either{<:Any, R}
Base.promote_typejoin(::Type{Either{L, R}}, ::Type{Const}) where {L, R} = Either{<:Any, R}

Base.promote_typejoin(::Type{Const{L}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}
Base.promote_typejoin(::Type{Either{L, R}}, ::Type{Const{L}}) where {L, R} = Either{L, R}

Base.promote_typejoin(::Type{Const{L2}}, ::Type{Either{L1, R}}) where {L1, R, L2} = Either{<:Any, R}
Base.promote_typejoin(::Type{Either{L1, R}}, ::Type{Const{L2}}) where {L1, R, L2} = Either{<:Any, R}

Base.promote_typejoin(::Type{Const}, ::Type{Either{L, <:Any}}) where {L} = Either
Base.promote_typejoin(::Type{Either{L, <:Any}}, ::Type{Const}) where {L} = Either

Base.promote_typejoin(::Type{Const{L}}, ::Type{Either{L, <:Any}}) where {L} = Either{L, <:Any}
Base.promote_typejoin(::Type{Either{L, <:Any}}, ::Type{Const{L}}) where {L} = Either{L, <:Any}

Base.promote_typejoin(::Type{Const{L1}}, ::Type{Either{L2, <:Any}}) where {L1, L2} = Either
Base.promote_typejoin(::Type{Either{L2, <:Any}}, ::Type{Const{L1}}) where {L1, L2} = Either

Base.promote_typejoin(::Type{Const}, ::Type{Either}) = Either
Base.promote_typejoin(::Type{Either}, ::Type{Const}) = Either

# promote_typejoin Identity & Either
# -----------------------------

Base.promote_typejoin(::Type{Identity}, ::Type{Either{L, R}}) where {L, R} = Either{L, <:Any}
Base.promote_typejoin(::Type{Either{L, R}}, ::Type{Identity}) where {L, R} = Either{L, <:Any}

Base.promote_typejoin(::Type{Identity{R}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}
Base.promote_typejoin(::Type{Either{L, R}}, ::Type{Identity{R}}) where {L, R} = Either{L, R}

Base.promote_typejoin(::Type{Identity{R2}}, ::Type{Either{L, R1}}) where {L, R1, R2} = Either{L, <:Any}
Base.promote_typejoin(::Type{Either{L, R1}}, ::Type{Identity{R2}}) where {L, R1, R2} = Either{L, <:Any}

Base.promote_typejoin(::Type{Identity}, ::Type{Either{<:Any, R}}) where {R} = Either
Base.promote_typejoin(::Type{Either{<:Any, R}}, ::Type{Identity}) where {R} = Either

Base.promote_typejoin(::Type{Identity{R}}, ::Type{Either{<:Any, R}}) where {R} = Either{<:Any, R}
Base.promote_typejoin(::Type{Either{<:Any, R}}, ::Type{Identity{R}}) where {R} = Either{<:Any, R}

Base.promote_typejoin(::Type{Identity{R2}}, ::Type{Either{<:Any, R1}}) where {R1, R2} = Either
Base.promote_typejoin(::Type{Either{<:Any, R1}}, ::Type{Identity{R2}}) where {R1, R2} = Either

Base.promote_typejoin(::Type{Identity}, ::Type{Either}) = Either
Base.promote_typejoin(::Type{Either}, ::Type{Identity}) = Either

# promote_type Either & Either
# -----------------------------

# we need to be cautious here, as we cannot dispatch on Type{<:Either{<:Any, R}} or similar, because R might not be defined
Base.promote_typejoin(::Type{Either{L, R}}, ::Type{Either{L, R}}) where {L, R} = Either{L, R}

Base.promote_typejoin(::Type{Either{L1, R}}, ::Type{Either{L2, R}}) where {L1, R, L2} = Either{<:Any, R}

Base.promote_typejoin(::Type{Either{L, R1}}, ::Type{Either{L, R2}}) where {L, R1, R2} = Either{L, <:Any}

Base.promote_typejoin(::Type{Either{L1, R1}}, ::Type{Either{L2, R2}}) where {L1, R1, L2, R2} = Either

Base.promote_typejoin(::Type{Either{L, <:Any}}, ::Type{Either{L, <:Any}}) where {L} = Either{L, <:Any}

Base.promote_typejoin(::Type{Either{L1, <:Any}}, ::Type{Either{L2, <:Any}}) where {L1, L2} = Either

Base.promote_typejoin(::Type{Either{<:Any, R}}, ::Type{Either{<:Any, R}}) where {L, R} = Either{<:Any, R}

Base.promote_typejoin(::Type{Either{<:Any, R1}}, ::Type{Either{<:Any, R2}}) where {R1, R2} = Either

Base.promote_typejoin(::Type{Either}, ::Type{Either}) = Either


text:
using Test

include("../src/01.Julia.Introduction.jl")
include("../src/02.VlasovAmpere.jl")
include("../src/03.RotationGPU.jl")
include("../src/04.ParticleGroup.jl")
include("../src/05.GEMPIC.jl")
include("../src/06.Conclusions.jl")

ext: jl


text:
function ann(w1, w2, w3, x1)
    xx2 = w1 * x1
    x2 = log.(1. + exp.(xx2))   # soft RELU unit
    xx3 = w2 * x2
    x3 = log.(1. + exp.(xx3))   # soft RELU unit
    x4 = w3 * x3
    sum(1. ./ (1. + exp.(-x4)))  # sigmoid output
end

test_compare(ann; w1=rand(10,10), w2=rand(10,10), w3=rand(1,10), x1=rand(10))

f = ann
inputs = [:w1=>rand(10,10), :w2=>rand(10,10), :w3=>rand(1,10), :x1=>rand(10)]


text:
export Services, ShellExecReq

"List of fjåge services."
module Services
  const SHELL = "org.arl.fjage.shell.Services.SHELL"
end

"Shell command execution request message."
const ShellExecReq = MessageClass(@__MODULE__, "org.arl.fjage.shell.ShellExecReq")


text:
module WebIOServeTest
    using WebIO
    using Test

    using Mux
    # Test that webio_serve was defined
    @test isdefined(WebIO, :webio_serve)
    # Test that it was also exported
    @test isdefined(WebIOServeTest, :webio_serve)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1AzureDiskVolumeSource <: SwaggerModel
    cachingMode::Any # spec type: Union{ Nothing, String } # spec name: cachingMode
    diskName::Any # spec type: Union{ Nothing, String } # spec name: diskName
    diskURI::Any # spec type: Union{ Nothing, String } # spec name: diskURI
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiCoreV1AzureDiskVolumeSource(;cachingMode=nothing, diskName=nothing, diskURI=nothing, fsType=nothing, kind=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("cachingMode"), cachingMode)
        setfield!(o, Symbol("cachingMode"), cachingMode)
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("diskName"), diskName)
        setfield!(o, Symbol("diskName"), diskName)
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("diskURI"), diskURI)
        setfield!(o, Symbol("diskURI"), diskURI)
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1AzureDiskVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiCoreV1AzureDiskVolumeSource

const _property_map_IoK8sApiCoreV1AzureDiskVolumeSource = Dict{Symbol,Symbol}(Symbol("cachingMode")=>Symbol("cachingMode"), Symbol("diskName")=>Symbol("diskName"), Symbol("diskURI")=>Symbol("diskURI"), Symbol("fsType")=>Symbol("fsType"), Symbol("kind")=>Symbol("kind"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiCoreV1AzureDiskVolumeSource = Dict{Symbol,String}(Symbol("cachingMode")=>"String", Symbol("diskName")=>"String", Symbol("diskURI")=>"String", Symbol("fsType")=>"String", Symbol("kind")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1AzureDiskVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1AzureDiskVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1AzureDiskVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1AzureDiskVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1AzureDiskVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1AzureDiskVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1AzureDiskVolumeSource)
    (getproperty(o, Symbol("diskName")) === nothing) && (return false)
    (getproperty(o, Symbol("diskURI")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1AzureDiskVolumeSource }, name::Symbol, val)
end


text:
using TimeZones: ParseNextError, TimeZone, localzone, parse_tz_format
using TimeZones: _path_tz_name

# Parse the TZ environment variable format
# Should mirror the behaviour of running:
# `date -u; env TZ="..." date`

@test_throws ParseNextError parse_tz_format("AB")
@test_throws ParseNextError parse_tz_format("+")
@test_throws ParseNextError parse_tz_format("12")

@test parse_tz_format("") == FixedTimeZone("UTC", 0, 0)
@test_throws ParseNextError parse_tz_format("ABC")
@test_throws ParseNextError parse_tz_format("ABC+")
@test_throws ParseNextError parse_tz_format("ABC-")

@test parse_tz_format("ABC1") == FixedTimeZone("ABC", -1 * 3600, 0)
@test parse_tz_format("ABC+1") == FixedTimeZone("ABC", -1 * 3600, 0)
@test parse_tz_format("ABC-1") == FixedTimeZone("ABC", 1 * 3600, 0)

@test parse_tz_format("ABC-24")  == FixedTimeZone("ABC", 24 * 3600, 0)
@test_throws ParseNextError parse_tz_format("ABC-25")
@test_throws ParseNextError parse_tz_format("ABC-100")

@test parse_tz_format("ABC-00:59")  == FixedTimeZone("ABC", 59 * 60, 0)
@test_throws ParseNextError parse_tz_format("ABC-00:99")
@test_throws ParseNextError parse_tz_format("ABC-00:100")

@test parse_tz_format("ABC-00:00:59")  == FixedTimeZone("ABC", 59, 0)
@test_throws ParseNextError parse_tz_format("ABC-00:00:99")
@test_throws ParseNextError parse_tz_format("ABC-00:00:100")

@test_throws ParseNextError parse_tz_format("ABC+00:")
@test_throws ParseNextError parse_tz_format("ABC+12:")

@testset "localzone" begin
    @testset "_path_tz_name" begin
        @test _path_tz_name("") === nothing
        @test _path_tz_name("/tmp/UTC/file") === nothing

        for name in ("UTC", "Europe/Warsaw", "America/Indiana/Indianapolis")
            # Use eval to improve readability of tests when failures occur
            @eval begin
                @test _path_tz_name($name) == $name
                @test _path_tz_name("/usr/share/zoneinfo/$($name)") == $name
                @test _path_tz_name("/var/db/timezone/zoneinfo/$($name)") == $name
            end
        end

        # Verify that the algorithm employed extracts the full time zone name.
        # e.g. Passing in "/usr/share/zoneinfo/Etc/GMT0" may return "GMT0" instead of
        # "Etc/GMT0" if we process from right-to-left.
        @testset "alternative match" begin
            @test _path_tz_name("GMT0") == "GMT0"  # A valid time zone name
            @test _path_tz_name("Etc/GMT0") == "Etc/GMT0"
        end
    end

    # Ensure that the current system's local time zone is supported. If this test fails make
    # sure to report it as an issue.
    @testset "system time zone supported" begin
        @test isa(localzone(), TimeZone)
    end

    # Note: Be careful not to have the tests rely on time zone information being
    # pre-installed on the system. Some minimal systems, such as Docker containers, will not
    # have any system time zone information.
    Sys.isunix() && @testset "TZ environmental variable" begin
        @testset "invalid" begin
            # Bad TZ environment variable formats
            withenv("TZ" => "+12:00") do
                @test_throws ParseNextError localzone()
            end

            # Unable to locate time zone on system
            withenv("TZ" => ":") do
                @test_throws ErrorException localzone()
            end
            withenv("TZ" => ":Etc/Foo") do
                @test_throws ErrorException localzone()
            end

            # Bad TZ which is not a file and not a direct representation. Under these
            # conditions we should throw a error message from parsing the direct
            # representation.
            withenv("TZ" => "A") do
                try
                    localzone()
                catch e
                    @test e isa ParseNextError
                    @test occursin("Unhandled TZ environment variable", e.msg)
                end
            end
        end

        @testset "utc" begin
            withenv("TZ" => "") do
                @test localzone() == utc
            end
            withenv("TZ" => "UTC") do
                @test localzone() == utc
            end
        end

        # When TZ starts with a `:` then the remainder of the string should be treated as a
        # path. For TimeZones.jl we will only load tzfiles when this is specified.
        @testset "force file use" begin
            utc = TimeZone("UTC+0")

            withenv("TZ" => ":UTC+0") do
                @test_throws ErrorException localzone()
            end
            withenv("TZ" => "UTC+0") do
                @test localzone() == utc
            end
        end

        # Use a time zone unrecognized by IANA or TimeZones.jl to verify that the TZDIR
        # environmental variable is being respected.
        @testset "TZDIR environmental variable" begin
            try
                mkdir(joinpath(TZFILE_DIR, "Test"))
                cp(joinpath(TZFILE_DIR, "Etc", "UTC"), joinpath(TZFILE_DIR, "Test", "UTC"))

                @test_throws ArgumentError TimeZone("Test/UTC")
                test_utc = open(joinpath(TZFILE_DIR, "Test", "UTC")) do f
                    TimeZones.read_tzfile(f, "Test/UTC")
                end
                withenv("TZ" => ":Test/UTC", "TZDIR" => TZFILE_DIR) do
                    @test localzone() == test_utc
                end
                withenv("TZ" => "Test/UTC", "TZDIR" => TZFILE_DIR) do
                    @test localzone() == test_utc
                end
            finally
                rm(joinpath(TZFILE_DIR, "Test"), recursive=true, force=true)
            end
        end

        @testset "absolute path" begin
            warsaw_path = joinpath(TZFILE_DIR, "Europe", "Warsaw")
            warsaw_from_file = open(warsaw_path) do f
                TimeZones.read_tzfile(f, "Europe/Warsaw")
            end

            withenv("TZ" => ":" * abspath(warsaw_path)) do
                @test localzone() == warsaw_from_file
            end
            withenv("TZ" => abspath(warsaw_path)) do
                @test localzone() == warsaw_from_file
            end
        end

        @testset "relative path" begin
            warsaw_path = joinpath(TZFILE_DIR, "Europe", "Warsaw")
            warsaw_from_file = open(warsaw_path) do f
                TimeZones.read_tzfile(f, "Europe/Warsaw")
            end
            warsaw = TimeZone("Europe/Warsaw")

            withenv("TZ" => ":Europe/Warsaw", "TZDIR" => TZFILE_DIR) do
                @test localzone() == warsaw_from_file
            end
            withenv("TZ" => "Europe/Warsaw") do
                @test localzone() == warsaw
            end
        end
    end
end


text:
using Statistics
using LinearAlgebra
using Printf

#function nnComputeCosts(nn,datax,datay;print=true,dIdx=1:size(datax,2))
#  c = [ nnCost(datay[:,d],nnForward(nn,datax[:,d])) for d=dIdx ];
#
#  #Print costs
#  if print
#    #for d=1:length(c)
#    #    @printf("dataset %3d: cost=%10.8f\n",d,c[d]);
#    #end
#    @printf("  mean costs: %10.8f\n",mean(c));
#    @printf("median costs: %10.8f\n",median(c));
#    @printf("2-norm costs: %10.8f\n",norm(c,2)/length(c));
#  end
#
#  return c;
#end
#In this version we get the neural net and input data
function nnComputeCosts(nn,datax,datay;print=true,dIdx=1:size(datax,2))
  outputs = nnComputeOutputs(nn,datax;dIdx=dIdx);
  #c = [ nnCost(datay[:,d],outputs[:,d]) for d=dIdx ];

  ##Print costs
  #if print
  #  #for d=1:length(c)
  #  #    @printf("dataset %3d: cost=%10.8f\n",d,c[d]);
  #  #end
  #  @printf("  mean costs: %10.8f\n",mean(c));
  #  @printf("median costs: %10.8f\n",median(c));
  #  @printf("2-norm costs: %10.8f\n",norm(c,2)/length(c));
  #end

  #return c;
  return nnComputeCosts(outputs,datay;print=print,dIdx=dIdx);
end

#In this version we just get the outputs
function nnComputeCosts(outputs,datay;print=true,dIdx=1:size(datay,2))
  c = [ nnCost(datay[:,d],outputs[:,d]) for d=dIdx ];

  #Print costs
  if print
    #for d=1:length(c)
    #    @printf("dataset %3d: cost=%10.8f\n",d,c[d]);
    #end
    @printf("  mean costs: %10.8f\n",mean(c));
    @printf("median costs: %10.8f\n",median(c));
    @printf("2-norm costs: %10.8f\n",norm(c,2)/length(c));
  end

  return c;
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
module GJK

export closest_points
export minimum_distance
export tolerance_verification
export collision_detection

using LinearAlgebra
using StaticArrays

include("support.jl")
include("simplex.jl")
include("proximity.jl")

end # module


text:
module GalleryNLEVP
    # A module which handles interaction with the berlin
    # manchester NLEVP collection

    using MATLAB
    using NonlinearEigenproblems.NEPCore
    using NonlinearEigenproblems.NEPTypes
    using NonlinearEigenproblems.Gallery

    # We have to explicitly specify functions that we want "overload"
    import .NEPCore.compute_Mder
    import .NEPCore.size
    import .NEPCore.compute_Mlincomb

    export NLEVP_NEP

    import .Gallery.nep_gallery
    export nep_gallery

    function fetch_nlevp_path()
        # Try to find NLEVP
        nlevp_path=joinpath(ENV["HOME"],"src","nlevp"); # default ~/src/nlevp
        try
            nlevp_path=ENV["NLEVP_PATH"]
        catch
            # Environment variables was not set
        end

        return nlevp_path;

    end

    """
         NLEVP_NEP represents a NEP in the NLEVP-toolbox
    Example usage: nep=NLEVP_NEP("gun")
    """
    struct NLEVP_NEP <: NEP
        n::Integer
        name::String
        Ai::Array
        installation_path::String
        function NLEVP_NEP(name)
            nlevp_path=fetch_nlevp_path();
            return NLEVP_NEP(name,nlevp_path)
        end
        function NLEVP_NEP(name,nlevp_path)
            if (!isfile(joinpath(nlevp_path,"nlevp.m")))
                error("Unable to find NLEVP when looking in path=",nlevp_path,"  The directory containing  `nlevp.m` should be available in this path. Downloag from https://github.com/ftisseur/nlevp")
            end

            mat"""
            addpath($nlevp_path);
            [$Ai,funs] = nlevp($name);
        """

            this=new(size(Ai[1],1),name,Ai,nlevp_path);
        end
    end

    function compute_Mder(nep::NLEVP_NEP,λ::Number,i::Integer=0)
        lambda=Complex{Float64}(λ)  # avoid type conversion problems
        #println("type",typeof(lambda))
        ## The following commented code is calling nlevp("eval",...)
        ## directly and does not work. We use functions instead
        #        nep_name::String=nep.name
        #        @mput lambda nep_name
        #        if (i==0)
        #            println(λ)
        #            @matlab begin
        #                ll=1+0.1i
        #                M=nlevp("eval",nep_name,lambda
        #            @matlab end
        #            @mget M
        #            return M
        #        else
        #            @matlab begin
        #                (M,Mp)=nlevp("eval",nep_name,lambda)
        #            @matlab end
        #            @mget Mp
        #            return Mp
        #        end
        #    return f,fp
        D=call_current_fun(lambda,i)
        f=D[i+1,:]
        M=zero(nep.Ai[1]);
        for i=1:length(nep.Ai)
            M=M+nep.Ai[i]*f[i]
        end
        return M
    end

    compute_Mlincomb(nep::NLEVP_NEP,λ::Number,V::Union{AbstractMatrix,AbstractVector}, a::Vector) = compute_Mlincomb_from_Mder(nep,λ,V,a)
    compute_Mlincomb(nep::NLEVP_NEP,λ::Number,V::Union{AbstractMatrix,AbstractVector}) = compute_Mlincomb(nep,λ,V, ones(eltype(V),size(V,2)))

    # Return function values and derivatives of the current matlab session "funs"
    # stemming from a previous call to [Ai,funs]=nlevp(nepname).
    # The returned matrix containing derivatives has (maxder+1) rows
    function call_current_fun(lambda,maxder::Integer=0)
        l::ComplexF64=ComplexF64(lambda)  # avoid type problems
        mat"""
    C=cell($maxder+1,1);
    [C{:}]=funs($l);
    $D=cell2mat(C);
    """
        return D
    end

    # size for NLEVP_NEPs
    function size(nep::NLEVP_NEP)
        return (nep.n,nep.n)
    end
    function size(nep::NLEVP_NEP,dim)
        return nep.n
    end

    """
    nep_gallery(NLEVP_NEP, name)
    nep_gallery(NLEVP_NEP, name, nlevp_path)

Loads a NEP from the Berlin-Manchester collection of nonlinear
eigenvalue problems.
"""
    function nep_gallery(::Type{T},name::String) where {T<:NLEVP_NEP}
        nep=NLEVP_NEP(name)
        return nep
    end
    function nep_gallery(::Type{T},name::String,nlevp_path::String) where {T<:NLEVP_NEP}
        nep=NLEVP_NEP(name,nlevp_path)
        return nep
    end

end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__(true)

module AutoHashEquals

export @auto_hash_equals

function auto_hash(name, names)
    _tuple(x) = ntuple(i -> getfield(x, i), length(names))

    quote
        function Base.hash(a::$(name), h::UInt)
            hash($(_tuple)(a), hash(typeof(a), h))
        end
    end
end

function auto_equals(name, names)
    _tuple(x) = ntuple(i -> getfield(x, i), length(names))

    quote
        function Base.:(==)(a::$(name), b::$(name))
            ==($(_tuple)(a), $(_tuple)(b))
        end

        function Base.isequal(a::$(name), b::$(name))
            isequal($(_tuple)(a), $(_tuple)(b))
        end
    end
end

struct UnpackException <: Exception
    msg
end

unpack_name(node::Symbol) = node

function unpack_name(node::Expr)
    if node.head == :macrocall
        unpack_name(node.args[2])
    else
        i = node.head == :type || node.head == :struct ? 2 : 1   # skip mutable flag
        if length(node.args) >= i && isa(node.args[i], Symbol)
            node.args[i]
        elseif length(node.args) >= i && isa(node.args[i], Expr) && node.args[i].head in (:(<:), :(::))
            unpack_name(node.args[i].args[1])
        elseif length(node.args) >= i && isa(node.args[i], Expr) && node.args[i].head == :curly
            unpack_name(node.args[i].args[1])
        else
            throw(UnpackException("cannot find name in $(node)"))
        end
    end
end

macro auto_hash_equals(typ)

    @assert typ.head == :type || typ.head == :struct
    name = unpack_name(typ)

    names = Vector{Symbol}()
    for field in typ.args[3].args
        try
            push!(names, unpack_name(field))
        catch ParseException
            # not a field
        end
    end

    quote
        Base.@__doc__($(esc(typ)))
        $(esc(auto_hash(name, names)))
        $(esc(auto_equals(name, names)))
    end
end

end


text:
fsa = false
all_methods = false

using StaticArrays
@static if fsa
    using FixedSizeArrays
end

if !isdefined(:N)
    N = 4
end

M_f = div(10^9, N^3)
M_g = div(2*10^8, N^2)

# Size

A = rand(Float64,N,N)
As = SMatrix{N,N}(A)
Am = MMatrix{N,N}(A)
@static if fsa
    Af = Mat(ntuple(j -> ntuple(i->A[i,j], N), N)) # there is a bug in FixedSizeArrays Mat constructor (13 July 2016)
end

if !isdefined(:f_mut_marray) || !isdefined(:benchmark_suite) || benchmark_suite == false
    @generated f(n::Integer, A) = :(@inbounds (C = A; for i = 1:n; C = C*A; end; return C))
    @generated f_unrolled{M}(n::Integer, A::Union{SMatrix{M,M},MMatrix{M,M}}) = :(@inbounds (C = A; for i = 1:n; C = StaticArrays.A_mul_B_unrolled(C,A); end; return C))
    @generated f_unrolled_chunks{M}(n::Integer, A::Union{SMatrix{M,M},MMatrix{M,M}}) = :(@inbounds (C = A; for i = 1:n; C = StaticArrays.A_mul_B_unrolled_chunks(C,A); end; return C))
    @generated f_loop{M}(n::Integer, A::Union{SMatrix{M,M},MMatrix{M,M}}) = :(@inbounds (C = A; for i = 1:n; C = StaticArrays.A_mul_B_loop(C,A); end; return C))
    @generated f_via_sarray{M}(n::Integer, A::MMatrix{M,M})= :(@inbounds (C = A; for i = 1:n; C = MMatrix{M,M}(SMatrix{M,M}(C)*SMatrix{M,M}(A)); end; return C))
    @generated f_mut_array(n::Integer, A) = :(@inbounds (C = copy(A); tmp = similar(A); for i = 1:n;  A_mul_B!(tmp, C, A); map!(identity, C, tmp); end; return C))
    @generated f_mut_marray(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B!(tmp, C, A); C.data = tmp.data; end; return C))
    @generated f_mut_unrolled(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_unrolled!(tmp, C, A); C.data = tmp.data; end; return C))
    @generated f_mut_chunks(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_unrolled_chunks!(tmp, C, A); C.data = tmp.data; end; return C))
    @generated f_blas_marray(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; tmp = similar(A); for i = 1:n; StaticArrays.A_mul_B_blas!(tmp, C, A); C.data = tmp.data; end; return C))

    @generated g(n::Integer, A) = :(@inbounds (C = A; for i = 1:n; C = C + A; end; return C))
    @generated g_mut(n::Integer, A) = :(@inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; @inbounds map!(+, C, C, A); end; return C))
    @generated g_via_sarray{M}(n::Integer, A::MMatrix{M,M}) = :(@inbounds (C = similar(A); C[:] = A[:]; for i = 1:n; C = MMatrix{M,M}(SMatrix{M,M}(C) + SMatrix{M,M}(A)); end; return C))

    # Notes: - A[:] = B[:] is allocating in Base, unlike `map!`
    #        - Also, the same goes for Base's implementation of broadcast!(f, A, B, C) (but map! is OK).
    #        - I really need to implement copy() in StaticArrays... (and maybe a special instance of setindex!(C, :))
    #        - I tried generated functions here to see if they force the compiler to specialize more.

end

# Test and compile f()

C = f(2, A)
C_mut = f_mut_array(2, A)
if N <= 4;  @assert C_mut ≈ C; end

println("=====================================")
println("    Benchmarks for $N×$N matrices")
println("=====================================")
@static if all_methods
    print("SMatrix * SMatrix compilation time (unrolled):         ")
    @time eval(quote
        Cs_unrolled = f_unrolled(2, As)
        Cs_unrolled::SMatrix
        if N <= 4; @assert Cs_unrolled ≈ C; end
    end)

    print("SMatrix * SMatrix compilation time (chunks):           ")
    @time eval(quote
        Cs_chunks = f_unrolled_chunks(2, As)
        Cs_chunks::SMatrix
        if N <= 4; @assert Cs_chunks ≈ C; end
    end)

    print("SMatrix * SMatrix compilation time (loop):           ")
    @time eval(quote
        Cs_loop = f_loop(2, As)
        Cs_loop::SMatrix
        if N <= 4; @assert Cs_loop ≈ C; end
    end)

    print("MMatrix * MMatrix compilation time (unrolled):         ")
    @time eval(quote
        Cm_unrolled = f_unrolled(2, Am)
        Cm_unrolled::MMatrix
        if N <= 4; @assert Cm_unrolled ≈ C; end
    end)

    print("MMatrix * MMatrix compilation time (chunks):           ")
    @time eval(quote
        Cm_chunks = f_unrolled_chunks(2, Am)
        Cm_chunks::MMatrix
        if N <= 4; @assert Cm_chunks ≈ C; end
    end)

    print("MMatrix * MMatrix compilation time (loop):           ")
    @time eval(quote
        Cm_loop = f_loop(2, Am)
        Cm_loop::MMatrix
        if N <= 4; @assert Cm_loop ≈ C; end
    end)

    Cm_via_sarray = f_via_sarray(2, Am)
    Cm_via_sarray::MMatrix
    if N <= 4; @assert Cm_via_sarray ≈ C; end
else
    #Cs_unrolled = f_unrolled(2, As)
    #Cs_unrolled::SMatrix
    #if N <= 4; @assert Cs_unrolled ≈ C; end

    #Cs_chunks = f_unrolled_chunks(2, As)
    #Cs_chunks::SMatrix
    #if N <= 4; @assert Cs_chunks ≈ C; end

    #Cm_unrolled = f_unrolled(2, Am)
    #Cm_unrolled::MMatrix
    #if N <= 4; @assert Cm_unrolled ≈ C; end

    #Cm_chunks = f_unrolled_chunks(2, Am)
    #Cm_chunks::MMatrix
    #if N <= 4; @assert Cm_chunks ≈ C; end
end

# Warmup and some checks
Cs = f(2, As)
Cm = f(2, Am)

Cs::SMatrix
if N <= 4; @assert Cs ≈ C; end
Cm::MMatrix
if N <= 4; @assert Cm ≈ C; end

@static if fsa
    @static if all_methods
        print("Mat * Mat compilation time:                            ")
        @time eval(quote
            Cf = f(2, Af)
            Cf::Mat
            if N <= 4; @assert Cf == C; end
        end)
    else
        Cf = f(2, Af)
        Cf::Mat
        if N <= 4; @assert Cf == C; end
    end
end

# Mutating versions

Cm_mut = f_mut_marray(2, Am)
Cm_mut::MMatrix
if N <= 4; @assert Cm_mut ≈ C; end

@static if all_methods
    println()
    print("A_mul_B!(MMatrix, MMatrix) compilation time (unrolled):")
    @time eval(quote
        Cm_mut_unrolled = f_mut_unrolled(2, Am)
        Cm_mut_unrolled::MMatrix
        if N <= 4; @assert Cm_mut_unrolled ≈ C; end
    end)

    print("A_mul_B!(MMatrix, MMatrix) compilation time (chunks):  ")
    @time eval(quote
        Cm_mut_chunks = f_mut_chunks(2, Am)
        Cm_mut_chunks::MMatrix
        if N <= 4; @assert Cm_mut_chunks ≈ C; end
    end)

    print("A_mul_B!(MMatrix, MMatrix) compilation time (BLAS):    ")
    @time eval(quote
        Cm_blas = f_blas_marray(2, Am)
        Cm_blas::MMatrix
        if N <= 4; @assert Cm_blas ≈ C; end
    end)
else
    #Cm_mut_unrolled = f_mut_marray(2, Am)
    #Cm_mut_unrolled::MMatrix
    #if N <= 4; @assert Cm_mut_unrolled ≈ C; end

    #Cm_mut_chunks = f_mut_chunks(2, Am)
    #Cm_mut_chunks::MMatrix
    #if N <= 4; @assert Cm_mut_chunks ≈ C; end

    #Cm_blas = f_blas_marray(2, Am)
    #Cm_blas::MMatrix
    #if N <= 4; @assert Cm_blas ≈ C; end
end

# Test and compile g()

C = g(2, A)
C_mut = g_mut(2, A)

if N <= 4; @assert C_mut ≈ C; end

Cs = g(2, As)
Cm = g(2, Am)
Cm_mut = g_mut(2, Am)

Cs::SMatrix
if N <= 4; @assert Cs == C; end
Cm::MMatrix
if N <= 4; @assert Cm == C; end
Cm_mut::MMatrix
if N <= 4; @assert Cm_mut == C; end

@static if all_methods
    Cm_via_sarray = g_via_sarray(2, Am)
    Cm_via_sarray::MMatrix
    if N <= 4; @assert Cm_via_sarray == C; end
end

@static if fsa
    Cf = g(2, Af)
    Cf::Mat
    if N <= 4; @assert Cf == C; end
end

println()

# Do the performance tests

println("Matrix multiplication")
println("---------------------")
begin
   print("Array               ->"); @time f(M_f, A)
   @static if fsa
       print("Mat                 ->"); @time f(M_f, Af)
   end
   print("SArray              ->"); @time f(M_f, As)
   print("MArray              ->"); @time f(M_f, Am)
   @static if all_methods
       print("SArray (unrolled)   ->"); @time f_unrolled(M_f, As)
       print("SArray (chunks)     ->"); @time f_unrolled_chunks(M_f, As)
       print("SArray (loop)       ->"); @time f_loop(M_f, As)
       print("MArray (unrolled)   ->"); @time f_unrolled(M_f, Am)
       print("MArray (chunks)     ->"); @time f_unrolled_chunks(M_f, Am)
       print("MArray (loop)       ->"); @time f_loop(M_f, Am)
       print("MArray (via SArray) ->"); @time f_via_sarray(M_f, Am)
   end
end
println()

println("Matrix multiplication (mutating)")
println("--------------------------------")
begin
   print("Array               ->"); @time f_mut_array(M_f, A)
   print("MArray              ->"); @time f_mut_marray(M_f, Am)
   @static if all_methods
       print("MArray (unrolled)   ->"); @time f_mut_unrolled(M_f, Am)
       print("MArray (chunks)     ->"); @time f_mut_chunks(M_f, Am)
       print("MArray (BLAS gemm!) ->"); @time f_blas_marray(M_f, Am)
   end
end
println()

println("Matrix addition")
println("---------------")
begin
   print("Array               ->"); @time g(M_g, A)
   @static if fsa
       print("Mat                 ->"); @time g(M_g, Af)
   end
   print("SArray              ->"); @time g(M_g, As)
   print("MArray              ->"); @time g(M_g, Am)
   @static if all_methods
       print("MArray (via SArray) ->"); @time g_via_sarray(M_g, Am)
   end
end
println()

println("Matrix addition (mutating)")
println("--------------------------")
begin
   print("Array  ->"); @time g_mut(M_g, A) # broadcast! seems to be broken!
   print("MArray ->"); @time g_mut(M_g, Am)
end
println()


text:
Base.@kwdef struct LocalDIMParameters
    compression=Matrix
    atol_near::Float64=0.1
    rtol_near::Float64=5
    nb_multiplier::Int=3
    r_multiplier::Float64=10
    aux_surf_multiplier::Float64=5
end

function assemble_ldim(iop::IntegralOperator,p::LocalDIMParameters=LocalDIMParameters())
    X = target_surface(iop)
    Y = source_surface(iop)
    out        = p.compression(iop)
    correction = assemble_ldim_correction(iop,p)
    axpy!(1,correction,out) # out <-- out + correction
end

function assemble_ldim_correction(iop::IntegralOperator,p::LocalDIMParameters)
    X    = target_surface(iop)
    Y    = source_surface(iop)
    pde  = iop.kernel.op
    dict_near = near_interaction_list(qnodes(X),Y;atol=p.atol_near)
    T = eltype(iop)
    # vectors to store nonzero indices (i,j) and value
    Is = Int[]
    Js = Int[]
    Vs = T[]
    for (E,list_near) in dict_near
        _assemble_ldim_correction!(Is,Js,Vs,iop,E,list_near,p)
    end
    Sp = sparse(Is,Js,Vs,size(iop)...)
    return Sp
end

@noinline function _assemble_ldim_correction!(Is,Js,Vs,iop,E,list_near,p::LocalDIMParameters)
    X    = target_surface(iop)
    Y    = source_surface(iop)
    T    = eltype(iop)
    k    = iop.kernel
    α,β = combined_field_coefficients(k)
    op   = iop.kernel.op
    if kernel_type(iop) isa Union{SingleLayer,DoubleLayer}
        K0 = SingleLayerKernel(op)
        K1 = DoubleLayerKernel(op)
        need_normal = false
    elseif kernel_type(iop) isa Union{AdjointDoubleLayer,HyperSingular}
        K0 = AdjointDoubleLayerKernel(op)
        K1 = HyperSingularKernel(op)
        need_normal = true
    else
        error()
    end
    el2qnodes = elt2dof(Y,E)
    els  = elements(Y,E)
    @assert length(els) == length(list_near)
    ni      = size(el2qnodes,1)
    nb      = p.nb_multiplier*ni
    M       = Matrix{T}(undef,2*ni,nb)
    I       = zeros(T,nb)
    for nel in 1:length(els)
        el    = els[nel]
        tags  = @view el2qnodes[:,nel]
        # interpolation nodes on source element
        yi = @view Y.qnodes[tags]
        νi = @view Y.qnormals[tags]
        wi = @view Y.qweights[tags]
        # basis
        l       = sum(wi)
        h       = l/length(wi) # typical space of a grid
        γ₀B,γ₁B = basis_ldim(op,el,p.r_multiplier*l,nb)
        # interpolation matrix
        interpolation_matrix!(M,yi,νi,γ₀B,γ₁B)
        # F   = qr!(M)
        D   = [diagm(α*ones(length(yi))) ; diagm(β*ones(length(yi)))]
        # now loop over target points in X which need correction when integrated
        # over el
        for (iglob,jloc) in list_near[nel]
            x   = X.qnodes[iglob]
            y₀  = yi[jloc]
            ν₀  = νi[jloc]
            r   = x-y₀
            d   = norm(r)
            d > p.rtol_near*h && continue
            side = dot(r,ν₀)
            location = d == 0 ? :onsurface : side > 0 ? :outside : :inside
            if location == :onsurface
                Γ = auxiliary_domain_ldim(el,-p.aux_surf_multiplier*h)
                σ = -0.5
            elseif location == :outside
                Γ = auxiliary_domain_ldim(el,-h)
                σ = 0.0
            elseif location == :inside
                Γ = auxiliary_domain_ldim(el,h)
                σ = 0.0
            end
            if K0 isa SingleLayerKernel
                _integrate_basis_ldim!(I,Γ,K0,K1,γ₀B,γ₁B,x,σ) # modify r
            elseif K0 isa AdjointDoubleLayerKernel
                νx  = X.qnormals[iglob]
                _integrate_basis_ldim!(I,Γ,K0,K1,γ₀B,γ₁B,x,νx,σ) # modify r
            else
                error()
            end
            w   = (transpose(I)/M)*D
            v   = transpose(w) - iop[iglob,tags] # remove the "naive" approximation
            append!(Is,fill(iglob,length(w)))
            append!(Js,tags)
            append!(Vs,v)
        end
    end
end

"""
    singular_weights_ldim(k,τ,yi,x)

Given a kernel `k`, generate quadrature weights `w` such that
```math
\\int_\\tau k(x,y) \\varphi(y) ds\\_y \\approx \\sum_q \\varphi(y_q) w_q
```
"""
@inline function singular_weights_ldim(k::AbstractKernel,el,yi,νi,x,side=:onsurface)
    h = 0.1
    # pre-compute things which do not depend on x
    op  = k.op
    T   = return_type(k)
    K0  = SingleLayerKernel(op)
    K1  = DoubleLayerKernel(op)
    nb  = 3*length(yi)
    γ₀B,γ₁B = basis_ldim(op,el,10*h,nb)
    M       = Matrix{T}(undef,2*length(yi),nb)
    interpolation_matrix!(M,yi,νi,γ₀B,γ₁B)
    F   = qr(M)
    # compute things which depend on target point x
    if side == :onsurface
        Γ = auxiliary_domain_ldim(el,-h)
        σ = -0.5
    elseif side == :outside
        Γ = auxiliary_domain_ldim(el,-h)
        σ = 0.0
    elseif side == :inside
        Γ = auxiliary_domain_ldim(el,h)
        σ = 0.0
    end
    I   = zeros(T,nb)
    _integrate_basis_ldim!(I,Γ,K0,K1,γ₀B,γ₁B,x,σ) # modify r
    α,β = combined_field_coefficients(k)
    D   = [diagm(α*ones(length(yi))) ; diagm(β*ones(length(yi)))]
    w   = ((transpose(I)/F.R)*adjoint(F.Q))*D
    return w
end

function _singular_weights_ldim_kernel(el,K0,K1,γ₀B,γ₁B,x,side)
    nb = length(γ₀B)
    if side == :onsurface
        Γ = auxiliary_domain_ldim(el,-h)
        σ = -0.5
    elseif side == :outside
        Γ = auxiliary_domain_ldim(el,-h)
        σ = 0.0
    elseif side == :inside
        Γ = auxiliary_domain_ldim(el,h)
        σ = 0.0
    end
    I   = zeros(T,nb)
    _integrate_basis_ldim!(I,Γ,K0,K1,γ₀B,γ₁B,x,σ) # modify r
    α,β = combined_field_coefficients(k)
    D   = [diagm(α*ones(length(yi))) ; diagm(β*ones(length(yi)))]
    w   = ((transpose(I)/F.R)*adjoint(F.Q))*D
    return w
end

function interpolation_matrix!(M,yi,νi,γ₀B,γ₁B)
    nb = length(γ₀B)
    for j in 1:nb
        for i in 1:length(yi)
            M[i,j]            = γ₀B[j](yi[i])
            M[length(yi)+i,j] = γ₁B[j](yi[i],νi[i])
        end
    end
    return M
end

function _integrate_basis_ldim!(r,Γ,K0::SingleLayerKernel,K1::DoubleLayerKernel,γ₀B,γ₁B,x,σ)
    nb = length(r)
    for n in 1:nb
        r[n] = sum(Γ) do b
            integrate(ReferenceLine) do v
                y  = b(v)
                νy = normal(b,v)
                μ  = measure(b,v)
                (K0(x,y)*γ₁B[n](y,νy) - K1(x,y,νy)*γ₀B[n](y))*μ
            end
        end
        r[n] = r[n] + σ*γ₀B[n](x)
    end
    return r
end

function _integrate_basis_ldim!(r,Γ,K0::AdjointDoubleLayerKernel,K1::HyperSingularKernel,γ₀B,γ₁B,x,νx,σ)
    nb = length(r)
    for n in 1:nb
        r[n] = sum(Γ) do b
            integrate(ReferenceLine) do v
                y  = b(v)
                νy = normal(b,v)
                μ  = measure(b,v)
                (K0(x,y,νx)*γ₁B[n](y,νy) - K1(x,y,νx,νy)*γ₀B[n](y))*μ
            end
        end
        r[n] = r[n] + σ*γ₁B[n](x,νx)
    end
    return r
end

function _integrate_basis_ldim_static!(r,Γ,K0,K1,γ₀B,γ₁B,x,σ,q,qs)
    nb = length(r)
    for n in 1:nb
        x̂,ŵ = q()
        b = Γ[2]
        r[n] = sum(zip(x̂,ŵ)) do (x̂,ŵ)
                y  = b(x̂)
                νy = normal(b,x̂)
                μ  = measure(b,x̂)
                (K0(x,y)*γ₁B[n](y,νy) - K1(x,y,νy)*γ₀B[n](y))*μ*ŵ
        end
        b = Γ[1]
        x̂,ŵ = qs()
        r[n] += sum(zip(x̂,ŵ)) do (x̂,ŵ)
                y  = b(x̂)
                νy = normal(b,x̂)
                μ  = measure(b,x̂)
                (K0(x,y)*γ₁B[n](y,νy) - K1(x,y,νy)*γ₀B[n](y))*μ*ŵ
        end
        b = Γ[3]
        x̂,ŵ = qs()
        x̂ = -x̂ .+ 1 |> reverse
        ŵ = reverse(ŵ)
        r[n] += sum(zip(x̂,ŵ)) do (x̂,ŵ)
                y  = b(x̂)
                νy = normal(b,x̂)
                μ  = measure(b,x̂)
                (K0(x,y)*γ₁B[n](y,νy) - K1(x,y,νy)*γ₀B[n](y))*μ*ŵ
        end
        r[n] = r[n] + σ*γ₀B[n](x)
    end
    return r
end

function basis_ldim(pde,el,r,nb)
    ν  = normal(el,0.5)
    xc = el(0.5) # midpoint of the element
    xs = _circle_sources(nsources=nb, radius=r, center=xc)
    G   = SingleLayerKernel(pde)
    dG  = DoubleLayerKernel(pde)
    γ₀B = map(xs) do x
        y -> G(x,y)
    end
    γ₁B = map(xs) do x
        (y,ν) -> dG(x,y,ν)
    end
    return γ₀B, γ₁B
end

function auxiliary_domain_ldim(el::AbstractElement,h)
    @assert domain(el) == ReferenceLine()
    a,b = el(0),el(1) # the endpoints
    ν   = normal(el,0.5)
    Γ₁  = LagrangeLine(b,b+h*ν)
    Γ₂  = LagrangeLine(b+h*ν,a+h*ν)
    Γ₃  = LagrangeLine(a+h*ν,a)
    return Γ₁,Γ₂,Γ₃
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module DistSystem
using DataFrames

export Substation, DG, System, factory_system, null_der, DER

include("Tools.jl")
using .Tools: make_Y_bus
import DataFrames: DataFrame

struct DG
    bus
    P_limit
    Q_limit
    Cost
end

struct DER
    bus::Any
    S_limit::Float64
    energy::Float64
    alpha::Float64
    beta::Float64
    P_limit::Vector{Float64}
    Q_limit::Vector{Float64}
    scenario::Vector{Float64}
    Cost::Vector{Float64}
end

struct Substation
    nominal_voltage
    bus
    voltage
    P_limit
    Q_limit
    Cost
end

mutable struct System
    Y
    nbuses
    VL
    VH
    PL
    QL
    amp
    m_load
    m_new_dg
    time_curr
    substation
    dgs
    Bsh
    buses
    System() = new()
end

function null_der(buses::Vector)::DER
    return DistSystem.DER(buses[2],
        0.0, 0.0, 0.0, 1.0, [0.0, 0.0], [0.0, 0.0], [0.0], [0.0, 0.0, 0.0]
    )
end

function factory_system(data::DataFrame, VL::Float64, VH::Float64, sub::Substation)
    sys = System()

    ### System's Data
    Sᴺ = 1e6
    Iᴺ = Sᴺ / (sub.nominal_voltage * sqrt(3))
    Zᴺ = (sub.nominal_voltage^2) / Sᴺ
    Yᴺ = 1 / Zᴺ

    sys.time_curr = 0.0
    sys.Y = make_Y_bus(data, sub.nominal_voltage) / Yᴺ
    sys.nbuses = size(sys.Y)[1]
    sys.VL = VL
    sys.VH = VH
    sys.m_load = [1.0]
    sys.m_new_dg = [[0.0]]
    sys.PL = collect(skipmissing(data.P_MW))
    sys.QL = collect(skipmissing(data.Q_MVAr))
    sys.buses = collect(skipmissing(data.Bus))
    sys.Bsh = (-(collect(skipmissing(data.Bshunt_MVAr)) * 1e6) ./ (sub.nominal_voltage^2)) / Yᴺ

    amp_pu = collect(skipmissing(data.Amp_pu))
    fb = collect(skipmissing(data.FB))
    tb = collect(skipmissing(data.TB))
    sys.amp = Dict((fb[k], tb[k]) => amp_pu[k] for k = 1:length(amp_pu))
    null_dg = null_der(sys.buses)
    sys.dgs = [null_dg]

    sys.substation = sub

    return sys

end

end


text:
using Bonsai
using Bonsai: combine_middleware, Middleware
using URIs
using Dates
using HTTP: Request

t = false
c = false

@testset "combine_middleware" begin

	function timer(stream, next)
		x = now()
		next(stream)
		elapsed = x - now()
		global t
		t = true
	end

	function cors(stream, next)
		next(stream)
		global c
		c = true
	end

	fn = combine_middleware([timer, cors ])
	fn(nothing)
	@test c && t
	@test combine_middleware([])(true)
end

app = App()

app.get("**") do stream, next
end

req = Request()
req.method = "GET"
req.target = "/"
match(app, req)

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
abstract type AbstractSimulatedAnnealing <: AbstractAlgorithm end
abstract type AbstractCooling end

"""
    struct FixedCooling <: AbstractCooling
        temperature::Float64
    end

Using this cooling scheduler, the algorithm accepts a worser solution with probability ``\\alpha^{-\\Delta}``, where
``\\alpha`` is the parameter temperature and ``\\Delta`` is the absolute difference of the fitness functions between the parent and child.

Through `FixedCooling`, you are basicaully using `Metropolis` algorithm.
"""

struct FixedCooling <: AbstractCooling
    temperature::Real
    function FixedCooling(temperature::Real = 0.0)
        new(temperature)
    end
end

function temperature(cooling::FixedCooling, iter::Int)
    return cooling.temperature
end

struct FixedCoolingTag <: AbstractCooling
    temperature::Real
    order::LaTeXString
    function FixedCoolingTag(temperature::Real = 0.0, order::LaTeXString = L"nothing")
        new(temperature, order)
    end
end

function temperature(cooling::FixedCoolingTag, temp::Real)
    return cooling.temperature
end

"""
    struct SimulatedAnnealing <: AbstractSimulatedAnnealing
        cooling::AbstractCooling
        stop::AbstractStop
        name::LaTeXString
        temperature::Float64
    end
"""
struct SimulatedAnnealing <: AbstractSimulatedAnnealing
    cooling::AbstractCooling
    stop::AbstractStop
    name::LaTeXString
    temperature::Real
    mutation::Mutation

    function SimulatedAnnealing(
        cooling::AbstractCooling;
        stop::AbstractStop = FixedBudget(10^11),
        name::LaTeXString = L"Simulated-Annealing",
        temperature::Real = -1.0,
        mutation::Mutation = KBitFlip(1)
    )
        new(cooling, stop, name, temperature, mutation)
    end
end

function optimize(x, setting::SimulatedAnnealing)
    trace = Trace(setting, x)
    x = initial(x)
    n = length(x)

    for iter ∈ 1:niterations(setting.stop)
        α = float(temperature(setting.cooling, iter))
        y = mutation(x, setting.mutation)
        Δ = fitness(y) - fitness(x)

        if Δ ≥ 0
            if Δ > 0
                if isoptimum(y)
                    record(trace, y, iter, isoptimum(y))
                end
            end
            x = y
            if isoptimum(x)
                return trace
            end
        elseif rand() ≤ α^(Δ)
            # record(trace, y, iter, isoptimum(y))
            x = y
        end
    end
    trace
end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
# Parents: use findall in adjacency matrix, instead of `Graphs.innerneigbhors` method:
children(mat::SparseMatrixCSC, node::Int64) = findall(x -> x != 0, mat[node, :])
children(net::AbstractSimpleWeightedGraph, node::Int64) = Graphs.inneighbors(net, node)

# Parents: use findall in adjacency matrix, instead of `Graphs.outneigbhors` method:
parents(mat::SparseMatrixCSC, node::Int64) = findall(x -> x != 0, mat[:, node])
parents(net::AbstractSimpleWeightedGraph, node::Int64) = Graphs.outneighbors(net, node)

# Parents: use findall in adjacency matrix, instead of `Graphs.all_neigbhors` method:
neighbors(mat::SparseMatrixCSC, node::Int64) = union(children(mat, node), parents(mat, node))
neighbors(net::AbstractSimpleWeightedGraph, node::Int64) = Graphs.all_neighbors(net, node)

function toplogicalsort(net::AbstractSimpleWeightedGraph; method = "graphs")
    if Graphs.is_directed(net)
        if method == "dfs"
            L, acyclic = _topologicalsort_dfs(net)
        elseif method == "kahn"
            L, acyclic = _topologicalsort_kahn(net)
        else
            try # Todo: ooh, that's not valid in all error cases!
                acyclic = true
                L = Graphs.topological_sort_by_dfs(net)
            catch e
                L = []
                acyclic = false
            end
        end
        return L, acyclic
    else
        print("undirected graph or cyclic graph. No topological order")
        return [], false
    end
end

function _topologicalsort_kahn(g::AbstractSimpleWeightedGraph)

    # if there exists a topological order of directed graph g, then it is acyclic
    acyclic = true

    adj = g.weights
    # get set $S$ of all nodes with no incoming edge
    n = nv(g)
    S = Int64[]
    num_parents = zeros(Int, 1, n) # no incoming EDGELABELSIZE
    for i = 1:n
        num_parents[i] = length(WeightedNetwork.parents(adj, i))
        if num_parents[i] == 0
            push!(S, i)
        end
    end

    L = [] # return value

    while !isempty(S)
        node = pop!(S)
        push!(L, node)
        childs = WeightedNetwork.children(adj, node)
        for kid in childs
            adj[node, kid] = 0
            if length(WeightedNetwork.parents(adj, kid)) == 0
                push!(S, kid)
            end
        end
    end

    if sum(adj) == 0.0 || sum(adj) == 0
        acyclic = true
        return L, acyclic
    else
        acyclic = false
        print("No tolpological order. This is not a DAG. The direct graph is cyclic!")
        return [], acyclic
    end
end

function _topologicalsort_dfs(g::AbstractSimpleWeightedGraph)

    # if there exists a topological order of directed graph g, then it is acyclic
    acyclic = true

    mat = g.weights
    L = Int64[]
    n = nv(g)
    marked = zeros(Int, 1, n)

    # choose node with most children:
    unmarked_node = marked[1]

    function visit(node)
        if marked[node] == 2 # 2 = permanent_mark
            return
        elseif marked[node] == 1 # 1 = temporary mark
            acyclic = false # then this unacyclic and not a DAG.
            print("No tolpological order. This is not a DAG. The direct graph is cyclic!")
            return [], acyclic
        end

        marked[node] = 1

        childs = WeightedNetwork.children(g, node)

        for kid in childs
            visit(kid)
        end

        marked[node] = 2
        push!(L, node)
    end

    i = 1
    while !isnothing(unmarked_node)
        unmarked_node = findfirst(isequal(0), marked)
        if !isnothing(unmarked_node)
            visit(unmarked_node[2])
        end
    end

    return reverse(L), acyclic
end

# 8-May-1998	 4:44 PM	ATC	Created under MATLAB 5.1.0.421
# ATC = Ali Taylan Cemgil,
# SNN - University of Nijmegen, Department of Medical Physics and Biophysics
# modified for Julia, Röttgermann, 12/2021
# Test-only
function _toposort(g::AbstractSimpleWeightedGraph)
    adj = g.weights
    N = nv(g)
    indeg = sum(g, 2)
    outdeg = sum(g, 1)

    seq = Int64[]

    for i = 1:N
        idx = findall(x -> x == 0, indeg)

        if isempty(idx)
            seq = []
            break
        end

        max, idx_max = findmax(outdeg[idx])
        indx = idx[idx_max]
        push!(seq, indx)

        indeg[indx] = -1
        idx = findall(x -> x > 0, adj[indx, :])
        !isempty(idx) ? indeg[idx] = indeg[idx] .- 1 : nothing
        @show seq
    end

    return vec(seq)'

end

function is_tree(mat::SparseMatrixCSC)
    # zusammenhängend
    error("Is not implemented yet!")
    #
end

is_tree(g::SimpleWeightedDiGraph) = is_tree(g.weights)
is_tree(g::SimpleWeightedGraph) = is_tree(g.weights)

function is_symmetric(M::SparseMatrixCSC)
    sz = size(M)

    if sz[1] != sz[2]
        return false
    else

        if M == M'
            return true
        else
            return false
        end
    end
end

is_symmetric(g::AbstractSimpleWeightedGraph) = is_symmetric(g.weights)

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Frames

This module contains functions for **frames** that is coordinate systems in 3D.
The orientation of a frame is described either with a 3x3 **rotation matrix**
or with a **quaternion vector** and its origin is described with a **Vector3D**:

- `const `[`Modia3D.RotationMatrix`](@ref) = SMatrix{3,3,Float64,9}`:
  Type of a Rotation matrix to rotate from a frame 1 into a frame 2.

- `const `[`Modia3D.Quaternion`](@ref) = SVector{4,Float64}`:
  Type of a Quaternion vector to rotate from a frame 1 into a frame 2.

- `const `[`Modia3D.Vector3D`](@ref) = SVector{3,Float64}`:
  Type of a vector in 3D space (e.g. position vector of the origin of a frame).

The following constants are defined

- `const `[`Modia3D.NullRotation`](@ref):
  RotationMatrix with no rotation from a frame 1 into a frame 2.

- `const `[`Modia3D.NullQuaternion`](@ref):
  Quaternion vector with no rotation from a frame 1 into a frame 2.

- `const `[`Modia3D.ZeroVector3D`](@ref):
  Vector3D with only zero elements.

If an angle is given as an argument to one of the functions below, it might be a
number (interpreted as having unit `rad`) or a number with a unit
(for example: `using Unitful; angle = 90u"°"`).

# Constructors for a RotationMatrix R

The following functions return a [`Modia3D.RotationMatrix`](@ref)` R`
to rotate a frame 1 into a frame 2.

| Function                                         | Description                                        |
|:-------------------------------------------------|:---------------------------------------------------|
| [`Modia3D.rot1`](@ref)(angle)                    | Rotate around `angle` along x-axis                 |
| [`Modia3D.rot2`](@ref)(angle)                    | Rotate around `angle` along y-axis                 |
| [`Modia3D.rot3`](@ref)(angle)                    | Rotate around `angle` along z-axis                 |
| [`Modia3D.rot123`](@ref)(angle1, angle2, angle3) | Rotate around angles along x,y,z-axes              |
| [`Modia3D.rotAxis`](@ref)(axis,angle)            | Rotate around `angle` along `axis` (= 1,2,3)       |
| [`Modia3D.rotAxis`](@ref)(axis,positive,angle)   | Rotate around `angle` if `positive`, else `-angle` |
| [`Modia3D.rot_e`](@ref)(e, angle)                | Rotate around `angle` along unit vector `e`        |
| [`Modia3D.rot_nxy`](@ref)(nx, ny)                | `nx`/`ny` are in x/y-direction of frame 2          |
| [`Modia3D.from_q`](@ref)(q)                      | Return `R` from Quaternion `q`                     |

# Constructors for a Quaternion q

The following functions return a [`Modia3D.Quaternion`](@ref)` q`
to rotate a frame 1 into a frame 2.
Since `q` and `-q` define the same rotation the constructor functions have
a keyword argument `q_guess::Quaternion = NullQuaternion`.
From the two possible solutions `q` the one is returned that is closer to `q_guess`.

| Function                                          | Description                                |
|:--------------------------------------------------|:-------------------------------------------|
| [`Modia3D.qrot1`](@ref)(angle)                    | Rotate around `angle` along x-axis         |
| [`Modia3D.qrot2`](@ref)(angle)                    | Rotate around `angle` along y-axis         |
| [`Modia3D.qrot3`](@ref)(angle)                    | Rotate around `angle` along z-axis         |
| [`Modia3D.qrot123`](@ref)(angle1, angle2, angle3) | Rotate around angles along x,y,z-axes      |
| [`Modia3D.qrot_e`](@ref)(e, angle)                | Rotate around `angle` along unit vector `e`|
| [`Modia3D.qrot_nxy`](@ref)(nx, ny)                | `nx`/`ny` are in x/y-direction of frame 2  |
| [`Modia3D.from_R`](@ref)(R)                       | Return `q` from RotationMatrix `R`         |

# Operations on Frames

The following functions provide operations on frames. The orientation of a frame is
defined with argument `Rq` meaning it can be either a
[`Modia3D.RotationMatrix`](@ref) ` R` or a
[`Modia3D.Quaternion`](@ref) ` q` (to rotate a frame 1 into a frame 2).

| Function                                       | Description                                   |
|:-----------------------------------------------|:----------------------------------------------|
| [`Modia3D.resolve1`](@ref)(Rq, v2)             | Transform vector `v` from frame 2 to frame 1  |
| [`Modia3D.resolve2`](@ref)(Rq, v1)             | Transform vector `v` from frame 1 to frame 2  |
| [`Modia3D.absoluteRotation`](@ref)(Rq01, Rq12) | Return rotation 0->2 from rot. 0->1 and 1->2  |
| [`Modia3D.relativeRotation`](@ref)(Rq01, Rq02) | Return rotation 1->2 from rot. 0->1 and 0->2  |
| [`Modia3D.inverseRotation`](@ref)(Rq01)        | Return rotation 1->0 from rot, 0->1           |
| [`Modia3D.planarRotationAngle`](@ref)(e,v1,v2) | Return angle of planar rotation along `e`     |
| [`Modia3D.eAxis`](@ref)(axis)                  | Return unit vector `e` in direction of `axis` |
| [`Modia3D.skew`](@ref)(v)                      | Return skew-symmetric matrix of vector v      |

# Interpolation of Frames

Given a set of frames by a vector `r` of position vectors to their origins and
and an optional vector `q` of Quaternions of their absolute orientations, then
the following functions interpolate linearly in these frames:

| Function                                        | Description                                    |
|:------------------------------------------------|:-----------------------------------------------|
| [`Modia3D.Path`](@ref)(r,q)                     | Return path defined by a vector of frames      |
| [`Modia3D.t_pathEnd`](@ref)(path)               | Return path parameter `t_end` of last frame    |
| [`Modia3D.interpolate`](@ref)(path,t)           | Return `(rt,qt)` of Path at path parameter `t` |
| [`Modia3D.interpolate_r`](@ref)(path,t)         | Return `rt` of Path at path parameter `t`      |

# Examples

```julia
using Modia3D
using Unitful

# R1,R2,R3 are the same RotationMatrices
R1 = Modia3D.rot1(pi/2)
R2 = Modia3D.rot1(90u"°")
R3 = Modia3D.rot_e([1,0,0], pi/2)
```

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)

*The functions of this module are mostly a mapping of some of the functions
of the Modelica Standard Library from Modelica
([Modelica.Mechanics.MultiBody.Frames](https://doc.modelica.org/help/Modelica_Mechanics_MultiBody_Frames.html#Modelica.Mechanics.MultiBody.Frames))
to Julia (taking advantage of Julia features
such as multiple dispatch and unit package Unitful)*.
"""
module Frames

export RotationMatrix, NullRotation  , assertRotationMatrix
export Quaternion    , NullQuaternion, assertQuaternion
export Vector3D      , ZeroVector3D  , axisValue

export  rot1,  rot2,  rot3,  rot123, rotAxis,  rot_e,  rot_nxy, from_q
export qrot1, qrot2, qrot3, qrot123,          qrot_e, qrot_nxy, from_R

export resolve1, resolve2, absoluteRotation, relativeRotation, inverseRotation
export planarRotationAngle, eAxis

export Path, t_pathEnd, interpolate, interpolate_r
export skew

using StaticArrays
using LinearAlgebra

include("vector3D.jl")
include("rotationMatrix.jl")
include("quaternion.jl")
include("interpolation.jl")

end


text:
addqueen(queens::Array{Vector{Int}}, queen::Vector{Int}) = push!(copy(queens), queen)

hitsany(queen::Vector{Int}, queens::Array{Vector{Int}}) = any(map((x) -> hits(queen, x), queens))
hits(a::Array{Int}, b::Array{Int}) = any(a .== b) || abs(a-b)[1] == abs(a-b)[2]

function solve(x::Int, y::Int, n::Int, d::Array{Vector{Int}})
  if n == 0
    return d
  end
  for px = 1:x
    for py = 1:y
      if !hitsany([px, py], d)
        s = solve(x, y, n-1, addqueen(d, [px, py]))
        if !isequal(s, None)
          return s
        end
      end
    end
  end
  return None
end

solve(x, y, n) = solve(x, y, n, Array(Vector{Int}, 0))

for i = 1:8
  print("Solve for $i\n")
  print(solve(8, 8, i))
  print("\n")
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using BenchmarkTools
using LinearAlgebra

SUITE = BenchmarkGroup()
for n in [100, 400, 800]
    A = randn(n, n)
    B = randn(n, n)
    C = zero(A)
    SUITE["n=$n"] = @benchmarkable mul!($C, $A, $B)
end


text:
immutable LineVertex{N}
    position::Vec{N, Float32}
    thickness::Float32
    color::Vec4f0
end

immutable LineUniforms
    model::Mat4f0
    thickness::Float32
    pattern_length::Float32
end

immutable Vert2Geom
    position::Vec4f0
    color::Vec4f0
    thickness::Float32
end

immutable Geom2Fragment
    thickness::Float32
    color::Vec4f0
    uv::Vec2f0
end
to_vec4(v::Vec3f0) = return Vec4f0(v[1], v[2], v[3], 1f0)
to_vec4(v::Vec2f0) = return Vec4f0(v[1], v[2], 0f0, 1f0)

get_position(x::Vec4f0) = x
get_position(x::LineVertex) = x.position

get_color(x::LineVertex, uniforms) = x.color
get_color(x, uniforms) = uniforms.color

get_thickness(x::LineVertex, uniforms) = x.thickness
get_thickness(x, uniforms) = uniforms.thickness

function screen_space(vertex::Vec4f0, canvas)
    return (vertex[Vec(1, 2)] / vertex[4]) .* canvas.resolution
end

function vert_linesegments(vertex, canvas, uniforms)
    # m =  * uniforms.model
    geomout = Vert2Geom(
        canvas.projection * to_vec4(get_position(vertex)),
        get_color(vertex, uniforms),
        get_thickness(vertex, uniforms)
    )
    return geomout
end

function emit_vertex2(emit!, geom_in, canvas, uniforms, position::Vec2f0, uv::Vec2f0, index)
    inpos = geom_in[index].position
    x = (position ./ canvas.resolution) * inpos[4]
    outpos = Vec4f0(x[1], x[2], inpos[3], inpos[4])
    fragout = Geom2Fragment(
        geom_in[index].thickness + 2f0,
        geom_in[index].color,
        uv,
    )
    emit!(outpos, fragout)
    return
end

function geom_linesegments(emit!, vertex_out, canvas, uniforms)
    # get the four vertices passed to the shader:
    p0 = screen_space(vertex_out[1].position, canvas) # start of previous segment
    p1 = screen_space(vertex_out[2].position, canvas) # end of previous segment, start of current segment

    thickness0 = vertex_out[1].thickness
    thickness1 = vertex_out[2].thickness
    thickness_aa0 = thickness0 + 2f0
    thickness_aa1 = thickness1 + 2f0

    # determine the direction of each of the 3 segments (previous, current, next)
    vun0 = p1 - p0
    v0 = normalize(vun0)
    # determine the normal of each of the 3 segments (previous, current, next)
    n0 = Vec2f0(-v0[2], v0[1])
    l = norm(p1 - p0)

    uv0 = thickness_aa0 / thickness0
    uv1 = thickness_aa1 / thickness_aa1
    emit_vertex2(emit!, vertex_out, canvas, uniforms, p0 + thickness_aa0 * n0, Vec2f0(0f0, -uv0), 1)
    emit_vertex2(emit!, vertex_out, canvas, uniforms, p0 - thickness_aa0 * n0, Vec2f0(0f0, uv0), 1)
    emit_vertex2(emit!, vertex_out, canvas, uniforms, p1 + thickness_aa1 * n0, Vec2f0(l, -uv1), 2)
    emit_vertex2(emit!, vertex_out, canvas, uniforms, p1 - thickness_aa1 * n0, Vec2f0(l, uv1), 2)
    return
end

function frag_linesegments(geom_out, canvas, uniforms)
    uv = geom_out.uv; color = geom_out.color
    xy = Vec2f0(0.5f0, uv[2])
    alpha = aastep(0f0, xy[1])
    alpha2 = aastep(-1f0, 1f0, xy[2])
    outcolor = Vec4f0(color[1], color[2], color[3], color[4] * alpha * alpha2)
    (outcolor, )
end


text:
immutable Linear <: Degree{1} end

"""
Assuming uniform knots with spacing 1, the `i`th peice of linear b-spline
implemented here is defined as follows.

    y_i(x) = c p(x) + cp p(1-x)

where

    p(δx) = x

and the values `cX` for `X ∈ {_, p}` are the coefficients.

Linear b-splines are naturally interpolating, and require no prefiltering;
there is therefore no need for boundary conditions to be provided.

Also, although the implementation is slightly different in order to re-use
the framework built for general b-splines, the resulting interpolant is just
a piecewise linear function connecting each pair of neighboring data points.
"""
Linear

"""
`define_indices_d` for a linear b-spline calculates `ix_d = floor(x_d)` and
`fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring for
`Linear`), as well as the auxiliary quantity `ixp_d`
"""
function define_indices_d(::Type{BSpline{Linear}}, d, pad)
    symix, symixp, symfx, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("fx_",d), Symbol("x_",d)
    quote
        $symix = clamp(floor(Int, $symx), 1, size(itp, $d)-1)
        $symixp = $symix + 1
        $symfx = $symx - $symix
    end
end

"""
In `coefficients` for a linear b-spline we assume that `fx_d = x-ix-d` and
we define `cX_d` for `X ∈ {_, p}` such that

    c_d = p(fx_d)
    cp_d = p(1-fx_d)

where `p` is defined in the docstring entry for `Linear` and `fx_d` in the
docstring entry for `define_indices_d`.
"""
function coefficients(::Type{BSpline{Linear}}, N, d)
    sym, symp, symfx = Symbol("c_",d), Symbol("cp_",d), Symbol("fx_",d)
    quote
        $sym = 1 - $symfx
        $symp = $symfx
    end
end

"""
In `gradient_coefficients` for a linear b-spline we assume that `fx_d = x-ix_d`
and we define `cX_d` for `X ⋹ {_, p}` such that

    c_d  = p'(fx_d)
    cp_d = p'(1-fx_d)

where `p` is defined in the docstring entry for `Linear`, and `fx_d` in the
docstring entry for `define_indices_d`.
"""
function gradient_coefficients(::Type{BSpline{Linear}}, d)
    sym, symp, symfx = Symbol("c_",d), Symbol("cp_",d), Symbol("fx_",d)
    quote
        $sym = -1
        $symp = 1
    end
end

"""
In `hessian_coefficients` for a linear b-spline we assume that `fx_d = x-ix_d`
and we define `cX_d` for `X ⋹ {_, p}` such that

    c_d  = p''(fx_d)
    cp_d = p''(1-fx_d)

where `p` is defined in the docstring entry for `Linear`, and `fx_d` in the
docstring entry for `define_indices_d`. (These are both ≡ 0.)
"""
function hessian_coefficients(::Type{BSpline{Linear}}, d)
    sym, symp = Symbol("c_",d), Symbol("cp_",d)
    quote
        $sym = $symp = 0
    end
end

# This assumes fractional values 0 <= fx_d <= 1, integral values ix_d and ixp_d (typically ixp_d = ix_d+1,
#except at boundaries), and an array itp.coefs
function index_gen{IT<:DimSpec{BSpline}}(::Type{BSpline{Linear}}, ::Type{IT}, N::Integer, offsets...)
    if length(offsets) < N
        d = length(offsets)+1
        sym = Symbol("c_", d)
        symp = Symbol("cp_", d)
        return :($sym * $(index_gen(IT, N, offsets..., 0)) + $symp * $(index_gen(IT, N, offsets..., 1)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end


text:
import PAPI

info("PAPI Inner Product Test")
info("Using flops")

@printf("%12s %12s %12s %12s %12s %12s\n", "n", "ops", "2n", "difference", "% error", "mflops")

for n = 500:500:5000
    a = rand(1,n)
    x = rand(n, 1)

    ops, mflops = @PAPI.flops begin
        c = a * x
    end

    @printf("%12d %12d %12d %12d %12.2f %12.2f\n",
            n, ops, 2*n, ops-2n, (1.0 - (2n / ops))*100, mflops)
end


text:
# Utility function for classifier head of vision transformer-like models
_seconddimmean(x) = mean(x, dims = 2)[:, 1, :]

"""
    weights(model)

Load the pre-trained weights for `model` using the stored artifacts.
"""
function weights(model)
  try
    path = joinpath(@artifact_str(model), "$model.bson")
    return BSON.load(path, @__MODULE__)[:weights]
  catch e
    throw(ArgumentError("No pre-trained weights available for $model."))
  end
end

"""
    loadpretrain!(model, name)

Load the pre-trained weight artifacts matching `<name>.bson` into `model`.
"""
loadpretrain!(model, name) = Flux.loadparams!(model, weights(name))

function _maybe_big_show(io, model)
  if isdefined(Flux, :_big_show)
    if isnothing(get(io, :typeinfo, nothing)) # e.g. top level in REPL
      Flux._big_show(io, model)
    else
      show(io, model)
    end
  else
    show(io, model)
  end
end


text:
function copperplatebalance(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, time_periods::Int64)

    devices_netinjection = remove_undef!(netinjection.var_active)
    timeseries_netinjection = sum(netinjection.timeseries_active, dims=1)

    cpn = JuMP.JuMPArray(Array{ConstraintRef}(undef,time_periods), 1:time_periods)

    for t in 1:time_periods
        # TODO: Check is sum() is the best way to do this in terms of speed.
        cpn[t] = @constraint(m, sum(netinjection.var_active[:,t]) == timeseries_netinjection[t])
    end

    JuMP.register_object(m, :CopperPlateBalance, cpn)

    return m
end


text:
# compute an adjoint-based error estimate

"""
  Options structure for mesh adapation

  **Mesh Adaptation options**

   * strategy: strategy used by [`getTargetSizes`](@ref), default 1
   * free_mesh: if true, the old mesh object will be freed after adaptation
                is complete, Only set this to false if you know what you
                are doing.  Default true.
   * callback: a function that will be called by [`solveAdapative`](@ref)
               every for each new mesh.  The solution in `eqn.q` will be
               the solution to the PDE on that mesh.  The callback must have
               the signature:

               ```
                 function callback(mesh, sbp, eqn, opts, iter)
               ```

               where `iter` is the current iteration number

   * fixed_fraction: a number between 0 and 1 that determines what fraction
                     of the elements are refined for fixed-fraction based
                     strategy.  Default 0.1
  **SolveAdaptive Options**

   * solve_itermax: maximum number of times h-adaptation/solve PDE cycles
                    to do.  Default no limit
   * element_limit: maximum number of elements to generate.  The adaptive
               workflow will abort the first time the mesh has more
               elements than this, so the returned mesh will have
               (possibly many) more elements than this.  Default no
               limit.

  **Strategy**

  The values for the strategy are:

   * 1: attempt to equidistribute the error by refining or coarsening every
        element.  This strategy attempts to reach the error tolerance in a
        single iteration, which can lead to over-refinment
   * 2: refine a fixed fraction of the elements each iteration, reducing
        their size by half.  The elements with the greatest error are selected.
        This strategy does not consider the error tolerance.  Uses the
        `fixed_fraction` field.
"""
mutable struct AdaptOpts
  strategy::Int
  #TODO: wrap more MeshAdapt options
  free_mesh::Bool
  callback::Function  # callback for solveAdapative

  # options used by different strategies
  fixed_fraction::Float64

  # solveAdaptive options
  solve_itermax::Int
  el_limit::Int

  # private options
  free_ls::Bool
end

"""
  Returns object with default options
"""
function AdaptOpts()
  # mesh adaptation options
  strategy = 1
  free_mesh = true
  callback = defaultCallback

  fixed_fraction = 0.1
  #TODO: need to figure out about load balancing post adaptation

  # solveAdapative options
  solve_itermax = typemax(Int)
  element_limit = typemax(Int)

  # private options
  free_ls = false

  return AdaptOpts(strategy, free_mesh, callback, fixed_fraction, solve_itermax,
                   element_limit, free_ls)
end

"""
  Default (do-nothing) callback
"""
function defaultCallback(mesh::AbstractMesh, sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts, idx::Integer)

  # do nothing

  return nothing
end

"""
  Writes a visualization file with the error estimate.  The output file
  is called "error_estimate".

  **Inputs**

   * mesh
   * el_err: elementwise error estimate
"""
function writeErrorField(mesh, el_err::AbstractVector)

  fshape_ptr = apf.getConstantShapePtr(mesh.dim)
  f = apf.createPackedField(mesh, "error estimate", 1, fshape_ptr)

  val = zeros(Float64, 1)
  for i=1:mesh.numEl
    val[1] = el_err[i]
    apf.setComponents(f, mesh.elements[i], 0, val)
  end

  writeVisFiles(mesh, "error_estimate")

  apf.destroyField(mesh, f)

  return nothing
end

import NonlinearSolvers: NewtonBJacobiPC

"""
  Computes an adjoint-based error estimate and an estimate of each elements's
  contribution to the total error.

  **Inputs**
   * adapt_opts: an [`AdaptOpts`](@ref) object
   * mesh
   * sbp
   * eqn
   * opts
   * func: the functional
   * ls: a LinearSolver used to solve for the adjoint

  **Outputs**

   * err: the estimated error in the functional
   * el_error: the estimated contribution of each element to `err`

  **AdaptOpts Usage**

   * free_ls: if true, free the linear solver after solving the coarse
              space adjoint.  This should free the memory before any
              fine-space data structures are allocated

"""
function calcErrorEstimate(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                           sbp::AbstractOperator{Tsbp},
                           eqn::AbstractSolutionData{Tsol, Tres},
                           opts,
                           func::AbstractFunctional,
                           ls::LinearSolver) where {Tsol, Tres, Tsbp}

  # compute adjoint on coarse space
  psi_H = zeros(Tres, mesh.numDof)
  calcAdjoint(mesh, sbp, eqn, opts, ls, func, psi_H; recalc_jac=true,
              recalc_pc=true, start_comm=true)
  if adapt_opts.free_ls
    free(ls)
    #TODO: gc here to make sure the memory is really freed?
  end

  # construct objects on fine space
  mesh_h, sbp_h, eqn_h, opts_h = createEnrichedObjects(mesh, sbp, eqn, opts)

  # interpolate adjoint, solution to fine space
  psi_h = zeros(Tres, mesh_h.numDof)
  interpField(mesh, sbp, psi_H, mesh_h, sbp_h, psi_h)

  # start parallel communication
  setParallelData(eqn_h.shared_data, opts["parallel_data"])
  startSolutionExchange(mesh_h, sbp_h, eqn_h, opts_h)

  # smooth adjoint on fine space
  #TODO: verbose=false, res_tol < 0
  pc = NewtonBJacobiPC(mesh_h, sbp_h, eqn_h, opts_h, itermax=5, verbose=true, res_tol=1e-8)
  calcPC(pc, mesh_h, sbp_h, eqn_h, opts_h, (evalResidual,), 0.0)

  dJdu_h_arr = zeros(Tsol, mesh_h.numDofPerNode, mesh_h.numNodesPerElement,
                           mesh_h.numEl)
  dJdu_h = zeros(Tsol, mesh_h.numDof)
  evalFunctionalDeriv_q(mesh_h, sbp_h, eqn_h, opts_h, func, dJdu_h_arr,
                        start_comm=false)
  array3DTo1D(mesh_h, sbp_h, eqn_h, opts_h, dJdu_h_arr, dJdu_h)
  scale!(dJdu_h, -1)  # the adjoint problem is dR/dq^T psi = -dJ/du
  applyPCTranspose(pc, mesh_h, sbp_h, eqn_h, opts_h, 0.0, dJdu_h, psi_h)
  free(pc)

  # compute residual on fine space
  # parallel communication was alredy done by calcPC()
  evalResidual(mesh_h, sbp_h, eqn_h, opts_h)
  array3DTo1D(mesh_h, sbp_h ,eqn_h, opts_h, eqn_h.res, eqn_h.res_vec)

  err, el_error = localizeErrorEstimate(mesh_h, sbp_h, eqn_h, opts_h,
                                        eqn_h.res_vec, psi_h)

  finalize(mesh_h)

  if opts["write_error_estimate"]
    writeErrorField(mesh, el_error)
  end

  return err, el_error

end

"""
  Computes the contribution of each element to the total error

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * R_h: fine grid residual evaluated using the interpolated solution
   * psi_h: fine grid adjoint

  **Outputs**

   * err: the estimate of the total error
   * el_error: the estimate of each element's contribution to the total error
"""
function localizeErrorEstimate(mesh::AbstractMesh,
                              sbp::AbstractOperator,
                              eqn::AbstractSolutionData{Tsol, Tres},
                              opts,
                              R_h::AbstractVector,
                              psi_h::AbstractVector) where {Tsol, Tres}
# the arguments are the fine-space objects
  # compute \sum | psi_j R_j |
  err = zero(Float64)
  @simd for i=1:length(R_h)
    err += real(R_h[i]*psi_h[i])
  end
  err = MPI.Allreduce(abs(err), MPI.SUM, eqn.comm)

  el_error = zeros(Float64, mesh.numEl)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        dof = mesh.dofs[k, j, i]
        el_error[i] += abs(R_h[dof]*psi_h[dof])
      end
    end
  end

  return err, el_error
end

"""
  Given the elemnt-local error, this function computes a new size field.

  **Inputs**

   * adapt_opts: an [`AdaptOpts`](@ref) object
   * mesh
   * sbp
   * eqn
   * opts
   * el_error: vector containing the element contribution to the total
               error, length `mesh.numEl`
   * err_target: the target total error (not per-element)

  **Outputs**

   * el_sizes: vector containing the target size of every element.

  **AdaptOpts Usage**

   * strategy: integer indicating what strategy to use for computing
               the sizes, currently:
     * 1: refine everything using theoretical convergence rate of p + 1
"""
function getTargetSizes(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                         sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts,
                         el_error::AbstractVector,
                         err_target::Number)

  el_sizes = getElementSizes(mesh)

  # currently we refine everything to equi-distribute the error.  This
  # likely leads to over-refinement

  err_target_el = err_target/mesh.numEl
  rate::Int = convert(Int, opts["order"] + 1)  # theoretical convergence rate

  if adapt_opts.strategy == 1
    getTargetSizes_1(adapt_opts, mesh, sbp, eqn, opts, el_error, err_target,
                   el_sizes)
  elseif adapt_opts.strategy == 2
    getTargetSizes_2(adapt_opts, mesh, sbp, eqn, opts, el_error, err_target,
                     el_sizes)
  else
    error("unrecognized sizing strategy: $strategy")
  end

  return el_sizes
end

import PdePumiInterface.adaptMesh

"""
  Wrapper for taking a new size field and adapting the mesh to it.
  Constructs new `AbstractSolutionData` object

  **Inputs**

   * adapt_opts: mesh adaptation options
   * mesh
   * sbp
   * eqn
   * opts
   * el_sizes: vector of length `mesh.numEl` containing the new size for
               each element.

  **Outputs**

   * newmesh
   * newsbp
   * neweqn
   * newopts

"""
function adaptMesh(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                   sbp::AbstractOperator,
                   eqn::AbstractSolutionData,
                   opts,
                   el_sizes::AbstractVector)

  #TODO: need options for mesh adaptation
  #      The main thing we want is load balancing

  # perform mesh adaptation (call PumiInterface function)
  newmesh, q_new = adaptMesh(mesh, sbp, opts, el_sizes, eqn.q_vec,
                             free_mesh=adapt_opts.free_mesh)

  # construct new objects
  newopts = deepcopy(opts)
  newmesh, newsbp, neweqn, newopts, pmesh = createObjects(newmesh, sbp,
                                                          newopts)
  copy!(neweqn.q_vec, q_new)

  # return new objects
  return newmesh, newsbp, neweqn, newopts

end

#------------------------------------------------------------------------------
# Interface functions

"""
  Computes an error estimate for the given functional.

  The PDE should have been solved before calling this function.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * func: the `AbstractFunctional`

  **Outputs**

   * err: the estimated functional error
   * el_error: the estimated contribution of each element to `err`
"""
function calcErrorEstimate(mesh::AbstractMesh,
                       sbp::AbstractOperator,
                       eqn::AbstractSolutionData,
                       opts,
                       func::AbstractFunctional)

  # make Newton linear solver
  adapt_opts = AdaptOpts()

  pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  adapt_opts.free_ls = true

  # call other method
  return calcErrorEstimate(adapt_opts, mesh, sbp, eqn, opts, func, ls)
end

"""
  Does one round of h-adaptation, using an adjoint-based error estimate.
  If the returned objects are used to
  solve the PDE, the functional error should be smaller, however it
  may not meet the supplied tolerance.  To continue doing h-adaptation until
  the tolerance is met, see TODO: other function

  On entry, the solution in `eqn.q_vec` should be the solution to the PDE.

  **Inputs**

   * adapt_opts: an [`AdaptOpts`](@ref) object
   * mesh
   * sbp
   * eqn
   * opts
   * func: functional to use for computing adjoint-based error estimate
   * ls: a LinearSolver object.  The linear operator and preconditioner
         will be re-calculated before the linear solve is done
   * err_target: target functional error

  **Outputs**

   * newmesh: new mesh object
   * newsbp: new sbp object
   * neweqn: new eqn object
   * newopts: new options dictonary
   * err: the estimated error at the *prior* to adaptation.  If this
          value is less than the `err_target`, then no adaptation is done
          and the `err` is the estimated error on the current mesh.

  Note that if the original mesh satisfies the error tolerance, the
  input objects are returned and no mesh adaptation is done.

  **Keyword Arguments**

   * free_ls: if true, the linear solver will be freed after use.  This
              can decrease the maximum memory usage of the algorithm.
"""
function doHAdaptation(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                       sbp::AbstractOperator,
                       eqn::AbstractSolutionData,
                       opts,
                       func::AbstractFunctional,
                       ls::LinearSolver,
                       err_target::Number;
                       free_ls::Bool=false)
  # do all the steps
  adapt_opts.free_ls = adapt_opts.free_ls || free_ls

  err, el_err = calcErrorEstimate(adapt_opts, mesh, sbp, eqn, opts, func, ls)

  if err > err_target  #TODO: some tolerance here?
    el_sizes = getTargetSizes(adapt_opts, mesh, sbp, eqn, opts, el_err,
                              err_target)
    newmesh, newsbp, neweqn, newopts = adaptMesh(adapt_opts, mesh, sbp, eqn,
                                                 opts, el_sizes)
  else
    newmesh = mesh
    newsbp = sbp
    neweqn = eqn
    newopts = opts
  end

  return newmesh, newsbp, neweqn, newopts, err
end

"""
  This method constructs its own linear solver, rather than using a provided
  one.  All other arguments are the same, as well as the outputs.

  **Inputs**

   * adapt_opts
   * mesh
   * sbp
   * eqn
   * opts
   * func
   * err_target
"""
function doHAdaptation(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                       sbp::AbstractOperator,
                       eqn::AbstractSolutionData,
                       opts,
                       func::AbstractFunctional,
                       err_target::Number)
  # make Newton linear solver
  pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

  # call other method
  return doHAdaptation(adapt_opts, mesh, sbp, eqn, opts, func, ls,
                       err_target, free_ls=true)

end

"""
  This function performs repeated h-adaptation until the error estimate
  meets the tolerance (or some other criteria causes adaptation to stop
  early.

  **Inputs**

   * adapt_opts: an [`AdaptOpts`](@ref).  In particular, see the
                 `SolveAdapative` section for options that affects this
                 function
   * mesh
   * sbp
   * eqn
   * opts
   * func: the [`AbstractFunctional`](@ref) to estimate the error of
   * err_target: the target error

  **Outputs**

   * mesh: the adapted mesh
   * sbp
   * eqn: the equation object for the adapted mesh
   * opts; the options dictionary.

"""
function solveAdaptive(adapt_opts::AdaptOpts, mesh::AbstractMesh,
                       sbp::AbstractOperator, eqn::AbstractSolutionData,
                       opts,
                       func::AbstractFunctional,
                       err_target::Number)
  # we assume the PDE was *not* previously solved, unlike doHAdaptation
  myrank = mesh.myrank

  # keep solving until tolerance met or some other criteria
  ic_orig = opts["IC_name"]
  exit_code = 1  # 1 = itermax, 0 = anything else
  for i=1:adapt_opts.solve_itermax

    solvePDE(mesh, sbp, eqn, opts)
    # for all future solves, use the previous solution
    opts["IC_name"] = "ICPassThrough"

    if opts["write_adapt_vis"]
      saveSolutionToMesh(mesh, eqn.q_vec)
      writeVisFiles(mesh, "adapt_$i")
    end

    # call the user-defined callback
    adapt_opts.callback(mesh, sbp, eqn, opts, i)

    # do another round of h-adaptation (or return the current objects if the
    # error meets the tolerance)
    old_numEl = mesh.numEl
    h_old = calcMeshH(mesh, sbp, eqn, opts)
    mesh, sbp, eqn, opts, err = doHAdaptation(adapt_opts, mesh, sbp, eqn,
                                              opts, func, err_target)

    # printing and convergence checks
    println("\n\nOn iteration ", i, ", error estimate = ", err, ", numEl = ", old_numEl, ", avg mesh size = ", h_old)

    numel_global = MPI.Allreduce(mesh.numEl, MPI.SUM, eqn.comm)
    if err < err_target
      @mpi_master println(BSTDOUT, "solveAdaptive met error tolerance with final error estimate of ", err, " < ", err_target, " with ", numel_global, " elements")
      exit_code = 0
      break
    end

    if numel_global > adapt_opts.el_limit
      @mpi_master println(BSTDOUT, "solveAdaptive exiting due to element limit.  Number of elements = ", numel_global, " > ", adapt_opts.element_limit)
      exit_code = 0
      break
    end

  end  # end loop i

  if exit_code == 1
    @mpi_master println(BSTDOUT, "solveAdaptive exited due to itermax, number of elements = ", mesh.numGlobalEl)
  end
  flush(BSTDOUT)

  opts["IC_name"] = ic_orig  # restore the initial condition
  return mesh, sbp, eqn, opts
end


text:
abstract type Absorber end

struct AbsorberProperties{T <: Real}
    μ_a :: T
    μ_s :: T
    g   :: T
    nt  :: T
    albedo :: T
end
AbsorberProperties(μ_a, μ_s, g, nt) = AbsorberProperties(μ_a, μ_s, g, nt, μ_s/(μ_s + μ_a))

struct CatesianVolume{T <: Real} <: Absorber
    corner1 :: SVector{3,T}
    corner2 :: SVector{3,T}
    d_xyz :: SVector{3,T}
    Δxyz :: SVector{3,T}
    Nxyz :: SVector{3,Integer}
    prop :: AbsorberProperties
    # data :: SizedArray{Tuple, T}
    data
    function CatesianVolume(corner1::Vector{T}, corner2::Vector{T}, d_xyz::Vector{T}, prop::Vector{T}) where {T<:Real}
        corner1 = SVector{3,T}(corner1)
        corner2 = SVector{3,T}(corner2)
        Δxyz = SVector{3,T}(corner2 - corner1)
        Nxyz = SVector{3,Integer}(Int.(round.(Δxyz ./ d_xyz)))
        prop = AbsorberProperties(prop[1], prop[2], prop[3], prop[4])
        data = SizedArray{Tuple{Nxyz...}, T}(zeros(Nxyz...))
        new{T}(corner1, corner2, d_xyz, Δxyz, Nxyz, prop, data)
    end
end
CatesianVolume(lower::T, upper::T, dxyz::T, prop::Vector) where T <: Real = CatesianVolume(
    lower*ones(3),
    upper*ones(3),
    dxyz* ones(3),
    prop
)
CatesianVolume(lower::Real, upper::Real, dxyz::Real, prop::Vector) = CatesianVolume(
    promote(lower*ones(3),  upper*ones(3),  dxyz* ones(3),  prop)...
)

struct RadialVolume{T <: Real} <: Absorber
    r1 :: T
    r2 :: T
    dr :: T
    Δr :: T
    Nr :: Integer
    data :: Vector{T}
    function RadialVolume(r1::T, r2::T, dr::T) where {T<:Real}
        Δr = r2 - r1
        Nr = Int(round(Δr / dr))
        data = zeros(Nr)
        new{T}(r1, r2, dr, Δr, Nr, data)
    end
end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d

include("init_pot_Hcoul.jl")
include("init_pot_harm_3d.jl")
include("../utils/ortho_gram_schmidt.jl")
include("../utils/orthonormalize.jl")
include("EnergiesT.jl")
include("op_nabla2.jl")
include("op_H.jl")
include("calc_Energies.jl")
include("calc_grad.jl")

function test_grad()
  # LF parameters
  NN = [25, 25, 25]
  AA = [0.0, 0.0, 0.0]
  BB = [6.0, 6.0, 6.0]

  Npoints = prod(NN)

  # Initialize LF
  LF = init_LF3d_p( NN, AA, BB )
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # Parameter for potential
  center = AA + 0.5*(BB-AA)
  # Potential
  ω = 2.0
  Vpot = init_pot_harm_3d( LF, ω, center )

  Ncols = 4

  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  println("sum(Vpot) = ", sum(Vpot))
  println("ΔV = ", ΔV)

  v = zeros(Npoints,Ncols)
  for i = 1:Ncols
    v[i,i] = 1.0/sqrt(ΔV)
  end

  grad = calc_grad( LF, Vpot, v )
  println(sum(grad))

  Energies = calc_Energies( LF, Vpot, v)
  println(Energies)
end

test_grad()