text:
# TODO: Add a more thorough explanation of the dependency confusion
# vulnerability and how this guideline mitigates it.

const guideline_dependency_confusion =
    Guideline("No UUID conflict with other registries.",
              data -> has_no_dependency_confusion(data.pkg,
                                                  data.registry_head,
                                                  data.public_registries))

# TODO: Needs a strategy to handle connection failures for the public
# registries. Preferably they should also be cloned only once and then
# just updated to mitigate the effect of them being temporarily
# offline. This could be implemented with the help of the Scratch
# package, but requires Julia >= 1.5.
function has_no_dependency_confusion(pkg, registry_head, public_registries)
    # We know the name of this package but not its uuid. Look it up in
    # the registry that includes the current PR.
    packages = TOML.parsefile(joinpath(registry_head, "Registry.toml"))["packages"]
    filter!(packages) do (key, value)
        value["name"] == pkg
    end
    # For Julia >= 1.4 this can be simplified with the `only` function.
    always_assert(length(packages) == 1)
    uuid = first(keys(packages))
    # Also need to find out the package repository.
    package_repo = TOML.parsefile(joinpath(registry_head, packages[uuid]["path"], "Package.toml"))["repo"]
    for repo in public_registries
        try
            registry = clone_repo(repo)
            registry_toml = TOML.parsefile(joinpath(registry, "Registry.toml"))
            packages = registry_toml["packages"]
            if haskey(packages, uuid)
                message = string("UUID $uuid conflicts with the package ",
                                 packages[uuid]["name"], " in registry ",
                                 registry_toml["name"], " at $repo. ",
                                 "This could be a dependency confusion attack.")
                # Conflict detected. This is benign if the package name
                # *and* the package URL matches.
                if packages[uuid]["name"] != pkg
                    return false, message
                end
                package_path = packages[uuid]["path"]
                other_package_repo = TOML.parsefile(joinpath(registry, package_path, "Package.toml"))["repo"]
                if package_repo != other_package_repo
                    return false, message
                end
            end
        catch
            message = string("Failed to clone public registry $(repo) for a check against dependency confusion.\n",
                             "This is an internal issue with the AutoMerge process and has nothing to do with ".
                             "the package being registered but requires manual intervention before AutoMerge ",
                             "can be resumed.")
            return false, message
        end
    end

    return true, ""
end


text:
using Test, JuLIP, JuLIPMaterials
using JuLIPMaterials: ForceConstantMatrix1

datapath = joinpath(dirname(pathof(JuLIP)), "..", "data")
eam_Fe = EAM(datapath * "/pfe.plt",
             datapath * "/ffe.plt",
             datapath * "/F_fe.plt")

# equilibrate a unit cell
fe1 = bulk(:Fe)
set_calculator!(fe1, eam_Fe)
variablecell!(fe1)
minimise!(fe1)

# get the force constants
fcm = ForceConstantMatrix1(eam_Fe, fe1, h = 1e-5)

println("compute exact hessian on a large cell")
at = fe1 * 12
fixedcell!(at)
set_calculator!(at, eam_Fe)
H = JuLIP.hessian_pos(eam_Fe, at)
Hvec = JuLIP.hessian(at)

println("random virtual displacement => compare H * u and fcm * u")
U = rand(JVecF, length(at))
Uvec = mat(U)[:]
V = vecs(Hvec * Uvec)
# sanity check
@test V[1] ≈ sum( H[1,n] * U[n] for n = 1:length(U) )

# now multiply FCM * U
Vfcm = fcm * (at, U)
err_HxU = maximum(norm.(Vfcm - V)) / maximum(norm.(V))
@show err_HxU
println(@test err_HxU < 1e-7)

println("Check that individual blocks match")
X = positions(at)
x̄ = sum(X) / length(X)
r = [norm(x-x̄) for x in X]
n0 = findall( r .== minimum(r) )[1]
all_found = true
all_passed = true
for n = 1:length(at)
   n == n0 && continue
   norm(H[n0,n]) < 1e-12 && continue
   R = X[n] - X[n0]
   found = false
   for m = 1:length(fcm.R)
      if norm(R - fcm.R[m]) < 1e-7
         found = true
         if norm(H[n0, n] - fcm.H[m]) < 1e-5
            print("+")
         else
            global all_passed = false
            print("-")
         end
         break
      end
   end
   if !found
      global all_found = false
      print("x")
   end
end
println(@test all_passed)
println(@test all_found )

ext: jl


text:
runClustering(dataArray,dataMask,denois;minSize=1,alpha=0.05,doChildrenFirst=false)=runClustering(CPU1(),dataArray,dataMask,denois,minSize=minSize,alpha=alpha,doChildrenFirst=doChildrenFirst)

function runClustering(resource::CPU1,dataArray,dataMask,denois;minSize=1,alpha=0.05,doChildrenFirst=false)

    dim = size(dataArray)
    ndim = ndims(dataArray)
    ntime = dim[ndim]
    coord = dim[1:(ndim-1)]
    fullLength = prod(coord)
    dataMaskInd = (LinearIndices(dataMask))[findall(dataMask)]
    nvox = length(dataMaskInd)
    dataArray = copy(transpose((reshape(dataArray,(fullLength,ntime))[dataMaskInd,:])))
    iter = floor(Integer,log2(ntime))-1
    from = round(Int,exp2(iter+1))
    realIter = floor(Integer,log2(from-1))-1
    quant = alpha/(iter+1)
    thrs = [quantile(Chisq(i),1-quant) for i=exp2.(0:(realIter-1))]

    infoDen = denois["infoDen"]
    dataProj = denois["dataProj"]
    dataVar = denois["var"]

    toCluster = ones(Int64,length(infoDen))

    if doChildrenFirst
         ## Build the children list of all denoised voxels
        childrenDict = pmap(1:length(infoDen)) do idx
            getChildren(idx,infoDen,dataProj,dataVar,thrs,toCluster)
        end
    end

    badVoxels = Int64[]
    clusterDict = Dict[]
    actualMinSize = Inf

    while any(toCluster.==1)
        idToCluster = (LinearIndices(toCluster))[findall(toCluster)]#find(toCluster)
        println(length(idToCluster), "---", actualMinSize)
        ## Find in the remaining voxels to clusterize
        ## that with the largest list of children (if doChildrenFirst=true) or neighbors (otherwise)
        ## remaining to clusterize
        if doChildrenFirst
            sizeto = map(childrenDict[idToCluster]) do x
                sum(toCluster[x])
            end
        else
            sizeto = map(infoDen[idToCluster]) do x
                sum(toCluster[x["Vx"]])
            end
        end

        idMax = argmax(sizeto)
        actualMinSize = sizeto[idMax]
        if actualMinSize<minSize
            break
        end
        idCl = idToCluster[idMax]

        ## build cluster from voxel defined by id.tocluster
        ## find the active children (still in the list of voxels to clusterize) of id.tocluster
        if doChildrenFirst
            childrenOk = childrenDict[idCl]
        else
            childrenOk = getChildren(idCl,infoDen,dataProj,dataVar,thrs,toCluster)
        end
        ## children are active children if toCluster[children] is true
        if length(childrenOk)==0
            toCluster[idCl] = 0
            println("nothing to do")
            continue
        end

        inChildrenLength = length(childrenOk)
        if inChildrenLength <= minSize
            ## not enough children, we just add voxel to closest cluster
            #addToClosestCluster(idCl,toCluster,clusterDict,badVoxels,dataProj,dataVar,quant,iter)
            toCluster[idCl] = 0
            ## add voxel to the closest cluster i.e. with largest p-value
            ind = checkClusterNew(dataProj[:,idCl],idCl,clusterDict,dataVar,quant,realIter,pvalueMax=true)

            if !isa(ind,Bool)
                clusterDict[ind]["cluster"] = [clusterDict[ind]["cluster"];idCl]
                clusterDict[ind]["center"] = mean(dataProj[:,clusterDict[ind]["cluster"]],dims=2)
            end
            badVoxels = [badVoxels;idCl]
            continue
        end

        check = 1
        connectedClusters = Dict[]
        newCluster = Dict[]

        while length(check)>0
           newCluster = robustMean(childrenOk,dataProj,dataVar)
            #println(length(newCluster["cluster"]))
           check = checkClusterNew(newCluster["center"],newCluster["cluster"],clusterDict,dataVar,quant,realIter,pvalueMax=false,useFdr=true)
           # println(check)
           if length(check)>0
           ## at least, one existing cluster is coherent with the new one
           ## get the connected clusters
               connectedClusters = [connectedClusters;clusterDict[check]]
           ## get the connected voxels id
               #println(length(clusterDict[check][1]["cluster"]))
               connectedIdx = mapreduce((x) ->  x["cluster"],vcat,clusterDict[check],init=[])
               #println(length(connectedIdx))
           ## remove the connected clusters from the cluster list
               deleteat!(clusterDict,sort(check))
           ## add the connected in the neighbors of idx.tocluster
               childrenOk = unique([childrenOk;connectedIdx])
           end
        end

        if length(newCluster["outliers"])<inChildrenLength
            ## we found a new cluster !!!
            ## remove the new cluster from pixel to clusterize
            toCluster[newCluster["cluster"]] = 0
            toCluster[newCluster["outliers"]] = 1
            ## update cluster list
            clusterDict = [clusterDict;Dict("cluster"=>newCluster["cluster"],"center"=>newCluster["center"])]
            #println(length(clusterDict))
        else
            ## TOO MANY OUTLIERS: we do not build new cluster
            ## Return cluster list to its original state
            ## Add id.cluster to the closest cluster
            ## return cluster list in its original state
            clusterDict = [clusterDict;connectedClusters]
            ## add voxel to closest cluster
            ## mark voxel as already clusterize
            toCluster[idCl] = 0

            ## add voxel to the closest cluster i.e. with largest p-value
            ind = checkClusterNew(dataProj[:,idCl],idCl,clusterDict,dataVar,quant,realIter,pvalueMax=true)
            if !isa(ind,Bool)
                clusterDict[ind]["cluster"] = [clusterDict[ind]["cluster"];idCl]
                clusterDict[ind]["center"] = mean(dataProj[:,clusterDict[ind]["cluster"]],dims=2)
            end
            badVoxels = [badVoxels;idCl]
        end

    end
    ## sort the cluster list in decreasing order
    sort!(clusterDict,by = ((x) -> length(x["cluster"])),rev=true)
    #println(length(clusterDict))
    ## compute the cluster centers
    for i in 1:length(clusterDict)
        clusterDict[i]["center"]=mean(dataArray[:,clusterDict[i]["cluster"]],dims=2)
    end

    ## extract centers and clusters
    centers = map((x) -> x["center"],clusterDict)
    clusters = map((x) ->  x["cluster"],clusterDict)

    return(Dict("clusters" => clusters, "centers" => centers, "badVoxels" => badVoxels))
end

function getClusteringResults(dataArray,denois,clust)
    clustMap = 0*Array{Int64}(undef,size(dataArray)[1:(end-1)])
    #outArray = dataArray
    for i in 1:length(clust["clusters"])
        for j in clust["clusters"][i]
            tmp = denois["infoDen"][j]
            if tmp == nothing
                continue
            end
            mid = tmp["Cx"]
     #       outArray[mid,:] = clust["centers"][i]
            clustMap[mid] = i
        end
    end
    # Dict("clustArray" => outArray,"clustMap"=>clustMap,"clustCenter"=>clust["centers"])
    Dict("clustMap"=>clustMap,"clustCenter"=>clust["centers"])
end


text:
@testset "produce" begin
    @testset "basic" begin
        @system SProduce begin
            a => produce(SProduce) ~ produce
        end
        @system SProduceController(Controller) begin
            s(context) ~ ::SProduce
        end
        sc = instance(SProduceController)
        s = sc.s
        @test length(s.a) == 0
        @test collect(s.a) == []
        update!(sc)
        @test length(s.a) == 1
        @test collect(s.a) == [s.a[1]]
        @test length(s.a[1].a) == 0
        update!(sc)
        @test length(s.a) == 2
        @test collect(s.a) == [s.a[1], s.a[2]]
        @test length(s.a[1].a) == 1
        @test length(s.a[2].a) == 0
    end

    @testset "single" begin
        @system SProduceSingle begin
            a => produce(SProduceSingle) ~ produce::SProduceSingle
        end
        @system SProduceSingleController(Controller) begin
            s(context) ~ ::SProduceSingle
        end
        sc = instance(SProduceSingleController)
        s = sc.s
        @test length(s.a) == 0
        @test collect(s.a) == []
        update!(sc)
        @test length(s.a) == 1
        a = s.a[1]
        @test collect(s.a) == [a]
        @test length(s.a[1].a) == 0
        update!(sc)
        @test length(s.a) == 1
        @test collect(s.a) == [a]
        @test a === s.a[1]
        @test length(s.a[1].a) == 1
        @test_throws BoundsError s.a[2]
    end

    @testset "kwargs" begin
        @system SProduceKwargs begin
            a => produce(SProduceKwargs) ~ produce
            i(t=context.clock.time) => t ~ preserve(u"hr")
        end
        @system SProduceKwargsController(Controller) begin
            s(context) ~ ::SProduceKwargs
        end
        sc = instance(SProduceKwargsController)
        s = sc.s
        @test length(s.a) == 0 && s.i' == 0u"hr"
        update!(sc)
        @test length(s.a) == 1 && s.i' == 0u"hr"
        @test length(s.a[1].a) == 0 && s.a[1].i' == 0u"hr"
        update!(sc)
        @test length(s.a) == 2 && s.i' == 0u"hr"
        @test length(s.a[1].a) == 1 && s.a[1].i' == 0u"hr"
        @test length(s.a[2].a) == 0 && s.a[2].i' == 1u"hr"
        @test length(s.a[1].a[1].a) == 0 && s.a[1].a[1].i' == 1u"hr"
    end

    @testset "nothing" begin
        @system SProduceNothing begin
            a => nothing ~ produce
        end
        @system SProduceNothingController(Controller) begin
            s(context) ~ ::SProduceNothing
        end
        sc = instance(SProduceNothingController)
        s = sc.s
        @test length(s.a) == 0
        update!(sc)
        @test length(s.a) == 0
    end

    @testset "query index" begin
        @system SProduceQueryIndex begin
            p => produce(SProduceQueryIndex) ~ produce::SProduceQueryIndex[]
            i(context.clock.tick) ~ preserve::int
            a(x=p["*"].i) => sum(x) ~ track
            b(x=p["**"].i) => sum(x) ~ track
        end
        @system SProduceQueryIndexController(Controller) begin
            s(context) ~ ::SProduceQueryIndex
        end
        sc = instance(SProduceQueryIndexController)
        s = sc.s
        @test length(s.p) == 0
        update!(sc)
        @test length(s.p) == 1
        @test s.a' == 0 # (0)
        @test s.b' == 0 # (0)
        update!(sc)
        @test length(s.p) == 2
        @test s.a' == 1 # (0 + 1)
        @test s.b' == 2 # ((0 ~ 1) + 1)
        update!(sc)
        @test length(s.p) == 3
        @test s.a' == 3 # (0 + 1 + 2)
        @test s.b' == 10 # ((0 ~ ((1 ~ 2) + 2) + (1 ~ 2) + 2)
    end

    @testset "query condition with flag" begin
        @system SProduceQueryConditionTrackBool begin
            p => produce(SProduceQueryConditionTrackBool) ~ produce::SProduceQueryConditionTrackBool[]
            i(context.clock.tick) ~ preserve::int
            f(i) => isodd(i) ~ flag
            a(x=p["*/f"].i) => sum(x) ~ track
            b(x=p["**/f"].i) => sum(x) ~ track
        end
        @system SProduceQueryConditionTrackBoolController(Controller) begin
            s(context) ~ ::SProduceQueryConditionTrackBool
        end
        sc = instance(SProduceQueryConditionTrackBoolController)
        s = sc.s
        @test length(s.p) == 0
        update!(sc)
        @test length(s.p) == 1
        @test s.a' == 0 # (#0)
        @test s.b' == 0 # (#0)
        update!(sc)
        @test length(s.p) == 2
        @test s.a' == 1 # (0 + #1)
        @test s.b' == 2 # (0 ~ #1) + #1)
        update!(sc)
        @test length(s.p) == 3
        @test s.a' == 1 # (0 + #1 + 2)
        @test s.b' == 2 # (0 ~ ((#1 ~ 2) + 2) + (#1 ~ 2) + 2)
        update!(sc)
        @test length(s.p) == 4
        @test s.a' == 4 # (0 + #1 + 2 + #3)
        @test s.b' == 26 # (0 ~ ((#1 ~ ((2 ~ #3) + #3)) + (2 ~ #3) + #3) + (#1 ~ ((2 ~ #3) + #3)) + (2 ~ #3) + #3)
    end

    @testset "adjoint" begin
        @system SProduceAdjoint begin
            p => produce(SProduceAdjoint) ~ produce::SProduceAdjoint[]
            i(context.clock.tick) ~ preserve::int
        end
        @system SProduceAdjointController(Controller) begin
            s(context) ~ ::SProduceAdjoint
        end
        sc = instance(SProduceAdjointController)
        s = sc.s
        update!(sc)
        @test length(s.p["*"]') == 1
        @test length(s.p["**"]') == 1
        @test s.p["*"].i' == [0]
        @test s.p["**"].i' == [0]
        update!(sc)
        @test length(s.p["*"]') == 2
        @test length(s.p["**"]') == 3
        @test s.p["*"].i' == [0, 1]
        @test s.p["**"].i' == [0, 1, 1]
    end

    @testset "when" begin
        @system SProduceWhen begin
            t(context.clock.tick) ~ track::int
            w(t) => isodd(t) ~ flag
            a => produce(SProduceWhen) ~ produce(when=w)
        end
        @system SProduceWhenController(Controller) begin
            s(context) ~ ::SProduceWhen
        end
        sc = instance(SProduceWhenController)
        s = sc.s
        @test length(s.a) == 0
        update!(sc)
        @test length(s.a) == 0
        update!(sc)
        @test length(s.a) == 1
        @test length(s.a[1].a) == 0
        update!(sc)
        @test length(s.a) == 1
        @test length(s.a[1].a) == 0
        update!(sc)
        @test length(s.a) == 2
        @test length(s.a[1].a) == 1
        @test length(s.a[2].a) == 0
        update!(sc)
        @test length(s.a) == 2
        @test length(s.a[1].a) == 1
        @test length(s.a[2].a) == 0
    end
end


text:
module FastUnitfulOperations

using LinearAlgebra
using Unitful

export *ᵘ, uproduct, uproduct!, usort, usort!, usortperm, usortperm!

"""
    uproduct(a::AbstractArray{T1}, b::AbstractArray{T2}) where {T1 and/or T2 <:Quantity}
    *ᵘ(a::AbstractArray{T1}, b::AbstractArray{T2}) where {T1 and/or T2 <:Quantity}

Fast matrix multiplication for Unitful arrays.

This method reinterprets the unitful arrays as arrays with the number's type to make use of
optimized matrix multiplication for integers or floats.
The method will throw an error if the array contains Unitful elements that do not have a fixed
value or unit type.

Type `*\\^u` followed by the tab key to input the operator `*ᵘ`.
"""
function uproduct(a::AbstractArray{T1}, b::AbstractArray{T2}) where {T1<:Quantity,T2<:Quantity}
    c = ustrip(a) * ustrip(b)
    u = unit(T1) * unit(T2)
    return reinterpret(Quantity{eltype(c),dimension(u),typeof(u)}, c)
end

function uproduct(a::AbstractArray{T1}, b::AbstractArray) where {T1<:Quantity}
    c = ustrip(a) * b
    return reinterpret(T1, c)
end

function uproduct(a::AbstractArray, b::AbstractArray{T2}) where {T2<:Quantity}
    c = a * ustrip(b)
    return reinterpret(T2, c)
end

"""
    uproduct!(c::AbstractArray, a::AbstractArray{T1}, b::AbstractArray{T2}) where {T1 and/or T2 <:Quantity}

Fast in-place matrix multiplication for Unitful arrays.

This method reinterprets the unitful arrays as arrays with the number's type to make use of
optimized matrix multiplication for integers or floats.
The method will throw an error if the array contains Unitful elements that do not have a fixed
value or unit type.
"""
function uproduct!(c::AbstractArray{T1}, a::AbstractArray{T2}, b::AbstractArray{T3}) where {T1<:Quantity,T2<:Quantity,T3<:Quantity}
    mul!(ustrip(c), ustrip(a), ustrip(b))
    return c
end

function uproduct!(c::AbstractArray, a::AbstractArray{T2}, b::AbstractArray{T3}) where {T2<:Quantity,T3<:Quantity}
    mul!(c, ustrip(a), ustrip(b))
    return c
end

function uproduct!(c::AbstractArray{T1}, a::AbstractArray, b::AbstractArray{T3}) where {T1<:Quantity,T3<:Quantity}
    mul!(ustrip(c), a, ustrip(b))
    return c
end

function uproduct!(c::AbstractArray{T1}, a::AbstractArray{T2}, b::AbstractArray) where {T1<:Quantity,T2<:Quantity}
    mul!(ustrip(c), ustrip(a), b)
    return c
end

"""
    uproduct(a, b)
    *ᵘ(a, b)

Fallback to regular multiplication: `uproduct(a, b) = a * b`
"""
uproduct(a, b) = a * b

"""
    uproduct!(c, a, b)

Fallback to regular multiplication: `uproduct!(c, a, b) = mul!(c, a, b)`
"""
uproduct!(c, a, b) = mul!(c, a, b)

*ᵘ = uproduct

"""
    usort(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    usort(a; kwargs...)

Sort a unitful array as a reinterpreted array according to the number's type.

Calls [`Base.sort`](https://docs.julialang.org/en/v1/base/sort/#Base.sort) on the reinterpreted
array. This enables the use of sorting algorithms that only work for integers or floats
(e.g. RadixSort from
[SortingAlgorithms](https://github.com/JuliaCollections/SortingAlgorithms.jl).

The method falls back `Base.sort` for other array types.
"""
function usort(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    a_sorted = sort(ustrip(a); kwargs...)
    return reinterpret(T, a_sorted)
end

"""
    usort!(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    usort!(a; kwargs...)

Sort a unitful array in-place as a reinterpreted array according to the number's type.

Calls [`Base.sort!`](https://docs.julialang.org/en/v1/base/sort/#Base.sort!) on the reinterpreted
array. This enables the use of sorting algorithms that only work for integers or floats
(e.g. RadixSort from
[SortingAlgorithms](https://github.com/JuliaCollections/SortingAlgorithms.jl).

The method falls back `Base.sort!` for other array types.
"""
function usort!(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    sort!(ustrip(a); kwargs...)
    return a
end

"""
    usortperm(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    usortperm(a; kwargs...)

Return the permutation vector to sort a unitful array as a reinterpreted array according to the
number's type.

Calls [`Base.sortperm`](https://docs.julialang.org/en/v1/base/sort/#Base.sortperm) on the reinterpreted
array. This enables the use of sorting algorithms that only work for integers or floats
(e.g. RadixSort from
[SortingAlgorithms](https://github.com/JuliaCollections/SortingAlgorithms.jl).

The method falls back `Base.sortperm` for other array types.
"""
function usortperm(a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    return sortperm(ustrip(a); kwargs...)
end

"""
    usortperm!(ix, a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    usortperm!(ix, a; kwargs...)

Return the permutation vector in the preallocated vector `ix` to sort a unitful array as a
reinterpreted array according to the number's type.

Calls [`Base.sortperm!`](https://docs.julialang.org/en/v1/base/sort/#Base.sortperm!) on the reinterpreted
array. This enables the use of sorting algorithms that only work for integers or floats
(e.g. RadixSort from
[SortingAlgorithms](https://github.com/JuliaCollections/SortingAlgorithms.jl).

The method falls back `Base.sortperm!` for other array types.
"""
function usortperm!(ix, a::AbstractArray{T}; kwargs...) where {T<:Quantity}
    return sortperm!(ix, ustrip(a); kwargs...)
end

usort(a; kwargs...) = sort(a; kwargs...)
usort!(a; kwargs...) = sort!(a; kwargs...)
usortperm(a; kwargs...) = sortperm(a; kwargs...)
usortperm!(ix, a; kwargs...) = sortperm!(ix, a; kwargs...)

end


text:
struct ThisPlaceholder end

const this = ThisPlaceholder()

Base.show(io::IO, ::ThisPlaceholder) = print(io, :this)
Base.show(io::IO, ::MIME"text/plain", ::ThisPlaceholder) = print(io, :this)

mutable struct FullFieldInfo
    name::Union{Symbol, Nothing}
    rawtype::Any
    line::Union{LineNumberNode, Missing}
    type::DataType
    construct::Construct
end

struct ParameterInfo
    name::Symbol
    rawtype::Any
    type::DataType
end

mutable struct ConstructContext
    structname::Any
    thisfields::Vector{FullFieldInfo}
    parameters::Vector{ParameterInfo}
end

const _lck = ReentrantLock()
const _context = ConstructContext(nothing, FullFieldInfo[], ParameterInfo[])

function withcontext(action::Function, structname)
    lock(_lck)
    _context.structname = structname
    _context.thisfields = FullFieldInfo[]
    _context.parameters = ParameterInfo[]
    try
        action(_context)
    finally
        _context.structname = nothing
        _context.thisfields = FullFieldInfo[]
        _context.parameters = ParameterInfo[]
        unlock(_lck)
    end
end
abstract type ContextualExpr{T} end

struct InstantValue{T} <: ContextualExpr{T}
    value::T
end

struct ThisProperty{T} <: ContextualExpr{T}
    prop::Symbol
end
ThisProperty(prop::Symbol) = let idx = findfirst(field => field.name == prop, _context.thisfields)
    if idx isa Int
        ThisProperty{_context.thisfields[idx].type}(prop)
    else
        error("type $(_context.structname) has no field $prop.")
    end
end

Base.getproperty(::ThisPlaceholder, name::Symbol) = ThisProperty(name)

struct Property{S, T} <: ContextualExpr{T}
    obj::ContextualExpr{S}
    prop::Symbol
end
Property(obj::ContextualExpr{S}, prop::Symbol) where {S} = Property{S, Base.fieldtype(S, prop)}(obj, prop)
Base.getproperty(obj::ContextualExpr, name::Symbol) = Property(obj, name)

struct FunctionCall{T} <: ContextualExpr{T}
    func::Function
    args::Vector{ContextualExpr}
end

macro construct(structdef::Expr)
    construct_impl(__module__, __source__, gensym("CustomConstruct"), structdef)
end

macro construct(constructname::Symbol, structdef::Expr)
    construct_impl(__module__, __source__, constructname, structdef)
end

struct FieldInfo
    name::Union{Symbol, Nothing}
    constype::Union{Construct, DataType}
    line::Union{LineNumberNode, Missing}
end

struct OtherFieldInfo
    expr::Any
    line::Union{LineNumberNode, Missing}
end

function construct_impl(mod::Module, source::LineNumberNode, constructname::Symbol, structdef::Expr)
    if structdef.head == :struct
        fields = dumpfields(mod, structdef)
        typedstructdef = replacestructdef(structdef, fields)
        constructdef = generateconstructdef(constructname, getdefname(structdef))
        Expr(:block,
            source, typedstructdef,
            source, constructdef)
    else
        error("invalid syntax: @construct must be used with a struct type definition.")
    end
end

function dumpfields(mod::Module, structdef::Expr)
    @assert structdef.head == :struct
    stfields = structdef.args[3].args
    fields = Vector{Union{FieldInfo, OtherFieldInfo}}()
    sizehint!(fields, length(stfields))
    for i in eachindex(stfields)
        field = stfields[i]
        if field isa LineNumberNode
            continue
        else
            line = stfields[i-1] isa LineNumberNode ? stfields[i-1] : missing
            if field isa Symbol
                push!(fields, FieldInfo(field, Any, line))
            elseif field isa Expr && field.head == :(::)
                constype = Core.eval(mod, field.args[end])
                if length(field.args) == 2 # x::Int
                    push!(fields, FieldInfo(field.args[1], constype, line))
                elseif length(field.args) == 1 # ::Padding(4)
                    push!(fields, FieldInfo(nothing, constype, line))
                else
                    push!(fields, OtherFieldInfo(field, line))
                end
            else
                push!(fields, OtherFieldInfo(field, line))
            end
        end
    end
    fields
end

function replacestructdef(structdef::Expr, fields::Vector{Union{FieldInfo, OtherFieldInfo}})
    @assert structdef.head == :struct
    stfields = Vector()
    sizehint!(stfields, length(structdef.args[3].args))
    for field in fields
        if field isa FieldInfo
            if isnothing(field.name) # omit field without name
                continue
            else
                if !ismissing(field.line)
                    push!(stfields, field.line)
                end
                push!(stfields, Expr(:(::), field.name, constructtype(field.constype)))
            end
        else
            if !ismissing(field.line)
                push!(stfields, field.line)
            end
            push!(stfields, field.expr)
        end
    end
    Expr(:struct, structdef.args[1], structdef.args[2], Expr(:block, stfields...))
end

function generateconstructdef(constructname::Symbol, structname)
    Expr(:block,
        Expr(:struct,
            false,
            Expr(:(<:), esc(constructname), Expr(:curly, GlobalRef(Constructs, :Construct), esc(structname))),
            Expr(:block),
        ),
        Expr(:function,
            Expr(:call,
                GlobalRef(Constructs, :Construct),
                Expr(:(::), Expr(:curly, GlobalRef(Core, :Type), esc(structname)))
            ),
            Expr(:block,
                Expr(:call, esc(constructname))
            ),
        )
    )
end

getdefname(name::Symbol) = name
getdefname(sym::GlobalRef) = sym.name
function getdefname(expr::Expr)
    if expr.head == :struct
        getdefname(expr.args[2])
    elseif expr.head in [:abstract, :(<:), :curly, :call, :function, :(=), :where]
        getdefname(expr.args[1])
    else
        error("syntax error: invalid definition $expr.")
    end
end

hasthis(sym::Symbol) = sym == :this
hasthis(sym::GlobalRef) = hasthis(sym.name)
hasthis(ex::Expr) = any(hasthis, ex.args)
hasthis(arr::AbstractArray) = any(hasthis, arr)
hasthis(tuple::Tuple) = any(hasthis, tuple)
hasthis(_) = false

iscontextfree(ex) = !hasthis(ex)


text:
#numshell.jl - testing if number shells working

type shell
  n::Int64
end

import Base.+, Base.-, Base.*, Base./
+(x::shell, y::shell) = shell(x.n + y.n)
-(x::shell, y::shell) = shell(x.n - y.n)
*(x::shell, y::shell) = shell(x.n * y.n)
/(x::shell, y::shell) = shell(x.n / y.n)

add!(x, y, d) = (d.n = x.n + y.n; nothing)
sub!(x, y, d) = (d.n = x.n - y.n; nothing)
mul!(x, y, d) = (d.n = x.n - y.n; nothing)
div!(x, y, d) = (d.n = x.n - y.n; nothing)

macro shellparse(expr)

end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using FFTW

struct Not{T}
  idx::T
end
import Base: to_indices, uncolon, tail, _maybetail
import LinearAlgebra: norm, dot, mul!

@inline to_indices(A, inds, I::Tuple{Not, Vararg{Any}}) =
   (setdiff(uncolon(inds, (:, tail(I)...)), I[1].idx), to_indices(A, _maybetail(inds), tail(I))...)

@testset "Grid Routines" begin

  # size
  nx = 12; ny = 12

  # sample point
  i = 5; j = 7
  cellzero = Nodes(Dual,(nx,ny))
  nodezero = Nodes(Primal,cellzero)
  facezero = Edges(Primal,cellzero)
  dualfacezero = Edges(Dual,cellzero)

  @test typeof(cellzero) <: ScalarGridData

  @test typeof(nodezero) <: ScalarGridData

  @test typeof(facezero) <: VectorGridData

  @test typeof(dualfacezero) <: VectorGridData

  @test typeof(facezero.u) <: XEdges{Primal}
  @test typeof(facezero.v) <: YEdges{Primal}

  @test typeof(dualfacezero.u) <: XEdges{Dual}
  @test typeof(dualfacezero.v) <: YEdges{Dual}

  @test celltype(cellzero) == Dual
  @test celltype(facezero) == Primal
  @test griddatatype(dualfacezero) == Edges

  @test griddatatype(NodePair(Dual,dualfacezero)) == NodePair
  @test griddatatype(EdgeGradient(Dual,cellzero)) == EdgeGradient

  cellzero2 = typeof(cellzero)()
  @test typeof(cellzero2) == typeof(cellzero)

  cellzero2 = zero(cellzero)
  @test typeof(cellzero2) == typeof(cellzero)

  cellunit = deepcopy(cellzero)
  cellunit[i,j] = 1.0

  nodeunit = deepcopy(nodezero)
  nodeunit[i,j] = 1.0

  facexunit = deepcopy(facezero)
  facexunit.u[i,j] = 1.0

  faceyunit = deepcopy(facezero)
  faceyunit.v[i,j] = 1.0

  dualfacexunit = deepcopy(dualfacezero)
  dualfacexunit.u[i,j] = 1.0

  dualfaceyunit = deepcopy(dualfacezero)
  dualfaceyunit.v[i,j] = 1.0

  @testset "Basic array operations" begin
    w = zero(cellunit)
    w .= cellunit
    @test w[i,j] == 1.0
    q = similar(facexunit)
    q .= facexunit
    @test q.u[i,j] == 1.0
    @test iszero(q.v)

    # providing a vector of data to GridData
    data = zeros(Float64,length(w))
    w2 = typeof(w)(data)
    w2[1,1] = 2.0
    @test vec(w2) == data

    data = zeros(Float64,length(q))
    q2 = typeof(q)(data)
    @test all(q2 .== 0.0)
    @test vec(q2) == data

    q2.u .= 2.0
    @test all(data[1:length(q2.u)] .== 2.0)

  end

  @testset "Arithmetic" begin

    w = 2*cellunit
    @test w[i,j] == 2
    w2 = w + w
    @test w2[i,j] == 4

    q = Edges(Dual,w)
    q.u .= 1.0
    q2 = 2*q
    @test all(q2.u .== 2.0)
    q2 .= q + q2
    @test all(q2.u .== 3.0)

    t = EdgeGradient(Dual,w)
    t .= 1.0
    t2 = 2*t
    @test all(t2 .== 2.0)
    t2 = -t
    @test all(t2 .== -1.0)

  end

  @testset "Inner products and norms" begin
    w = zero(cellunit)
    i0, j0 = rand(2:nx-1), rand(2:ny-1)
    w[i0,j0] = 1.0
    @test norm(w)*sqrt((nx-2)*(ny-2)) == 1.0
    w .= 1.0
    @test norm(w) == 1.0

    @test 2*w == w*2

    p = Nodes(Primal,w)
    p .= 1.0
    @test norm(p) == 1.0
    p2 = deepcopy(p)
    @test dot(p,p2) == 1.0
    @test norm(p-p2) == 0.0

    q = Edges(Dual,w)
    q.u .= 1.0
    q2 = deepcopy(q)
    @test dot(q,q2) == 1.0

    q = Edges(Primal,w)
    q.u .= 1.0
    q2 = deepcopy(q)
    @test dot(q,q2) == 1.0

    @test integrate(w) == 1.0

    @test integrate(p) == 1.0

    q .= 1
    @test norm(2*q) == sqrt(8)

    v = Edges(Primal,(100,200))
    v.u[3:end-3,3:end-3] .= rand(size(v.u[3:end-3,3:end-3])...)
    v.v[3:end-3,3:end-3] .= rand(size(v.v[3:end-3,3:end-3])...)

    dq = EdgeGradient(Primal,v)

    grid_interpolate!(dq,v)

    dq2 = similar(dq)
    dq2.dudx[3:end-3,3:end-3] .= rand(size(dq2.dudx[3:end-3,3:end-3])...)
    dq2.dudy[3:end-3,3:end-3] .= rand(size(dq2.dudy[3:end-3,3:end-3])...)
    dq2.dvdx[3:end-3,3:end-3] .= rand(size(dq2.dvdx[3:end-3,3:end-3])...)
    dq2.dvdy[3:end-3,3:end-3] .= rand(size(dq2.dvdy[3:end-3,3:end-3])...)

    v2 = similar(v)
    grid_interpolate!(v2,dq2)

    @test isapprox(dot(dq2,dq),dot(v2,v),atol=100*eps(1.0))

  end

  @testset "Dual cell center data Laplacian" begin
    lapcell = laplacian(cellunit)
    @test lapcell[i,j] == -4.0
    lapcell[i,j] = 0.0
    @test lapcell[i+1,j] == lapcell[i-1,j] == lapcell[i,j-1] == lapcell[i,j+1] == 1.0
    lapcell[i+1,j] = lapcell[i-1,j] = lapcell[i,j-1] = lapcell[i,j+1] = 0.0
    @test iszero(lapcell)
  end

  @testset "Dual cell center data curl" begin
    q = curl(cellunit)
    @test q.u[i,j-1] == 1.0 && q.u[i,j] == -1.0
    q.u[i,j-1] = q.u[i,j] = 0.0
    @test iszero(q.u)
    @test q.v[i-1,j] == -1.0 && q.v[i,j] == 1.0
    q.v[i-1,j] = q.v[i,j] = 0.0
    @test iszero(q.v)
  end

  @testset "Dual cell node gradient" begin
    q = grad(nodeunit)
    @test q.u[i,j] == 1.0 && q.u[i+1,j] == -1.0
    q.u[i,j] = q.u[i+1,j] = 0.0
    @test iszero(q.u)
    @test q.v[i,j] == 1.0 && q.v[i,j+1] == -1.0
    q.v[i,j] = q.v[i,j+1] = 0.0
    @test iszero(q.v)
  end

  @testset "Face data curl" begin
    cellcurl = curl(facexunit)
    @test cellcurl[i,j] == -1.0 && cellcurl[i,j+1] == 1.0
    cellcurl[i,j] = cellcurl[i,j+1] = 0.0
    @test iszero(cellcurl)
    cellcurl = curl(faceyunit)
    @test cellcurl[i,j] == 1.0 && cellcurl[i+1,j] == -1.0
    cellcurl[i,j] = cellcurl[i+1,j] = 0.0
    @test iszero(cellcurl)
  end

  @testset "Face data divergence" begin
    nodediv = divergence(facexunit)
    @test nodediv[i,j] == -1.0 && nodediv[i-1,j] == 1.0
    nodediv[i,j] = nodediv[i-1,j] = 0.0
    @test iszero(nodediv)
    nodediv = divergence(faceyunit)
    @test nodediv[i,j] == -1.0 && nodediv[i,j-1] == 1.0
    nodediv[i,j] = nodediv[i,j-1] = 0.0
    @test iszero(nodediv)
  end

  @testset "Face data Laplacian" begin
    lap = laplacian(facexunit)
    @test lap.u[i,j] == -4.0
    lap.u[i,j] = 0.0
    @test lap.u[i+1,j] == lap.u[i-1,j] == lap.u[i,j-1] == lap.u[i,j+1] == 1.0
    lap.u[i+1,j] = lap.u[i-1,j] = lap.u[i,j-1] = lap.u[i,j+1] = 0.0
    @test iszero(lap.u)
    @test iszero(lap.v)

    lap = laplacian(faceyunit)
    @test lap.v[i,j] == -4.0
    lap.v[i,j] = 0.0
    @test lap.v[i+1,j] == lap.v[i-1,j] == lap.v[i,j-1] == lap.v[i,j+1] == 1.0
    lap.v[i+1,j] = lap.v[i-1,j] = lap.v[i,j-1] = lap.v[i,j+1] = 0.0
    @test iszero(lap.u)
    @test iszero(lap.v)
  end

  @testset "Dual face data divergence" begin
    celldiv = divergence(dualfacexunit)
    @test celldiv[i,j] == 1.0 && celldiv[i+1,j] == -1.0
    celldiv[i,j] = celldiv[i+1,j] = 0.0
    @test iszero(celldiv)
    celldiv = divergence(dualfaceyunit)
    @test celldiv[i,j] == 1.0 && celldiv[i,j+1] == -1.0
    celldiv[i,j] = celldiv[i,j+1] = 0.0
    @test iszero(celldiv)
  end

  @testset "Face data shift to dual face" begin
    shiftx = Edges(Dual,cellzero)
    grid_interpolate!(shiftx,facexunit)
    @test shiftx.u[i,j] == shiftx.u[i-1,j] == shiftx.u[i,j+1] == shiftx.u[i-1,j+1] == 0.25
    shiftx.u[i,j] = shiftx.u[i-1,j] = shiftx.u[i,j+1] = shiftx.u[i-1,j+1] = 0.0
    @test iszero(shiftx.u)
    @test iszero(shiftx.v)
    shifty = Edges(Dual,cellzero)
    grid_interpolate!(shifty,faceyunit)
    @test shifty.v[i,j] == shifty.v[i,j-1] == shifty.v[i+1,j] == shifty.v[i+1,j-1] == 0.25
    shifty.v[i,j] = shifty.v[i,j-1] = shifty.v[i+1,j] = shifty.v[i+1,j-1] = 0.0
    @test iszero(shifty.u)
    @test iszero(shifty.v)
  end

  @testset "Dual cell center data shift to dual face" begin
    w = Edges(Dual,cellzero)
    grid_interpolate!(w,cellunit)
    @test w.u[i,j] == w.u[i-1,j] == 0.5
    w.u[i,j] = w.u[i-1,j] = 0.0
    @test iszero(w.u)
    @test w.v[i,j] == w.v[i,j-1] == 0.5
    w.v[i,j] = w.v[i,j-1] = 0.0
    @test iszero(w.v)
  end

  @testset "Face data shift to dual cell center" begin
    cellx = Nodes(Dual,cellzero)
    celly = Nodes(Dual,cellzero)
    grid_interpolate!((cellx,celly),facexunit)
    @test cellx[i,j] == 0.5 && cellx[i,j+1] == 0.5
    cellx[i,j] = cellx[i,j+1] = 0.0
    @test iszero(cellx)
    @test iszero(celly)

    celltmp = Nodes(Dual,cellzero)
    grid_interpolate!(celltmp,facexunit)
    @test celltmp[i,j] == 0.5 && celltmp[i,j+1] == 0.5
    celltmp[i,j] = celltmp[i,j+1] = 0.0
    @test iszero(celltmp)

    cellx = Nodes(Dual,cellzero)
    celly = Nodes(Dual,cellzero)
    grid_interpolate!((cellx,celly),faceyunit)
    @test celly[i,j] == 0.5 && celly[i+1,j] == 0.5
    celly[i,j] = celly[i+1,j] = 0.0
    @test iszero(cellx)
    @test iszero(celly)

    celltmp = Nodes(Dual,cellzero)
    grid_interpolate!(celltmp,faceyunit)
    @test celltmp[i,j] == 0.5 && celltmp[i+1,j] == 0.5
    celltmp[i,j] = celltmp[i+1,j] = 0.0
    @test iszero(celltmp)

  end

  @testset "Nonlinear operations" begin
    usq = Nodes(Primal,facexunit)
    magsq!(usq,facexunit)
    @test usq[i,j] == 0.5
    @test magsq(facexunit) == usq

    magsq!(usq,faceyunit)
    @test usq[i,j] == 0.5

    usq = Nodes(Dual,facexunit)
    magsq!(usq,dualfacexunit)
    @test usq[i,j] == 0.5

    magsq!(usq,dualfaceyunit)
    @test usq[i,j] == 0.5

    magu = Nodes(Primal,facexunit)
    mag!(magu,facexunit)
    @test magu[i,j] == sqrt(0.5)
    @test mag(facexunit) == magu

    magu = Nodes(Dual,facexunit)
    mag!(magu,dualfacexunit)
    @test magu[i,j] == sqrt(0.5)
    @test mag(dualfacexunit) == magu

  end

  @testset "div curl" begin
    @test iszero(divergence(curl(cellunit)))
    @test iszero(curl(grad(nodeunit)))
  end

  function randomize!(f::ScalarGridData;offset=0)
      f[offset:end-offset,offset:end-offset] .= randn(size(f[offset:end-offset,offset:end-offset]))
      return f
  end

  v = Edges(Primal,(100,200))
  offset = 4
  randomize!(v.u,offset=offset)
  randomize!(v.v,offset=offset)

  u = Edges(Primal,(100,200))
  randomize!(u.u,offset=offset)
  randomize!(u.v,offset=offset)

  vd = Edges(Dual,(100,200))
  randomize!(v.u,offset=offset)
  randomize!(v.v,offset=offset)

  ud = Edges(Dual,(100,200))
  randomize!(u.u,offset=offset)
  randomize!(u.v,offset=offset)

  f1 = Nodes(Primal,u)
  f2 = Nodes(Primal,u)
  randomize!(f1,offset=offset)
  randomize!(f2,offset=offset)

  w = Nodes(Dual,u)
  randomize!(w,offset=offset)

  dq = EdgeGradient(Primal,v)
  dv = EdgeGradient(Primal,v)
  dqd = EdgeGradient(Dual,v)
  dvd = EdgeGradient(Dual,v)

  lhsn = Nodes(Primal,f1)
  rhsn = Nodes(Primal,f1)
  lhsnd = Nodes(Dual,f1)
  rhsnd = Nodes(Dual,f1)

  lhse = Edges(Primal,f1)
  rhse = Edges(Primal,f1)
  lhsed = Edges(Dual,f1)

  qtmp = Edges(Primal,u)
  qtmpd = Edges(Dual,u)
  qtmp1d = Edges(Dual,u)
  qtmp2d = Edges(Dual,u)
  ftmp = Nodes(Primal,u)
  wtmp = Nodes(Dual,u)
  utmp = EdgeGradient(Primal,u)
  utmpd = EdgeGradient(Dual,u)

  @testset "Commutativity of operators" begin

    # Laplacian as divergence of the gradient (edge data)
    divergence!(lhse,grad!(dv,v))
    rhse = laplacian(v)
    @test isapprox(norm(lhse-rhse),0.0,atol=100.0*eps())

    # Node laplacian of the edge-to-node interpolation =
    # edge-to-node interpolation of the edge laplacian
    lhsn .= laplacian(grid_interpolate!(ftmp,v))
    rhsn .= 0.0
    grid_interpolate!(rhsn,laplacian(v))
    @test isapprox(norm(lhsn-rhsn),0.0,atol=100.0*eps())

    # Edge laplacian of the node-to-edge interpolation =
    # node-to-edge interpolation of the node laplacian
    lhse .= 0
    rhse .= 0
    qtmp .= 0
    grid_interpolate!(lhse,laplacian(f1))
    rhse .= laplacian(grid_interpolate!(qtmp,f1))
    @test isapprox(norm(lhse-rhse),0.0,atol=100.0*eps())

    # Curl of the edge laplacian = node laplacian of the curl
    lhsnd .= curl(laplacian(v))
    rhsnd .= laplacian(curl(v))
    @test isapprox(norm(lhsnd-rhsnd),0.0,atol=100.0*eps())

  end

  @testset "Discrete product rules" begin

    # Laplacian of curl = curl of Laplacian
    curllapv = curl(laplacian(v))
    lapcurlv = laplacian(curl(v))

    @test isapprox(norm(curllapv-lapcurlv),0.0,atol=100.0*eps())

    # grad(f1*f2) = f1*grad(f2) + grad(f1)*f2
    lhse .= 0
    rhse .= 0
    qtmp .= 0
    lhse .= grad(f1∘f2)
    rhse .= grad(f1)∘ grid_interpolate!(qtmp,f2) + grid_interpolate!(qtmp,f1) ∘ grad(f2)
    @test isapprox(norm(lhse-rhse),0.0,atol=100.0*eps())

    # div(f1*u) = grad(f1).u + f1*(div u)
    lhsn .= 0
    rhsn .= 0
    qtmp .= 0
    lhsn .= divergence(grid_interpolate!(qtmp,f1)∘u)
    grid_interpolate!(ftmp,grad(f1)∘u)
    rhsn .= ftmp + f1∘divergence(u)
    @test isapprox(norm(lhsn-rhsn),0.0,atol=100.0*eps())

    # div(u*w) = grad(w).u + w*(div u)
    lhsnd .= 0
    rhsnd .= 0
    qtmp1d .= 0
    qtmp2d .= 0
    lhsnd .= divergence(grid_interpolate!(qtmp1d,u)∘grid_interpolate!(qtmp2d,w))
    qtmp1d .= 0
    grid_interpolate!(rhsnd,grid_interpolate!(qtmp1d,u)∘grad(w))
    wtmp .= 0
    rhsnd .+= w∘grid_interpolate!(wtmp,divergence(u))
    @test isapprox(norm(lhsnd-rhsnd),0.0,atol=100.0*eps())

    # div(uv) = u.grad v + (div u)v
    divuv = zero(u)
    uv = EdgeGradient(Primal,u)
    divergence!(divuv,u*v)

    dv .= 0.0
    utmp .= 0.0
    dq .= 0
    ugradv = zero(u)
    product!(dq,transpose(grid_interpolate!(utmp,u)),grad!(dv,v))
    grid_interpolate!(ugradv,dq)

    lhse .= 0
    directional_derivative!(lhse,v,u)
    @test isapprox(norm(lhse - ugradv),0.0,atol=100.0*eps())

    qtmp .= 0.0
    vdivu = zero(u)
    product!(vdivu,grid_interpolate!(qtmp,divergence(u)),v)

    @test isapprox(norm(divuv - ugradv - vdivu),0.0,atol=100.0*eps())

    # Now on duals
    divuv = zero(ud)
    uv = EdgeGradient(Primal,ud)
    divergence!(divuv,ud*vd)

    dvd .= 0.0
    utmpd .= 0.0
    dqd .= 0.0
    ugradv = zero(ud)
    product!(dqd,transpose(grid_interpolate!(utmpd,ud)),grad!(dvd,vd))
    grid_interpolate!(ugradv,dqd)

    lhsed .= 0
    directional_derivative!(lhsed,vd,ud)
    @test isapprox(norm(lhsed - ugradv),0.0,atol=100.0*eps())

    qtmpd .= 0.0
    vdivu = zero(ud)
    product!(vdivu,grid_interpolate!(qtmpd,divergence(ud)),vd)

    @test isapprox(norm(divuv - ugradv - vdivu),0.0,atol=100.0*eps())

  end

  L = plan_laplacian(nx,ny;with_inverse=true)

  @testset "Laplacian of the LGF" begin
    ψ = L\cellunit
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    ψ = L\nodeunit
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    ψ = L\facexunit.u
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    ψ = L\faceyunit.v
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    ψ = L\dualfacexunit.u
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    ψ = L\dualfaceyunit.v
    lapψ = L*ψ
    @test lapψ[i,j]≈1.0
    @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=10.0*eps()) &&
            isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=10.0*eps())

    L2 = plan_laplacian(nodeunit,with_inverse=true)
    ψ = L2\nodeunit
    lapψ = L2*ψ
    @test lapψ[i,j]≈1.0

  end

  Lscale = plan_laplacian(nx,ny;with_inverse=true,factor=2.0)

  @testset "Laplacian of the LGF with factor" begin
     ψ = Lscale\cellunit
     lapψ = Lscale*ψ
     @test lapψ[i,j]≈1.0

     ψ = Lscale*cellunit
     @test ψ[i,j]≈-8.0 && ψ[i+1,j]≈2.0

     mul!(ψ,Lscale,cellunit)
     @test ψ[i,j]≈-8.0 && ψ[i-1,j]≈2.0

  end

  @testset "LGF for Helmholtz equation" begin
    alpha = 0.02
    LH(i,j,f::Function,α) = im*α*f(i,j,α)-(f(i-1,j,α)+f(i+1,j,α)+f(i,j+1,α)+f(i,j-1,α)-4*f(i,j,α))

    i0, j0 = rand(1:100), rand(0:100)
    @test abs(LH(i0,j0,CartesianGrids.lgf_helmholtz,alpha)) < 100.0*eps()

    @test isapprox(real(LH(0,0,CartesianGrids.lgf_helmholtz,alpha)),1.0;atol=100.0*eps())

  end

end

##### COMPLEX GRID #####

@testset "Complex Grid Routines" begin

# size
nx = 12; ny = 12

# sample point
i = 5; j = 7

cellzero = Nodes(Dual,(nx,ny),dtype=ComplexF64)
nodezero = Nodes(Primal,cellzero)
facezero = Edges(Primal,cellzero)
dualfacezero = Edges(Dual,cellzero)

a = 1.0+2.0im

cellunit = deepcopy(cellzero)
cellunit[i,j] = a

nodeunit = deepcopy(nodezero)
nodeunit[i,j] = a

facexunit = deepcopy(facezero)
facexunit.u[i,j] = a

faceyunit = deepcopy(facezero)
faceyunit.v[i,j] = a

dualfacexunit = deepcopy(dualfacezero)
dualfacexunit.u[i,j] = a

dualfaceyunit = deepcopy(dualfacezero)
dualfaceyunit.v[i,j] = a

@testset "Basic array operations" begin
  w = zero(cellunit)
  w .= cellunit
  @test w[i,j] == a
  q = similar(facexunit)
  q .= facexunit
  @test q.u[i,j] == a
  @test iszero(q.v.data)
end

@testset "Basic complex operations" begin
  w = conj(cellunit)
  @test typeof(w) == typeof(cellunit)
  @test w[i,j] == conj(a)

  w = real(cellunit)
  @test typeof(w) <: Nodes{Dual,nx,ny,Float64}
  @test real(w[i,j]) == real(a) && imag(w[i,j]) == 0.0

  w = imag(cellunit)
  @test typeof(w) <: Nodes{Dual,nx,ny,Float64}
  @test real(w[i,j]) == imag(a) && imag(w[i,j]) == 0.0

end

@testset "Inner products and norms" begin
  w = zero(cellunit)
  i0, j0 = rand(2:nx-1), rand(2:ny-1)
  w[i0,j0] = 1.0im
  @test norm(w)*sqrt((nx-2)*(ny-2)) == 1.0
  w .= 1.0im
  @test norm(w) == 1.0

  @test 2*w == w*2

  p = Nodes(Primal,w)
  p .= 1.0im
  @test norm(p) == 1.0
  p2 = deepcopy(p)
  @test dot(p,p2) == 1.0
  @test norm(p-p2) == 0.0

  q = Edges(Dual,w)
  q.u .= 1.0im
  q2 = deepcopy(q)
  @test dot(q,q2) == 1.0

  q = Edges(Primal,w)
  q.u .= 1.0im
  q2 = deepcopy(q)
  @test dot(q,q2) == 1.0

  @test integrate(w) == 1.0im

  @test integrate(p) == 1.0im

  q .= 1.0im
  @test norm(2*q) == sqrt(8)

end

@testset "Dual cell center data Laplacian" begin
  lapcell = laplacian(cellunit)
  @test lapcell[i,j] == -4.0*a
  lapcell[i,j] = 0.0
  @test lapcell[i+1,j] == lapcell[i-1,j] == lapcell[i,j-1] == lapcell[i,j+1] == a
  lapcell[i+1,j] = lapcell[i-1,j] = lapcell[i,j-1] = lapcell[i,j+1] = 0.0
  @test iszero(lapcell)
end

@testset "Dual cell center data curl" begin
  q = curl(cellunit)
  @test q.u[i,j-1] == 1.0*a && q.u[i,j] == -1.0*a
  q.u[i,j-1] = q.u[i,j] = 0.0
  @test iszero(q.u)
  @test q.v[i-1,j] == -1.0*a && q.v[i,j] == 1.0*a
  q.v[i-1,j] = q.v[i,j] = 0.0
  @test iszero(q.v)
end

@testset "Dual cell node gradient" begin
  q = grad(nodeunit)
  @test q.u[i,j] == 1.0*a && q.u[i+1,j] == -1.0*a
  q.u[i,j] = q.u[i+1,j] = 0.0
  @test iszero(q.u)
  @test q.v[i,j] == 1.0*a && q.v[i,j+1] == -1.0*a
  q.v[i,j] = q.v[i,j+1] = 0.0
  @test iszero(q.v)
end

@testset "Face data curl" begin
  cellcurl = curl(facexunit)
  @test cellcurl[i,j] == -1.0*a && cellcurl[i,j+1] == 1.0*a
  cellcurl[i,j] = cellcurl[i,j+1] = 0.0
  @test iszero(cellcurl)
  cellcurl = curl(faceyunit)
  @test cellcurl[i,j] == 1.0*a && cellcurl[i+1,j] == -1.0*a
  cellcurl[i,j] = cellcurl[i+1,j] = 0.0
  @test iszero(cellcurl)
end

@testset "Face data divergence" begin
  nodediv = divergence(facexunit)
  @test nodediv[i,j] == -1.0*a && nodediv[i-1,j] == 1.0*a
  nodediv[i,j] = nodediv[i-1,j] = 0.0
  @test iszero(nodediv)
  nodediv = divergence(faceyunit)
  @test nodediv[i,j] == -1.0*a && nodediv[i,j-1] == 1.0*a
  nodediv[i,j] = nodediv[i,j-1] = 0.0
  @test iszero(nodediv)
end

@testset "Face data Laplacian" begin
  lap = laplacian(facexunit)
  @test lap.u[i,j] == -4.0*a
  lap.u[i,j] = 0.0
  @test lap.u[i+1,j] == lap.u[i-1,j] == lap.u[i,j-1] == lap.u[i,j+1] == 1.0*a
  lap.u[i+1,j] = lap.u[i-1,j] = lap.u[i,j-1] = lap.u[i,j+1] = 0.0
  @test iszero(lap.u)
  @test iszero(lap.v)

  lap = laplacian(faceyunit)
  @test lap.v[i,j] == -4.0*a
  lap.v[i,j] = 0.0
  @test lap.v[i+1,j] == lap.v[i-1,j] == lap.v[i,j-1] == lap.v[i,j+1] == 1.0*a
  lap.v[i+1,j] = lap.v[i-1,j] = lap.v[i,j-1] = lap.v[i,j+1] = 0.0
  @test iszero(lap.u)
  @test iszero(lap.v)
end

@testset "Dual face data divergence" begin
  celldiv = divergence(dualfacexunit)
  @test celldiv[i,j] == 1.0*a && celldiv[i+1,j] == -1.0*a
  celldiv[i,j] = celldiv[i+1,j] = 0.0
  @test iszero(celldiv)
  celldiv = divergence(dualfaceyunit)
  @test celldiv[i,j] == 1.0*a && celldiv[i,j+1] == -1.0*a
  celldiv[i,j] = celldiv[i,j+1] = 0.0
  @test iszero(celldiv)
end

@testset "Face data shift to dual face" begin
  shiftx = Edges(Dual,cellzero)
  grid_interpolate!(shiftx,facexunit)
  @test shiftx.u[i,j] == shiftx.u[i-1,j] == shiftx.u[i,j+1] == shiftx.u[i-1,j+1] == 0.25*a
  shiftx.u[i,j] = shiftx.u[i-1,j] = shiftx.u[i,j+1] = shiftx.u[i-1,j+1] = 0.0
  @test iszero(shiftx.u)
  @test iszero(shiftx.v)
  shifty = Edges(Dual,cellzero)
  grid_interpolate!(shifty,faceyunit)
  @test shifty.v[i,j] == shifty.v[i,j-1] == shifty.v[i+1,j] == shifty.v[i+1,j-1] == 0.25*a
  shifty.v[i,j] = shifty.v[i,j-1] = shifty.v[i+1,j] = shifty.v[i+1,j-1] = 0.0
  @test iszero(shifty.u)
  @test iszero(shifty.v)
end

@testset "Dual cell center data shift to dual face" begin
  w = Edges(Dual,cellzero)
  grid_interpolate!(w,cellunit)
  @test w.u[i,j] == w.u[i-1,j] == 0.5*a
  w.u[i,j] = w.u[i-1,j] = 0.0
  @test iszero(w.u)
  @test w.v[i,j] == w.v[i,j-1] == 0.5*a
  w.v[i,j] = w.v[i,j-1] = 0.0
  @test iszero(w.v)
end

@testset "Face data shift to dual cell center" begin
  cellx = Nodes(Dual,cellzero)
  celly = Nodes(Dual,cellzero)
  grid_interpolate!((cellx,celly),facexunit)
  @test cellx[i,j] == 0.5*a && cellx[i,j+1] == 0.5*a
  cellx[i,j] = cellx[i,j+1] = 0.0
  @test iszero(cellx)
  @test iszero(celly)

  cellx = Nodes(Dual,cellzero)
  celly = Nodes(Dual,cellzero)
  grid_interpolate!((cellx,celly),faceyunit)
  @test celly[i,j] == 0.5*a && celly[i+1,j] == 0.5*a
  celly[i,j] = celly[i+1,j] = 0.0
  @test iszero(cellx)
  @test iszero(celly)
end

@testset "div curl" begin
  @test iszero(divergence(curl(cellunit)))
  @test iszero(curl(grad(nodeunit)))
end

L = plan_laplacian(nx,ny;with_inverse=true,dtype=ComplexF64)

@testset "Laplacian of the LGF" begin
  ψ = L\cellunit
  lapψ = L*ψ
  @test lapψ[i,j]≈1.0*a
  @test isapprox(maximum(abs.(lapψ[Not(i),:])),0.0;atol=100.0*eps()) &&
          isapprox(maximum(abs.(lapψ[:,Not(j)])),0.0;atol=100.0*eps())
end

Lscale = plan_laplacian(nx,ny;with_inverse=true,factor=2.0,dtype=ComplexF64)

@testset "Laplacian of the LGF with factor" begin
   ψ = Lscale\cellunit
   lapψ = Lscale*ψ
   @test lapψ[i,j]≈1.0*a

   ψ = Lscale*cellunit
   @test ψ[i,j]≈-8.0*a && ψ[i+1,j]≈2.0*a

   mul!(ψ,Lscale,cellunit)
   @test ψ[i,j]≈-8.0*a && ψ[i-1,j]≈2.0*a

end

α = 0.07
LH = plan_helmholtz(nx,ny,α;with_inverse=true)

@testset "Helmholtz of the LGF" begin
  ψ = LH\cellunit
  helmψ = LH*ψ
  @test helmψ[i,j]≈1.0*a
  @test isapprox(maximum(abs.(helmψ[Not(i),:])),0.0;atol=100.0*eps()) &&
          isapprox(maximum(abs.(helmψ[:,Not(j)])),0.0;atol=100.0*eps())
end

end

@testset "Fields" begin
    @testset "Hadamard Product" begin
        edges_p  = Edges(Primal,(30,40))
        edges_p.u .= rand(Float64,size(edges_p.u))

        # Should be safe for the output to be the same as the input
        edges_p2 = edges_p ∘ edges_p
        product!(edges_p, edges_p, edges_p)
        @test edges_p2.u == edges_p.u
        @test edges_p2.v == edges_p.v

        edges_d  = Edges{Dual, 30, 40, Float64}()
        @test_throws MethodError (edges_p ∘ edges_d)
    end

    @testset "Cartesian product of vectors" begin

      edges_p  = Edges(Primal,(30,40))
      edges_p.u .= rand(Float64,size(edges_p.u))

      edges_q  = Edges(Primal,(30,40))
      edges_q.u .= rand(Float64,size(edges_q.u))

      tensors_pq = edges_p * edges_q

      pt = EdgeGradient(Primal,Dual,(30,40))
      qt = EdgeGradient(Primal,Dual,(30,40))

      grid_interpolate!(pt,edges_p)
      grid_interpolate!(qt,edges_q)

      @test tensors_pq.dudx[12,25] == pt.dudx[12,25]*qt.dudx[12,25]
      @test tensors_pq.dvdx[12,25] == pt.dudy[12,25]*qt.dvdx[12,25]
      @test tensors_pq.dudy[12,25] == pt.dvdx[12,25]*qt.dudy[12,25]
      @test tensors_pq.dvdy[12,25] == pt.dvdy[12,25]*qt.dvdy[12,25]

    end

    @testset "Discrete Laplacian" begin
        s = Nodes(Dual,(30,40))
        s[3:end-2, 3:end-2] .= rand(26, 36)

        L = plan_laplacian(30, 40)

        @test L*s ≈ -curl(curl(s))

        @test_throws MethodError (L \ s)

        L = plan_laplacian(30, 40, with_inverse = true, fftw_flags = FFTW.PATIENT)
        @test L \ (L*s) ≈ s

        L! = plan_laplacian!(30, 40, with_inverse = true, fftw_flags = FFTW.PATIENT)
        sold = deepcopy(s)
        L! \ (L!*s)
        @test s ≈ sold
    end

    @testset "Integrating factor" begin
        s = Nodes(Dual,(30,40))
        s[15,15] = 1.0

        E1 = plan_intfact(1,s)
        E2 = plan_intfact(2,s)

        @test E1*(E1*s) ≈ E2*s

        E! = plan_intfact!(2,s)
        s2 = deepcopy(s)
        E! * s

        @test s ≈ E2*s2

        L = plan_laplacian(s,factor=2)
        EL = exp(L,1)
        @test EL*s ≈ E2*s

        p = XEdges(Primal,s)
        p[15,15] = 1.0
        @test E1*(E1*p) ≈ E2*p

        q = Edges(Primal,s)
        q.u[15,15] = 1.0
        @test E1*(E1*q) ≈ E2*q

        EmL = exp(L,-1)
        @test EmL\s == EL*s

    end

    @testset "Discrete Divergence" begin
        s = Nodes(Dual,(5,4))
        s .= rand(5, 4)

        @test norm(divergence(curl(s))) ≈ 0 atol=eps()

        s = Nodes(Primal,s)
        q′ = Edges(Primal,s)
        q′.u .= reshape(1:15, 5, 3)
        q′.v .= reshape(1:16, 4, 4)

        # Not sure if this is the behavior we want yet
        # Currently, the ghost cells are not affected
        # by the divergence operator
        #s .= 1.0
        divergence!(s, q′)
        @test s == [ 5.0  5.0  5.0
                     5.0  5.0  5.0
                     5.0  5.0  5.0
                     5.0  5.0  5.0 ]
    end

    @testset "Discrete Curl" begin
        s = Nodes(Dual,(5,4))
        s .= reshape(1:20, 4, 5)'

        q = curl(s)

        @test q.u == [ 1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0
                       1.0  1.0  1.0 ]

        @test q.v == [ -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0
                       -4.0  -4.0  -4.0  -4.0 ]
    end

    @testset "Shifting Primal Edges to Dual Edges" begin

        q = Edges(Primal,(5,4))
        q.u .= reshape(1:15, 5, 3)
        q.v .= reshape(1:16, 4, 4)
        Qq = Edges(Dual,q)

        grid_interpolate!(Qq,q)
        @test Qq.u == [ 0.0  4.0  9.0  0.0
                        0.0  5.0  10.0 0.0
                        0.0  6.0  11.0 0.0
                        0.0  7.0  12.0 0.0 ]

        @test Qq.v == [ 0.0  0.0   0.0
                        3.5  7.5  11.5
                        4.5  8.5  12.5
                        5.5  9.5  13.5
                        0.0  0.0   0.0 ]

    end

    @testset "Shifting Dual Edges to Primal Edges" begin

        q = Edges(Dual,(5,4))
        q.u .= reshape(1:16, 4, 4)
        q.v .= reshape(1:15, 5, 3)
        v = Edges(Primal,q)
        grid_interpolate!(v,q)

        @test v.u == [ 0.0  0.0   0.0
                       3.5  7.5  11.5
                       4.5  8.5  12.5
                       5.5  9.5  13.5
                       0.0  0.0   0.0 ]

        @test v.v == [ 0.0  4.0   9.0  0.0
                        0.0  5.0  10.0  0.0
                        0.0  6.0  11.0  0.0
                        0.0  7.0  12.0  0.0 ]

    end

    @testset "Shifting Dual Nodes to Dual Edges" begin

        w = Nodes(Dual,(5,4))
        w .= reshape(1:20, 5, 4)

        Ww = Edges(Dual,w)
        grid_interpolate!(Ww,w)

        @test Ww.u == [ 0.0  6.5  11.5  0.0
                        0.0  7.5  12.5  0.0
                        0.0  8.5  13.5  0.0
                        0.0  9.5  14.5  0.0 ]
        @test Ww.v == [ 0.0   0.0   0.0
                        4.5   9.5  14.5
                        5.5  10.5  15.5
                        6.5  11.5  16.5
                        0.0   0.0   0.0 ]
    end

    @testset "Physical grid" begin

        g = PhysicalGrid((-1.0,3.0),(-2.0,3.0),0.02)
        @test size(g) == (208,252)
        @test size(g,1) == 208
        @test size(g,2) == 252
        @test length(g) == 208*252
        @test origin(g) == (54,101)
        @test cellsize(g) == 0.02
        @test limits(g,1) == (-1.06,3.06)
        @test limits(g,2) == (-2.0,3.0)

    end
end


text:
using Test, FileIO
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

# Load in input data
rp = joinpath(dirname(@__FILE__), "../reference/SLM/")
in_inc = load(joinpath(rp, "solve_slm_system.jld2"))
in_dec = load(joinpath(rp, "solve_slm_system_decrease.jld2"))
in_sin = load(joinpath(rp, "solve_slm_system_sine.jld2"))
in_sinnoineq = load(joinpath(rp, "solve_slm_system_sinenoineq.jld2"))
noconstraints = load(joinpath(rp, "solve_slm_system_noconstraints.jld2"))

# Run tests
@testset "Coefficients of least-square spline" begin
    for in_data in [in_inc, in_dec, in_sin, in_sinnoineq]
        @test maximum(abs.(in_data["coef"] -
                           solve_slm_system(in_data["Mdes"], vec(in_data["rhs"]), in_data["Mreg"],
                                            vec(in_data["rhsreg"]), in_data["finalRP"], in_data["Meq"], vec(in_data["rhseq"]),
                                            in_data["Mineq"], vec(in_data["rhsineq"])))) < 5e-3
    end

    # Check no constraints doesn't error
    @test noconstraints["coefs"] ≈ solve_slm_system(in_sinnoineq["Mdes"], vec(in_sinnoineq["rhs"]), in_sinnoineq["Mreg"],
                                                    vec(in_sinnoineq["rhsreg"]), in_sinnoineq["finalRP"],
                                                    Matrix{Float64}(undef, 0, 0), Vector{Float64}(undef, 0),
                                                    in_sinnoineq["Mineq"], vec(in_sinnoineq["rhsineq"]))
end

#=
println("No sparse matrices")
@btime begin
    solve_slm_system(in_inc["Mdes"], vec(in_inc["rhs"]), in_inc["Mreg"],
                     vec(in_inc["rhsreg"]), in_inc["finalRP"], in_inc["Meq"], vec(in_inc["rhseq"]),
                     in_inc["Mineq"], vec(in_inc["rhsineq"]))
end

println("Sparse matrices")
@btime begin
    solve_slm_system(in_inc["Mdes"], vec(in_inc["rhs"]), in_inc["Mreg"],
                     vec(in_inc["rhsreg"]), in_inc["finalRP"], in_inc["Meq"], vec(in_inc["rhseq"]),
                     in_inc["Mineq"], vec(in_inc["rhsineq"]); use_sparse = true)
end
println("Adding an initial guess for coefficients")
@btime begin
    solve_slm_system(in_inc["Mdes"], vec(in_inc["rhs"]), in_inc["Mreg"],
                     vec(in_inc["rhsreg"]), in_inc["finalRP"], in_inc["Meq"], vec(in_inc["rhseq"]),
                     in_inc["Mineq"], vec(in_inc["rhsineq"]); use_sparse = true,
                     init = vec(in_inc["coef"]))
end

println("Use LLS")
@btime begin
    solve_slm_system(in_inc["Mdes"], vec(in_inc["rhs"]), in_inc["Mreg"],
                     vec(in_inc["rhsreg"]), in_inc["finalRP"], in_inc["Meq"], vec(in_inc["rhseq"]),
                     in_inc["Mineq"], vec(in_inc["rhsineq"]))
end

println("Use ADNLS")
@btime begin
    solve_slm_system(in_inc["Mdes"], vec(in_inc["rhs"]), in_inc["Mreg"],
                     vec(in_inc["rhsreg"]), in_inc["finalRP"], in_inc["Meq"], vec(in_inc["rhseq"]),
                     in_inc["Mineq"], vec(in_inc["rhsineq"]); use_lls = false)
end
=#


text:
export auto_init, auto_logenergy, logposterior
import OmegaCore

# # FIXME: Remoev soft predicates from here, but how?!?
# using SoftPredicates

"Randomly initialize `ω::ΩT`"
function auto_init(rng, x, ΩT)
  ω = ΩT()
  ω_ = OmegaCore.tagignorecondition(OmegaCore.tagrng(ω, rng))
  x(ω_)
  ω
end

"Exo prior soft constraints.  `T` is type of energy"
function auto_logenergy(x, ::Type{T}) where T
  function logenergyvar_(ω)
    ϵ = OmegaCore.condvarapply(x, ω, T) # need this in here for type stability - Julia limitation
    ℓ = OmegaCore.logenergyexo(ω)
    SoftPredicates.logerr(ϵ) + ℓ
  end
end

# le(x::Real) = x
# le(x::SoftPredicates.AbstractSoftBool) = SoftPredicates.logerr(x)

"Random variable that computes `logposterior(ω)` given `loglikelihood`"
function logposterior(likelihood)
  function logpost_(ω)
    ϵ = likelihood(ω)
    ℓ = OmegaCore.logenergyexo(ω)
    ϵ + ℓ
  end
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# include -
include("Include.jl")

# What index do I want to see?
plot_index = 41;

# Load the simulation result set -
result_set = process_simulation_files("./simple_ensemble")
for (index,data_array) in enumerate(result_set)

  @show index

  time_array = data_array[:,1];
  state_array = data_array[:,plot_index+1];

  plot(time_array,state_array);
end


text:
using NMRTools
using SimpleTraits
using Test

# generate some sample data
x1 = 0:.1:1;
x2 = 1:4;
x3 = 10:10:30;
y1 = [cos(n1) for n1 in x1];
y2 = [cos(n1)+n2 for n1 in x1,  n2 in x2];
y3 = [cos(n1)+n2/n3 for n1 in x1,  n2 in x2, n3 in x3];
md = Dict{Any,Any}(:test=>"hello data");
md1 = Dict{Any,Any}(:test=>"hello axis 1");
md2 = Dict{Any,Any}(:test=>"hello axis 2");
md3 = Dict{Any,Any}(:test=>"hello axis 3");
@test y2[1,1] == 2
@test y2[11,4] ≈ 4.5403023

@testset "Creating axes" begin
    for f in [X, Y, Z, Ti]
        ax = f(0:.1:1)
        @test ax[1] == 0.0
        #@test val(ax) == 0:.1:1
        #@test metadata(ax) isa Dict
    end
end

@testset "Creating NMRData" begin
    d = NMRData(y2, (X(x1), Y(x2)))
    @test metadata(d) isa Dict
    d = NMRData(y2, (X(x1), Y(x2)), metadata=md)
    @test d[1,1] == 2
    @test d == y2
    @test metadata(d) isa Dict
    @test metadata(d)[:test] == "hello data"
    @test size(d) == (11,4)
    @test length(d) == 44
    @test d[At(1.0),At(4)] ≈ 4.5403023
    @test d[Near(-0.5), Near(0.9)] == 2
    @test size(d[Between(0.25,0.65),Between(3,4)]) == (4,2)
    d = NMRData(y2, (X(0:.1:1, metadata=md1), Y(1:4, metadata=md2)), metadata=md)
    @test d[1,1] == 2
    @test d == y2
    @test metadata(d,X)[:test] == "hello axis 1"
    @test metadata(d,Y)[:test] == "hello axis 2"
end

@testset "Creating 1D NMRData" begin
    d = NMRData(y1, (X(x1), ))
    @test metadata(d) isa Dict
    @test d[1] == 1
    d = NMRData(y1, (X(x1), ), metadata=md)
    @test d[1] == 1
    @test d == y1
    @test metadata(d) isa Dict
    @test metadata(d)[:test] == "hello data"
end

@traitfn f(x::X) where {X; HasPseudoDimension{X}} = "Fake!"
@traitfn f(x::X) where {X; Not{HasPseudoDimension{X}}} = "Real!"

@testset "HasPseudoDim trait" begin
    d = NMRData(y1, (X(x1), ))
    @test f(d) == "Real!"
    d = NMRData(y2, (X(x1), Y(x2)))
    @test f(d) == "Real!"
    d = NMRData(y2, (X(x1), Ti(x2)))
    @test f(d) == "Fake!"
    d = NMRData(y3, (X(x1), Y(x2), Z(x3)))
    @test f(d) == "Real!"
    d = NMRData(y3, (X(x1), Y(x2), Ti(x3)))
    @test f(d) == "Fake!"
end


text:
# Bogumił Kamiński, 2022

# Codes for chapter 14

# Codes for section 14.1

using Plots
using Statistics

X = [1.0, 1.1, 1.3, 1.2, 1.2]
T = 1.0
m = 4
Y = mean(X)
K = 1.05
plot(range(0.0, T length=m+1), X;
           xlabel="T", legend=false, color="black")
hline!([Y], color="gray", lw=3, ls=:dash)
hline!([K], color="gray", lw=3, ls=:dot)
annotate!([(T, Y + 0.01, "Y"),
           (T, K + 0.01, "K"),
           (T, X[end] + 0.01, "X")])

using DataFrames
using Random
Random.seed!(1234);
X0, T, s, r, m = 1.0, 2.0, 0.2, 0.1, 4
gbm = DataFrame(X=X0, t=0.0)
for i in 1:m
    Z = randn()
    log_return = (r - s^2/2) * T/m + s * sqrt(T/m) * Z
    next_X = gbm.X[end] * exp(log_return)
    next_t = gbm.t[end] + T/m
    push!(gbm, (next_X, next_t))
end
gbm

# Codes for section 14.2

# start Julia with and additional -t4 command line switch

Threads.nthreads()

# Code for listing 14.1

function payoff_asian_sample(T, X0, K, r, s, m)::Float64
    X = X0
    sumX = X
    d = T / m
    for i in 1:m
        X *= exp((r - s^2 / 2) * d + s * sqrt(d) * randn())
        sumX += X
    end
    Y = sumX / (m + 1)
    return exp(-r * T) * max(Y - K, 0)
end

# Code for checking the results of the payoff simulation

payoff_asian_sample(1.0, 50.0, 55.0, 0.05, 0.3, 200)
payoff_asian_sample(1.0, 50.0, 55.0, 0.05, 0.3, 200)
payoff_asian_sample(1.0, 50.0, 55.0, 0.05, 0.3, 200)

# Benchmarking map

using BenchmarkTools
@btime map(i -> payoff_asian_sample(1.0, 50.0, 55.0, 0.05, 0.3, 200), 1:10_000);

using ThreadsX
@btime ThreadsX.map(i -> payoff_asian_sample(1.0, 50.0, 55.0, 0.05, 0.3, 200), 1:10_000);

# Codes for section 14.3

# Code for listing 14.2

using Statistics

function asian_value(T, X0, K, r, s, m, max_time)
    result = Float64[]
    start_time = time()
    while time() - start_time < max_time
        append!(result, ThreadsX.map(i -> payoff_asian_sample(T, X0, K, r, s, m), 1:10_000))
    end
    n = length(result)
    mv = mean(result)
    sdv = std(result)
    lo95 = mv - 1.96 * sdv / sqrt(n)
    hi95 = mv + 1.96 * sdv / sqrt(n)
    zero = mean(==(0), result)
    return (; n, mv, lo95, hi95, zero)
end

# Code for example of the mean function

mean(x -> x ^ 2, [1, 2, 3])

eq0 = ==(0)
eq0(1)
eq0(0)

# Code for shorthand NamedTuple notation

val1 = 10
val2 = "x"
(; val1, val2)

# Testing asian_value function

@time asian_value(1.0, 50.0, 55.0, 0.05, 0.3, 200, 0.25)
@time asian_value(1.0, 50.0, 55.0, 0.05, 0.3, 200, 0.25)
@time asian_value(1.0, 50.0, 55.0, 0.05, 0.3, 200, 0.25)

# Converting NamedTuple to JSON response

using Genie
Genie.Renderer.Json.json((firstname="Bogumił", lastname="Kamiński"))

# Codes for section 14.4

# Start a new Julia session

using HTTP
using JSON3
req = HTTP.post("http://127.0.0.1:8000",
                ["Content-Type" => "application/json"],
                JSON3.write((K=55.0, max_time=0.25)))
JSON3.read(req.body)

JSON3.write((K=55.0, max_time=0.25))

HTTP.post("http://127.0.0.1:8000",
          ["Content-Type" => "application/json"],
          JSON3.write((K="", max_time=0.25)))

using DataFrames
df = DataFrame(K=30:2:80, max_time=0.25)
df.data = map(df.K, df.max_time) do K, max_time
    @show K
    @time req = HTTP.post("http://127.0.0.1:8000",
                          ["Content-Type" => "application/json"],
                          JSON3.write((;K, max_time)))
    return JSON3.read(req.body)
end;

df

all(x -> x.status == "OK", df.data)

small_df = DataFrame(x=[(a=1, b=2), (a=3, b=4), (a=5, b=6)])
transform(small_df, :x => identity => AsTable)
transform(small_df, :x => AsTable)

df2 = select(df, :K, :data => ByRow(x -> x.value) => AsTable)

using Plots
plot(plot(df2.K, df2.mv; legend=false,
          xlabel="K", ylabel="expected value"),
     plot(df2.K, df2.zero; legend=false,
          xlabel="K", ylabel="probability of zero"))

max_stars_repo_name: bkamins/JuliaForDataAnalysis

max_issues_repo_name: bkamins/JuliaForDataAnalysis

max_forks_repo_name: bkamins/JuliaForDataAnalysis


text:
using UnitJuMP
using Test

import MutableArithmetics
const _MA = MutableArithmetics

@testset "Operators" begin
    include("test_operators.jl")
end

@testset "MutableArithmetics" begin
    include("test_ma.jl")
end

@testset "JuMP" begin
    include("test_jump.jl")
end