text:
export CompositionSystem, compose

"""
    CompositionSystem(G::AbstractSystem, F::AbstractSystem)

Construct the system ``G(F(x;p);p)``. Note that the parameters are passed to ``G`` and
``F`` are identical.
"""
struct CompositionSystem{S1<:AbstractSystem,S2<:AbstractSystem} <: AbstractSystem
    # The system is g ∘ f
    g::S2
    f::S1

    f_u::Vector{ComplexF64}
    f_ū::Vector{ComplexDF64}
    f_U::Matrix{ComplexF64}
    g_U::Matrix{ComplexF64}
    tu⁴::TaylorVector{5,ComplexF64}
    tu³::TaylorVector{4,ComplexF64}
    tu²::TaylorVector{3,ComplexF64}
    tu¹::TaylorVector{2,ComplexF64}
end
function CompositionSystem(g::AbstractSystem, f::AbstractSystem)
    size(g, 2) == size(f, 1) ||
        throw(ArgumentError("Cannot create composition G ∘ F since the number of variabels of `G` and the number polynomials of `F` doesn't match."))
    f_u = zeros(ComplexF64, size(f, 1))
    f_ū = zeros(ComplexDF64, size(f, 1))
    f_U = zeros(ComplexF64, size(f))
    g_U = zeros(ComplexF64, size(g))

    tu⁴ = TaylorVector{5}(ComplexF64, size(f, 1))
    tu¹ = TaylorVector{2}(tu⁴)
    tu² = TaylorVector{3}(tu⁴)
    tu³ = TaylorVector{4}(tu⁴)

    CompositionSystem(g, f, f_u, f_ū, f_U, g_U, tu⁴, tu³, tu², tu¹)
end

Base.size(C::CompositionSystem) = (size(C.g, 1), size(C.f, 2))

ModelKit.variables(F::CompositionSystem) = variables(F.f)
ModelKit.parameters(F::CompositionSystem) = parameters(F.f)
ModelKit.variable_groups(F::CompositionSystem) = variable_groups(F.f)

function Base.show(io::IO, C::CompositionSystem)
    println(io, "Composition G ∘ F:")
    println(io, "F: ")
    show(io, C.f)
    println(io, "\nG: ")
    show(io, C.g)
end

(C::CompositionSystem)(x, p = nothing) = C.g(C.f(x, p), p)
function evaluate!(u, C::CompositionSystem, x::AbstractVector{ComplexF64}, p = nothing)
    evaluate!(u, C.g, evaluate!(C.f_u, C.f, x, p), p)
end
function evaluate!(u, C::CompositionSystem, x::AbstractVector{ComplexDF64}, p = nothing)
    evaluate!(u, C.g, evaluate!(C.f_ū, C.f, x, p), p)
end

function evaluate_and_jacobian!(u, U, C::CompositionSystem, x, p = nothing)
    evaluate_and_jacobian!(C.f_u, C.f_U, C.f, x, p)
    evaluate_and_jacobian!(u, C.g_U, C.g, C.f_u, p)
    LA.mul!(U, C.g_U, C.f_U)
    nothing
end

_get_tu(C::CompositionSystem, ::Val{1}) = C.tu¹
_get_tu(C::CompositionSystem, ::Val{2}) = C.tu²
_get_tu(C::CompositionSystem, ::Val{3}) = C.tu³
_get_tu(C::CompositionSystem, ::Val{4}) = C.tu⁴

function taylor!(u, v::Val, C::CompositionSystem, tx, p = nothing)
    tu = _get_tu(C, v)
    taylor!(tu.data, v, C.f, tx, p)
    taylor!(u, v, C.g, tu, p)
end

"""
    compose(G::Union{AbstractSystem,System}, F::Union{AbstractSystem,System})

Construct the composition ``G(F(x))``. You can also use the infix operator
`∘` (written by \\circ).

### Example
```julia
julia> @var a b c x y z

julia> g = System([a * b * c]);

julia> f = System([x+y, y + z, x + z]);

julia> compose(g, f)
Composition G ∘ F:
F:
ModelKitSystem{(0xbb16b481c0808501, 1)}:
Compiled: System of length 3
 3 variables: x, y, z

 x + y
 y + z
 x + z

G:
ModelKitSystem{(0xf0a2384a42428501, 1)}:
Compiled: System of length 1
 3 variables: a, b, c

 a*b*c

julia> (g ∘ f)([x,y,z])
1-element Array{Expression,1}:
 (x + z)*(y + z)*(x + y)
```
"""
compose(g::AbstractSystem, f::AbstractSystem) = CompositionSystem(g, f)
compose(g::AbstractSystem, f::System) = CompositionSystem(g, ModelKitSystem(f))
function compose(g::System, f::System)
    pg = parameters(g)
    pf = parameters(f)
    (isempty(pf) || isempty(pg) || pf == pg) ||
        throw(ArgumentError("Cannot construct a composition of two system with different sets of parameters."))
    CompositionSystem(ModelKitSystem(g), ModelKitSystem(f))
end
import Base: ∘
∘(g::Union{AbstractSystem,System}, f::Union{AbstractSystem,System}) = compose(g, f)


text:
using OptiMode
using Test

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

ext: jl


text:
##===----------------------------------------------------------------------===##
##                                   ACCELERATE.JL                            ##
##                                                                            ##
##===----------------------------------------------------------------------===##
##                                                                            ##
##  This file provides access to Apple's Accelerate on OS X systems           ##
##                                                                            ##
##===----------------------------------------------------------------------===##
module Accelerate
import Vectorize: functions, addfunction

const libacc = "/System/Library/Frameworks/Accelerate.framework/Accelerate"

## This defines the mapping between Julia function names and Accelerate vecLib names.
## Each tuple is of the form (:julia_name, :accelerate_name).accelerate names should
## be stored without the "vv" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float64, and "f" suffix for Float32.
##
## This assumes arguments of the form (output_vector, input_vector, length)
##
const veclibfunctions =
    ((:ceil, :ceil), (:floor, :floor), (:sqrt, :sqrt), (:invsqrt, :rsqrt), (:inv, :rec),
     (:exp, :exp), (:exp2, :exp2), (:expm1, :expm1), (:log, :log), (:log1p, :log1p),
     (:log2, :log2), (:log10, :log10), (:sin, :sin), (:sinpi, :sinpi), (:cos, :cos),
     (:cospi, :cospi), (:tan, :tan), (:tanpi, :tanpi), (:asin, :asin), (:acos, :acos),
     (:atan, :atan), (:sinh, :sinh), (:cosh, :cosh), (:tanh, :tanh), (:asinh, :asinh),
     (:acosh, :acosh), (:atanh, :atanh), (:trunc,:int), (:round,:nint), (:exponent,:logb),
     (:abs,:fabs))

## For vecLibFunctions
for (T, suff) in ((Float64, ""), (Float32, "f"))

    for (f, fa) in veclibfunctions
        name = string(f)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T},)), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)})`
            Implements element-wise **$($name)** over a **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """
            function ($f)(X::Array{$T})
                out = similar(X)
                return ($f!)(out, X)
            end
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)})`
            Implements element-wise **$($name)** over a **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """
            function ($f!)(out::Array{$T}, X::Array{$T})
                ccall(($(string("vv", fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{Cint}),
                      out,X,Ref{Cint}(length(X)))
                return out
            end
        end
    end
    # 2 arg functions
    for f in (:copysign,)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T}, Array{T})), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T}, Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = similar(X)
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,Y,Ref{Cint}(length(X)))
                out
            end
        end
    end

    # for some bizarre reason, vvpow/vvpowf reverse the order of arguments.
    for f in (:pow,)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T}, Array{T})), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T}, Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = similar(X)
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,Y,X,Ref{Cint}(length(X)))
                out
            end
        end
    end

    # renamed 2 arg functions
    for (f,fa) in ((:rem,:fmod),(:fdiv,:div),(:atan,:atan2))
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T}, Array{T})), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T}, Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = similar(X)
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,Y,Ref{Cint}(length(X)))
                out
            end
        end
    end

    # two-arg return
    for f in (:sincos,)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T},)), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T}, Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            function ($f)(X::Array{$T})
                out1 = similar(X)
                out2 = similar(X)
                ($f!)(out1, out2, X)
            end
            function ($f!)(out1::Array{$T}, out2::Array{$T}, X::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out1,out2,X,Ref{Cint}(length(X)))
                out1, out2
            end
        end
    end

    # complex return
    for (f,fa) in ((:cis,:cosisin),)
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T},)), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{Complex{T}},Array{T})), "Vectorize.Accelerate.$(f!)")
        @eval begin
            function ($f)(X::Array{$T})
                out = Array{Complex{$T}}(undef, size(X))
                ($f!)(out, X)
            end
            function ($f!)(out::Array{Complex{$T}}, X::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{Complex{$T}},Ptr{$T},Ptr{Cint}),out,X,Ref{Cint}(length(X)))
                out
            end
        end
    end

end

## This defines the mapping between Julia function names and Accelerate vDSP names.
## Each tuple is of the form (:julia_name, :accelerate_name); accelerate names should
## be stored without the "vDSP" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float32, and "D" suffix for Float64.
##
## This assumes arguments of the form (input_vec_a, stride_a, input_vec_b, stride_b,
##                                                   output_vec, stride_out, length)
##
const vDSPfunctions = ((:add, :vadd, "addition"), (:sub, :vsub, "subtraction"),
                       (:div, :vdiv, "division"), (:mul, :vmul, "multiplication"),
                       (:max, :vmax, "maximum"), (:min, :vmin, "minimum"))

for (T, suff) in ((Float64, "D"), (Float32, ""))

    for (f, fa, name) in vDSPfunctions
        f! = Symbol("$(f)!")
        addfunction(functions, (f, (Array{T}, Array{T})), "Vectorize.Accelerate.$f")
        addfunction(functions, (f!, (Array{T}, Array{T}, Array{T})), "Vectorize.Accelerate.$f!")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`
            Implements element-wise **$($name)** over two **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """
            function ($f)(X::Array{$T}, Y::Array{$T})
                out = Array{$T}(undef, length(X))
                return ($f!)(out, Y, X)
            end
        end
        @eval begin
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)}, Y::Vector{$($T)})`
            Implements element-wise **$($name)** over two **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vDSP_", fa,suff)),libacc), Cvoid,
                      (Ptr{$T},Int64, Ptr{$T}, Int64, Ptr{$T}, Int64, Int64),
                      X, 1, Y, 1, out, 1, length(out))
                return out
            end
        end
    end
end

## This defines the mapping between Julia function names and Accelerate vDSP names that
## return a scalar value.
## Each tuple is of the form (:julia_name, :accelerate_name); accelerate names should
## be stored without the "vDSP" prefix. The accelerate names are stored without type
## suffix; it is assumed that there is no suffix for Float32, and "D" suffix for Float64.
##
## This assumes arguments of the form (input_vec_a, stride_a,  output_scalar, length)
const vDSPscalar = ((:sum,  :sve, "sum"), (:summag, :svemg, "sum-of-magnitudes"),
                    (:sumsqr, :svesq, "sum-of-squares"), (:maximum, :maxv, "maximum value"),
                    (:minimum, :minv, "minimum value"),  (:mean, :meanv, "mean value"))

for (T, suff) in ((Float64, "D"), (Float32, ""))

    for (f, fa, name) in vDSPscalar
        addfunction(functions, (f, (Array{T},)), "Vectorize.Accelerate.$f")
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`
            Computes the **$($name)** of a **Vector{$($T)}**.
            *Returns:* **$($T)**
            """
            function ($f)(X::Array{$T})
                out = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Cint, Ref{$T}, Cint),
                      X, 1, out, length(X))
                return out[]
            end
        end
    end

    for (f, fa, name) in ((:findmax, :maxvi, "maximum-value with index"),
                          (:findmin, :minvi, "minimum-value with index"))
        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)})`
            Computes the **$($name)** element-wise over a **Vector{$($T)}**.
            *Returns:* **$(($T, UInt64))**
            """
            function ($f)(X::Array{$T})
                index = Ref{Int}(0)
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Cvoid,
                      (Ptr{$T}, Int64,  Ref{$T}, Ref{Int}, UInt64),
                      X, 1, val, index, length(X))
                return (val[], index[]+1)
            end
        end
    end
end

end # module


text:
module Refill

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Refill" => Ahorn.EntityPlacement(
        Maple.Refill
    ),

    "Refill (Two Dashes)" => Ahorn.EntityPlacement(
        Maple.Refill,
        "point",
        Dict{String, Any}(
            "twoDash" => true
        )
    )
)

spriteOneDash = "objects/refill/idle00"
spriteTwoDash = "objects/refillTwo/idle00"

function getSprite(entity::Maple.Refill)
    twoDash = get(entity.data, "twoDash", false)

    return twoDash ? spriteTwoDash : spriteOneDash
end

function Ahorn.selection(entity::Maple.Refill)
    x, y = Ahorn.position(entity)
    sprite = getSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Refill, room::Maple.Room)
    sprite = getSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
export @energy

energy_dict = Dict{String, String}(
    "rhf" => "Fermi.HartreeFock.RHF()",
    "uhf" => "Fermi.HartreeFock.UHF()",
    "mp2" => "Fermi.MollerPlesset.RMP2()",
    "rmp2" => "Fermi.MollerPlesset.RMP2()",
    "ccsd" => "Fermi.CoupledCluster.RCCSD()",
    "ccsd(t)"=> "Fermi.CoupledCluster.RCCSDpT()",
)

"""
    Fermi.@energy

Macro to call functions to compute energy given current options. Arguments may be passed
using "=>" or "<="

# Examples

Generating a RHF wave function
```
rhf = @energy rhf
```

Using this wave function in a CCSD computations
```
@energy rhf => CCSD
```
or
```
@energy CCSD <= rhf
```

# Implemented methods:
    RHF            Restricted Hartree-Fock.
    UHF            Unrestricted Hartree-Fock
    RMP2           Restricted Moller-Plesset PT order 2.
    RMP3           Restricted Moller-Plesset PT order 3.
    CASCI          Restricted configuration interaction within given active space.
    CCSD           Restricted Coupled-Cluster with Single and Double substitutions
    CCSD(T)        Restricted Coupled-Cluster with Single and Double substitutions and perturbative triples.
    ecCCSD         Restricted externally corrected CCSD. External correction from CASCI computation.
    ecCCSD(T)      Restricted externally corrected CCSDpT. External correction from CASCI computation.
"""
macro energy(comm)
    # Clean spaces and colon from string
    clean_up(s) = String(filter(c->!occursin(c," :"),s))

    # Transform the input expression into a string
    A = clean_up(repr(comm))

    # If the command string is in between parenthesis, remove it
    while A[1] == '(' && A[end] == ')'
        A = A[2:end-1]
    end

    # If "=>" is present, then the second part is taken as the method and the first as the argument passed
    if occursin("=>", A)
        arg, method = split(A, "=>")
    elseif occursin("<=", A)
        method, arg = split(A, "<=")
    else
        method = A
        arg = ""
    end

    method = lowercase(String(method))
    arg = String(arg)   # Make sure arg is not a substring

    out = ""
    try
         out = replace(energy_dict[method], "()" => "("*arg*")")
    catch KeyError
        throw(FermiException("Invalid method for energy computation: \"$A\""))
    end

    expr_out = Meta.parse(out)

    return esc(expr_out)
end


text:
#=
Cohen:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

##A is the input matrix, Output: Triangle Count
function cohen(A)

    #extracting lower matrix from input matrix
    L = tril(A)

    #extracting upper matrix from input matrix
    U = triu(A')

    B = mul(L, U, Semirings.PLUS_TIMES)
    C = emul(B, A,  BinaryOps.TIMES)
    result = reduce(Monoids.PLUS_MONOID[Float64], C, dims=:)
    result = result /2
    print("RESULT ")
    print(result)

    return result

end


text:
using Mocking

create_file_system_1_patch = @patch function create_file_system(a...; b...)
    return Dict{String, Any}("PerformanceMode" => "generalPurpose", "CreationTime" => "1481841524.0", "NumberOfMountTargets" => 0, "SizeInBytes" => Dict{String, Any}("Value" => 0), "LifeCycleState" => "creating", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "CreationToken" => "tokenstring", "Encrypted" => true, "FileSystemId" => "fs-01234567", "OwnerId" => "012345678912")
end

create_mount_target_1_patch = @patch function create_mount_target(a...; b...)
    return Dict{String, Any}("MountTargetId" => "fsmt-12340abc", "IpAddress" => "192.0.0.2", "LifeCycleState" => "creating", "NetworkInterfaceId" => "eni-cedf6789", "FileSystemId" => "fs-01234567", "OwnerId" => "012345678912", "SubnetId" => "subnet-1234abcd")
end

describe_file_systems_1_patch = @patch function describe_file_systems(a...; b...)
    return Dict{String, Any}("FileSystems" => Any[Dict{String, Any}("PerformanceMode" => "generalPurpose", "CreationTime" => "1481841524.0", "NumberOfMountTargets" => 1, "SizeInBytes" => Dict{String, Any}("Value" => 6144), "LifeCycleState" => "available", "Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")], "CreationToken" => "tokenstring", "FileSystemId" => "fs-01234567", "OwnerId" => "012345678912", "Name" => "MyFileSystem")])
end

describe_lifecycle_configuration_1_patch = @patch function describe_lifecycle_configuration(a...; b...)
    return Dict{String, Any}("LifecyclePolicies" => Any[Dict{String, Any}("TransitionToIA" => "AFTER_30_DAYS")])
end

describe_mount_target_security_groups_1_patch = @patch function describe_mount_target_security_groups(a...; b...)
    return Dict{String, Any}("SecurityGroups" => Any["sg-4567abcd"])
end

describe_mount_targets_1_patch = @patch function describe_mount_targets(a...; b...)
    return Dict{String, Any}("MountTargets" => Any[Dict{String, Any}("MountTargetId" => "fsmt-12340abc", "IpAddress" => "192.0.0.2", "LifeCycleState" => "available", "NetworkInterfaceId" => "eni-cedf6789", "FileSystemId" => "fs-01234567", "OwnerId" => "012345678912", "SubnetId" => "subnet-1234abcd")])
end

describe_tags_1_patch = @patch function describe_tags(a...; b...)
    return Dict{String, Any}("Tags" => Any[Dict{String, Any}("Value" => "MyFileSystem", "Key" => "Name")])
end

put_lifecycle_configuration_1_patch = @patch function put_lifecycle_configuration(a...; b...)
    return Dict{String, Any}("LifecyclePolicies" => Any[Dict{String, Any}("TransitionToIA" => "AFTER_30_DAYS")])
end


text:
using Plots

f(x) = (1+0.5*x)*sin(x)

xs = range(3, 12, length=1000)
fs = f.(xs)
i_min = findall((fs[2:end-1] .< fs[1:end-2]) .& (fs[2:end-1] .< fs[3:end])) .- 1
i_max = findall((fs[2:end-1] .> fs[1:end-2]) .& (fs[2:end-1] .> fs[3:end])) .- 1

plot(xs, fs, label="", ticks=false)
scatter!(xs[i_min], fs[i_min], label="Local minimum")
scatter!(xs[i_max], fs[i_max], label="Local maximum")

file_name = joinpath("src", "lecture_07", "minmax.svg")
savefig(file_name)

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
# CodeTable64
# ===========

const CodeTable64 = CodeTable{64}

const BASE64_CODEPAD = 0x40  # PADding
const BASE64_CODEIGN = 0x41  # IGNore
const BASE64_CODEEND = 0x42  # END
const BASE64_CODEERR = 0xff  # ERRor

ignorecode(::Type{CodeTable64}) = BASE64_CODEIGN

"""
    CodeTable64(asciicode::String, pad::Char)

Create a code table for base64.
"""
function CodeTable64(asciicode::String, pad::Char)
    if !isascii(asciicode) || !isascii(pad)
        throw(ArgumentError("the code table must be ASCII"))
    elseif length(asciicode) != 64
        throw(ArgumentError("the code size must be 64"))
    end
    encodeword = Vector{UInt8}(undef, 64)
    decodeword = Vector{UInt8}(undef, 256)
    fill!(decodeword, BASE64_CODEERR)
    for (i, char) in enumerate(asciicode)
        bits = UInt8(i-1)
        code = UInt8(char)
        encodeword[bits+1] = code
        decodeword[code+1] = bits
    end
    padcode = UInt8(pad)
    decodeword[padcode+1] = BASE64_CODEPAD
    return CodeTable64(encodeword, decodeword, padcode)
end

@inline function encode(table::CodeTable64, byte::UInt8)
    return table.encodeword[Int(byte & 0x3f) + 1]
end

@inline function decode(table::CodeTable64, byte::UInt8)
    return table.decodeword[Int(byte)+1]
end

"""
The standard base64 code table (cf. Table 1 of RFC4648).
"""
const BASE64_STD = CodeTable64(
    "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/", '=')

"""
The url-safe base64 code table (cf. Table 2 of RFC4648).
"""
const BASE64_URLSAFE = CodeTable64(
    "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-_", '=')


text:
"""
    module Integration

Methods for integrating over instances of [`AbstractReferenceShape`](@ref).

Besides some standard quadrature rules used in the `FEM` module, this module also
defines singular integration routines useful for (weakly) singular integrands.

Exports:
$(EXPORTS)
"""
module Integration

using DocStringExtensions
using StaticArrays
using QuadGK
using LinearAlgebra

import QuadGK: quadgk
import HCubature: hcubature

using WaveProp.Utils
using WaveProp.Geometry
using WaveProp.Interpolation

import WaveProp.Geometry:
    domain,
    range,
    ambient_dimension,
    geometric_dimension,
    jacobian

import WaveProp.Interpolation: barycentric_lagrange_weights

export
    # abstract types
    AbstractQuadratureRule,
    AbstractSingularityHandler,
    # types
    Gauss,
    Trapezoidal,
    TrapezoidalP,
    Fejer,
    GaussLegendre,
    TensorProductQuadrature,
    Kress,
    KressP,
    IMT,
    Duffy,
    SingularQuadratureRule,
    # functions
    integrate,
    qnodes,
    qweights,
    qnormals,
    singular_weights,
    singular_quadrature,
    TensorProductSingularityHandler

include("quadrule.jl")
include("quadgk.jl")
include("hcubature.jl")
include("singularityhandler.jl")
include("singularquadrule.jl")

end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Test
import ClimaCore.DataLayouts: VIFH, slab, column

@testset "VIFH" begin
    TestFloatTypes = (Float32, Float64)
    for FT in TestFloatTypes
        S = Tuple{Complex{FT}, FT}
        Nv = 10 # number of vertical levels
        Ni = 4  # number of nodal points
        Nh = 10 # number of elements

        # construct a data object with 10 cells in vertical and
        # 10 elements in horizontal with 4 nodal points per element in horizontal
        array = rand(FT, Nv, Ni, 3, Nh)

        data = VIFH{S, Ni}(array)
        sum(x -> x[2], data)

        @test getfield(data.:1, :array) == @view(array[:, :, 1:2, :])
        @test getfield(data.:2, :array) == @view(array[:, :, 3:3, :])

        @test size(data) == (Ni, 1, 1, Nv, Nh)

        # test tuple assignment on columns
        val = (Complex{FT}(-1.0, -2.0), FT(-3.0))
        column(data, 1, 1)[1] = val
        @test array[1, 1, 1, 1] == -1.0
        @test array[1, 1, 2, 1] == -2.0
        @test array[1, 1, 3, 1] == -3.0

        # test value of assing tuple on slab
        sdata = slab(data, 1, 1)
        @test sdata[1] == val

        # sum of all the first field elements
        @test sum(data.:1) ≈
              Complex{FT}(sum(array[:, :, 1, :]), sum(array[:, :, 2, :]))
        @test sum(x -> x[2], data) ≈ sum(array[:, :, 3, :])
    end
end

@testset "VIFH type safety" begin
    Nv = 1 # number of vertical levels
    Ni = 1  # number of nodal points per element
    Nh = 1 # number of elements

    # check that types of the same bitstype throw a conversion error
    SA = (a = 1.0, b = 2.0)
    SB = (c = 1.0, d = 2.0)

    array = zeros(Float64, Nv, Ni, 2, Nh)
    data = VIFH{typeof(SA), Ni}(array)

    cdata = column(data, 1, 1)
    cdata[1] = SA
    @test cdata[1] isa typeof(SA)
    @test_throws MethodError cdata[1] = SB

    sdata = slab(data, 1, 1)
    @test sdata[1] isa typeof(SA)
    @test_throws MethodError sdata[1] = SB
end

@testset "broadcasting between VIFH data object + scalars" begin
    FT = Float64
    data1 = ones(FT, 2, 2, 2, 2)
    S = Complex{Float64}
    data1 = VIFH{S, 2}(data1)
    res = data1 .+ 1
    @test res isa VIFH{S}
    @test parent(res) ==
          FT[f == 1 ? 2 : 1 for i in 1:2, j in 1:2, f in 1:2, h in 1:2]
    @test sum(res) == Complex(16.0, 8.0)
    @test sum(Base.Broadcast.broadcasted(+, data1, 1)) == Complex(16.0, 8.0)
end

@testset "broadcasting between VF + IFH data object => VIFH" begin
    FT = Float64
    S = Complex{FT}
    data_vf = VF{S}(ones(FT, 3, 2))
    data_ifh = IFH{FT, 2}(ones(FT, 2, 1, 2))
    data_vifh = data_vf .+ data_ifh
    @test data_vifh isa VIFH{S}
    @test size(data_vifh) == (2, 1, 1, 3, 2)
    @test parent(data_vifh) ==
          FT[f == 1 ? 2 : 1 for v in 1:3, i in 1:2, f in 1:2, h in 1:2]

    @test parent(data_vifh .+ data_vf) ==
          FT[f == 1 ? 3 : 2 for v in 1:3, i in 1:2, f in 1:2, h in 1:2]
    @test parent(data_vifh .+ data_ifh) ==
          FT[f == 1 ? 3 : 1 for v in 1:3, i in 1:2, f in 1:2, h in 1:2]

end


text:
immutable UTF16String <: String
    data::Array{Uint16,1} # includes 16-bit NULL termination after string chars
    function UTF16String(data::Vector{Uint16})
        if length(data) < 1 || data[end] != 0
            throw(ArgumentError("UTF16String data must be NULL-terminated"))
        end
        new(data)
    end
end

utf16_is_lead(c::Uint16) = (c & 0xfc00) == 0xd800
utf16_is_trail(c::Uint16) = (c & 0xfc00) == 0xdc00
utf16_is_surrogate(c::Uint16) = (c & 0xf800) == 0xd800
utf16_get_supplementary(lead::Uint16, trail::Uint16) = char((lead-0xd7f7)<<10 + trail)

function endof(s::UTF16String)
    d = s.data
    i = length(d) - 1
    i == 0 && return i
    utf16_is_surrogate(d[i]) ? i-1 : i
end
function next(s::UTF16String, i::Int)
    if !utf16_is_surrogate(s.data[i])
        return char(s.data[i]), i+1
    elseif length(s.data)-1 > i && utf16_is_lead(s.data[i]) && utf16_is_trail(s.data[i+1])
        return utf16_get_supplementary(s.data[i], s.data[i+1]), i+2
    end
    error("invalid UTF-16 character index")
end

# TODO: optmize this
function encode16(s::String)
    buf = Uint16[]
    for c in s
        if c < 0x10000
            push!(buf, uint16(c))
        else
            push!(buf, uint16(0xd7c0 + (c>>10) & 0x3ff))
            push!(buf, uint16(0xdc00 + c & 0x3ff))
        end
    end
    push!(buf, 0) # NULL termination
    UTF16String(buf)
end

utf16(x) = convert(UTF16String, x)
convert(::Type{UTF16String}, s::UTF16String) = s
convert(::Type{UTF16String}, s::String) = encode16(s)
convert(::Type{Array{Uint16,1}}, s::UTF16String) = s.data
convert(::Type{Array{Uint16}}, s::UTF16String) = s.data

# TODO: optimize this
convert(::Type{UTF8String}, s::UTF16String) =
    sprint(length(s.data)-1, io->for c in s; write(io,c::Char); end)

sizeof(s::UTF16String) = sizeof(s.data) - sizeof(Uint16)
convert{T<:Union(Int16,Uint16)}(::Type{Ptr{T}}, s::UTF16String) =
    convert(Ptr{T}, pointer(s))

function is_valid_utf16(data::AbstractArray{Uint16})
    i = 1
    n = length(data) # this may include NULL termination; that's okay
    while i < n # check for unpaired surrogates
        if utf16_is_lead(data[i]) && utf16_is_trail(data[i+1])
            i += 2
        elseif utf16_is_surrogate(data[i])
            return false
        else
            i += 1
        end
    end
    return i > n || !utf16_is_surrogate(data[i])
end

is_valid_utf16(s::UTF16String) = is_valid_utf16(s.data)

function convert(::Type{UTF16String}, data::AbstractVector{Uint16})
    !is_valid_utf16(data) && throw(ArgumentError("invalid UTF16 data"))
    len = length(data)
    d = Array(Uint16, len + 1)
    d[end] = 0 # NULL terminate
    UTF16String(copy!(d,1, data,1, len))
end

convert(T::Type{UTF16String}, data::AbstractArray{Uint16}) =
    convert(T, reshape(data, length(data)))

convert(T::Type{UTF16String}, data::AbstractArray{Int16}) =
    convert(T, reinterpret(Uint16, data))

function convert(T::Type{UTF16String}, bytes::AbstractArray{Uint8})
    isempty(bytes) && return UTF16String(Uint16[0])
    isodd(length(bytes)) && throw(ArgumentError("odd number of bytes"))
    data = reinterpret(Uint16, bytes)
    # check for byte-order mark (BOM):
    if data[1] == 0xfeff        # native byte order
        d = Array(Uint16, length(data))
        copy!(d,1, data,2, length(data)-1)
    elseif data[1] == 0xfffe    # byte-swapped
        d = Array(Uint16, length(data))
        for i = 2:length(data)
            d[i-1] = bswap(data[i])
        end
    else
        d = Array(Uint16, length(data) + 1)
        copy!(d,1, data,1, length(data)) # assume native byte order
    end
    d[end] = 0 # NULL terminate
    !is_valid_utf16(d) && throw(ArgumentError("invalid UTF16 data"))
    UTF16String(d)
end

utf16(p::Ptr{Uint16}, len::Integer) = utf16(pointer_to_array(p, len))
utf16(p::Ptr{Int16}, len::Integer) = utf16(convert(Ptr{Uint16}, p), len)
function utf16(p::Union(Ptr{Uint16}, Ptr{Int16}))
    len = 0
    while unsafe_load(p, len+1) != 0; len += 1; end
    utf16(p, len)
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#function merit_function_anticipated_reduction(iter, dir)
#    return norm(comp(iter),Inf)^3 / (iter.point.mu)^2 + norm(dir.x,2)^2 * get_delta(iter)
#end

mutable struct Class_filter_ls <: abstract_ls_info
    step_size_P::Float64
    step_size_D::Float64
    num_steps::Int64
    frac_progress::Float64
    predict_red::Float64
    actual_red::Float64
    cur_merit::Float64
    do_ls::Bool

    function Class_filter_ls(iter::Class_iterate, dir::Class_point, pars::Class_parameters)
        this = new()
        this.step_size_D = 0.0
        this.step_size_P = 0.0
        this.num_steps = 0

        g = eval_grad_phi(iter, get_mu(iter))
        this.predict_red = -comp_merit(iter) + 0.5 * (dot(g,dir.x) - get_delta(iter) * LinearAlgebra.norm(dir.x,2)^2)
        #merit_function_predicted_reduction(iter, dir, 1.0);
        this.cur_merit = eval_merit_function(iter, pars)

        sufficient_descent = dot(g,dir.x) < 0.0 #-0.5 * LinearAlgebra.norm(g,2)^2 / LinearAlgebra.norm(get_lag_hess(iter),2)^2
        merit_reduce = merit_function_predicted_reduction(iter, dir, 1.0) < this.predict_red / 2.0
        #comp_sufficient = -get_mu(iter) < minimum(comp_predicted(iter,dir,1.0)) && maximum(comp_predicted(iter,dir,1.0)) < 100.0 * get_mu(iter)
        #comp_sufficient = -get_mu(iter) * pars.ls.comp_feas_agg < minimum(comp_predicted(iter,dir,1.0)) && maximum(comp_predicted(iter,dir,1.0)) < get_mu(iter) * (1.0 / pars.ls.comp_feas_agg - 1.0)
        this.do_ls = sufficient_descent #&& merit_reduce
        #&& comp_sufficient
        #-0.5 * merit_function_anticipated_reduction(iter, dir)

        return this
    end

    function Class_filter_ls(step_size_P::Float64,step_size_D::Float64,num_steps::Int64,frac_progress::Float64, predict_red::Float64)
        this = new(step_size_P, step_size_D, num_steps, frac_progress, predict_red)
        return this
    end
end

mutable struct Class_filter
    fval::Float64
    scaled_kkt_err::Float64
    mu::Float64

    function Class_filter(iter::Class_iterate, pars::Class_parameters)
      this = new()
      this.fval = eval_merit_function(iter, pars)

      kkt_err = LinearAlgebra.norm(eval_grad_lag(iter, get_mu(iter)),Inf)
      if pars.ls.kkt_include_comp
         kkt_err += LinearAlgebra.norm(comp(iter),Inf)
      end
      this.scaled_kkt_err = kkt_err * dual_scale(iter, pars) #eval_kkt_err(iter, pars)
      this.mu = iter.point.primal_scale

      return this
    end
end

function add!(ar::Array{Class_filter,1}, iter::Class_iterate, pars::Class_parameters )
    push!( ar, Class_filter(iter, pars) )
end

function satisfies_filter!(ar::Array{Class_filter,1}, can::Class_iterate, step_size::Float64, pars::Class_parameters, timer::class_advanced_timer)
    p = Class_filter(can, pars)

    failure = false

    for i = 1:length(ar)
        #tol = 1e-3 * max(1.0, p.scaled_kkt_err)
        expected_reduction = pars.ls.kkt_reduction_factor
        #
        kkt_reduction = (p.scaled_kkt_err / ar[i].scaled_kkt_err < (1.0 - pars.ls.kkt_reduction_factor * step_size))
        fval_reduction = p.fval < ar[i].fval - (p.scaled_kkt_err)^2
        fval_no_increase = p.fval < ar[i].fval + sqrt(p.scaled_kkt_err) #+ p.scaled_kkt_err^2
        net_reduction = p.scaled_kkt_err + p.fval < ar[i].fval + ar[i].scaled_kkt_err - (p.scaled_kkt_err)^2

        if pars.ls.filter_type == :default
            failure = !(p.mu < ar[i].mu || kkt_reduction)
        elseif pars.ls.filter_type == :test1
            failure = !(p.mu < ar[i].mu || kkt_reduction || fval_reduction)
        elseif pars.ls.filter_type == :test2
            failure = !(p.mu < ar[i].mu || (kkt_reduction && fval_no_increase)) # || fval_reduction)
        elseif pars.ls.filter_type == :test3
            failure = !(p.mu < ar[i].mu || net_reduction)
        else
            error("unknown filter type!!!")
        end

        if failure #&& p.fval < ar[i].fval + (p.scaled_kkt_err) ))
          # && p.fval < ar[i].fval + (p.scaled_kkt_err)))
            #println(p.scaled_kkt_err / ar[i].scaled_kkt_err)
            # p.fval < ar[i].fval - (p.scaled_kkt_err)^2 ||
            return false
        end
    end

    return true
end

function accept_func!(accept::Class_filter_ls, iter::Class_iterate, candidate::Class_iterate, dir::Class_point, step_size::Float64, filter::Array{Class_filter,1}, pars::Class_parameters, timer::class_advanced_timer)
    status, suggested_step_size_P = accept_func_stable!(accept, iter, candidate, dir, step_size, pars, timer)

    if status == :success || status == :infeasible
      return status, suggested_step_size_P
    else
      return accept_func_kkt!(accept, iter, candidate, dir, step_size, filter, pars, timer)
    end
end