text:
using Mimi
using OptiMimi
using Test

# From http://www.purplemath.com/modules/linprog3.htm

# You need to buy some filing cabinets. You know that Cabinet X costs
# $10 per unit, requires six square feet of floor space, and holds eight
# cubic feet of files. Cabinet Y costs $20 per unit, requires eight
# square feet of floor space, and holds twelve cubic feet of files. You
# have been given $140 for this purchase, though you don't have to spend
# that much. The office has room for no more than 72 square feet of
# cabinets. How many of which model should you buy, in order to maximize
# storage volume?

# Construct the Cabinets component
@defcomp Cabinets begin
    office = Index()

    # Two kinds of cabinet: model X and Y
    x = Parameter(index=[time, office])
    y = Parameter(index=[time, office])

    cost = Variable(index=[time, office])
    space = Variable(index=[time, office])
    volume = Variable(index=[time, office])
end

# Define all gradients
grad_cabinets_cost_x(m::Model) = roomdiagonal(m, :Cabinets, :cost, :x, (tt, ii) -> 10)
grad_cabinets_cost_y(m::Model) = roomdiagonal(m, :Cabinets, :cost, :y, (tt, ii) -> 20)
grad_cabinets_space_x(m::Model) = roomdiagonal(m, :Cabinets, :space, :x, (tt, ii) -> 6)
grad_cabinets_space_y(m::Model) = roomdiagonal(m, :Cabinets, :space, :y, (tt, ii) -> 8)
grad_cabinets_volume_x(m::Model) = roomdiagonal(m, :Cabinets, :volume, :x, (tt, ii) -> 8)
grad_cabinets_volume_y(m::Model) = roomdiagonal(m, :Cabinets, :volume, :y, (tt, ii) -> 12)

# And the constraints
constraintoffset_cabinets_cost(m::Model) = hallsingle(m, :Cabinets, :cost, (tt, ii) -> 140)
constraintoffset_cabinets_space(m::Model) = hallsingle(m, :Cabinets, :space, (tt, ii) -> 72)

# Create the model
m = Model()
set_dimension!(m, :time, [1])
set_dimension!(m, :office, collect(1:2))

cabinets = add_comp!(m, Cabinets)
cabinets[:x] = zeros(1, 2)
cabinets[:y] = zeros(1, 2)

# Create the constraint matrix
paramcomps = [:Cabinets, :Cabinets]
parameters = [:x, :y]
constcomps = [:Cabinets, :Cabinets]
constraints = [:cost, :space]

house = LinearProgrammingHouse(m, paramcomps, parameters, constcomps, constraints);

setobjective!(house, varsum(grad_cabinets_volume_x(m)))
setobjective!(house, varsum(grad_cabinets_volume_y(m)))

setconstraint!(house, grad_cabinets_cost_x(m))
setconstraint!(house, grad_cabinets_cost_y(m))
setconstraintoffset!(house, constraintoffset_cabinets_cost(m))

setconstraint!(house, grad_cabinets_space_x(m))
setconstraint!(house, grad_cabinets_space_y(m))
setconstraintoffset!(house, constraintoffset_cabinets_space(m))

# Solve the problem
@time sol = houseoptimize(house)
println(sol.sol)

@test sol.sol[1] ≈ 8.0
@test sol.sol[2] ≈ 8.0
@test sol.sol[3] ≈ 3.0
@test sol.sol[4] ≈ 3.0

# 100 cubic feet by buying eight of model X and three of model Y.

# What are the constraining elements?
println(constraining(house, sol.sol))


text:
using NBodyIPs, StaticArrays, BenchmarkTools, JuLIP, Test, Profile

# using JuLIP.Potentials: evaluate, evaluate_d
# using NBodyIPs: BondLengthDesc, BondAngleDesc, invariants, invariants_d, descriptor,
#                bo2angles, bo2edges
using NBodyIPs.Polys

TRANSFORM = PolyTransform(3, 2.9) # "r -> (2.9/r)^3"
DEGREES = [18, 14, 12, 10]
RCUT = [7.0, 5.80, 4.5, 4.1]

function random_ip(DT)
   DD = [ DT(TRANSFORM, CosCut(rcut-1.0, rcut))  for rcut in RCUT ]

   # 2-body potential
   B2 = nbpolys(2, DD[1], DEGREES[1])
   c = rand(length(B2))
   V2 = NBPoly(B2, c, DD[1])
   V2sp = StNBPoly(V2)

   # 3-body potential
   B3 = nbpolys(3, DD[2], DEGREES[2])
   c = rand(length(B3))
   V3 = NBPoly(B3, c, DD[2])
   V3sp = StNBPoly(V3)

   # 4-body potential
   B4 = nbpolys(4, DD[3], DEGREES[3])
   c = rand(length(B4))
   V4 = NBPoly(B4, c, DD[3])
   V4sp = StNBPoly(V4)

   # 5-body potential
   B5 = nbpolys(5, DD[4], DEGREES[4])
   c = rand(length(B5))
   V5 = NBPoly(B5, c, DD[4])
   V5sp = StNBPoly(V5)

   IP = NBodyIP([V2, V3, V4, V5])
   IPf = NBodyIP([V2sp, V3sp, V4sp, V5sp])

   return IP, IPf
end

_, IPf = random_ip(BondLengthDesc)
at = rattle!(bulk(:W, cubic=true) * 3, 0.01)
energy(IPf, at)
forces(IPf, at)
at = rattle!(bulk(:W, cubic=true) * 10, 0.01)
@info("Energy")
@time energy(IPf, at)
@time energy(IPf, at)
@info("Forces")
@time forces(IPf, at)
@time forces(IPf, at)

# test correctness first
IP_bl, _ = random_ip(BondLengthDesc)
IP_cl = NBodyIPs.Experimental.faster_blpot(IP_bl)

at = rattle!(bulk(:W, cubic=true) * 10, 0.01)
# set_pbc!(at, false)

V2_bl = IP_bl.components[1]
V2_cl = IP_cl.components[1]
@test energy(V2_bl, at) ≈ energy(V2_cl, at)

V3_bl = IP_bl.components[2]
V3_cl = IP_cl.components[2]
@test energy(V3_bl, at) ≈ energy(V3_cl, at)

V4_bl = IP_bl.components[3]
V4_cl = IP_cl.components[3]
@test energy(V4_bl, at) ≈ energy(V4_cl, at)

for DT in [BondLengthDesc, ClusterBLDesc]
   println("-----------------------------------------------------------")
   println(" Testing StaticPolynomials with D = $(DT)")
   println("-----------------------------------------------------------")

   DD = [ DT(TRANSFORM, CosCut(rcut-1, rcut))  for rcut in RCUT ]

   # 2-body potential
   B2 = nbpolys(2, DD[1], DEGREES[1])
   c = rand(length(B2))
   V2 = NBPoly(B2, c, DD[1])
   V2sp = StNBPoly(V2)

   # 3-body potential
   B3 = nbpolys(3, DD[2], DEGREES[2])
   c = rand(length(B3))
   V3 = NBPoly(B3, c, DD[2])
   V3sp = StNBPoly(V3)

   # 4-body potential
   B4 = nbpolys(4, DD[3], DEGREES[3])
   c = rand(length(B4))
   V4 = NBPoly(B4, c, DD[3])
   V4sp = StNBPoly(V4)

   IP = NBodyIP([V2, V3, V4])
   IPf = NBodyIP([V2sp, V3sp, V4sp])

   println("Neighbourlist")
   @time neighbourlist(at, cutoff(V2_bl))
   @time neighbourlist(at, cutoff(V2_bl))
   println("Dynamic Polynomials")
   @time energy(IP, at)
   @time energy(IP, at)
   @time forces(IP, at)
   @time forces(IP, at)
   println("Static Polynomials")
   @time energy(IPf, at)
   @time energy(IPf, at)
   @time forces(IPf, at)
   @time forces(IPf, at)
end

ext: jl

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
include("../DegeneracyHunter.jl")
using JuMP
#using Ipopt
using AmplNLWriter

#s = IpoptSolver(bound_relax_factor=0.0)
s = AmplNLSolver("ipopt")

m = Model(solver=s)

@variable(m, -10 <= a[1:5] <= 10, start=1.0)
@variable(m, -100 <= z <= 100, start=5.0)

# Fix this variable by setting upper bound = lower bound
@variable(m, 0 <= y <= 0)

@constraint(m, a[1] + a[2] - a[3] >= 10)
@NLconstraint(m, a[1]*a[4] - a[2]*y + a[2] >= 0)
@NLconstraint(m, a[3]*a[2] + a[1]*a[4] - a[5] == 0)
@NLconstraint(m, z == sum{a[i], i=1:5} - 0.001*sum{a[i]^2, i=1:5})
@constraint(m, a[3] + a[4] == 10*y)

@objective(m, Min, z)

solve(m)

DegeneracyHunter.printInfeasibleEquations(m, 1E-10)

println("Manually check constrain violations")

r = zeros(5)
a_ = zeros(5)

for i = 1:5
	a_[i] = getvalue(a[i])
end

z_ = getvalue(z)
y_ = getvalue(y)

r[1] = max(10 - a_[1] - a_[2] + a_[3], 0)
r[2] = max(-a_[2] + a_[2]*y_ - a_[1]*a_[4], 0)
r[3] = a_[3]*a_[2] + a_[1]*a_[4] - a_[5]
r[4] = z_ - sum(a_) + 0.001*sum(a_.^2)
r[5] = a_[3] + a_[4] - 10*y_

for i = 1:5
	println("r[",i,"] = ",r[i])
end


text:
module A

export AType

struct AType
	x::Int
end

struct SameNameType{T}
    x::T
end

end


text:
import Distributed

"""
## Summary

struct ExistingProcessManager <: Distributed.ClusterManager

## Fields

- wconfigs :: Vector{Distributed.WorkerConfig}

---

    ExistingProcessManager(wconfigs::Vector{Distributed.WorkerConfig})

Construct an `ExistingProcessManager` from a list of `WorkerConfig`s.

## Example

```jldoctest; setup = :(using ExistingProcessManagers; using Distributed)
julia> w1 = WorkerConfig();

julia> w1.host = "127.0.0.1";

julia> w1.port = 9601;

julia> w2 = WorkerConfig();

julia> w2.host = "127.0.0.1";

julia> w2.port = 9602;

julia> wconfigs = WorkerConfig[w1, w2];

julia> manager = ExistingProcessManager(wconfigs);
```

Now you can do:
```julia
julia> addprocs(manager)
```
"""
struct ExistingProcessManager <: Distributed.ClusterManager
    wconfigs::Vector{Distributed.WorkerConfig}
end


text:
#!/usr/bin/env julia

using PandocFilters

const defaults = Dict(
                      "Code" => Dict(
                                     "hide"=>true,
                                     "display"=>true
                                     ),
                      "CodeBlock" => Dict(
                                          "hide"=>true,
                                          "display"=>false
                                          )
                      )

function proccess_code(t,c,format,meta)
    if t=="Code" || t=="CodeBlock"
        #defaults for Code and CodeBlock
        options = defaults[t]
        (name,classes,keywords),code = c
        keywords = [k[1]=>k[2] for k in keywords]
        option = merge(options,keywords)
        response = []
        if !options["hide"]
          #append!(response, [Dict("t"=>t,"c"=>c)])
        end
        if options["display"]
            append!(response, [eval(parse(code))])
        end
        n = length(response)
        if n==1
            response = response[1]
        end
        response
    else
        nothing
    end
end

PandocFilters.filter(proccess_code)


text:
module api

import OpenCL.cl

depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("CLFFT not properly installed. Please run Pkg.build(\"CLFFT\") then restart Julia.")
end

macro clfft(func, arg_types)
    local args_in  = Symbol[Symbol("arg$i::$T")
                            for (i, T) in enumerate(arg_types.args)]
    local funcname = Symbol("clfft$func")
    @eval begin
        $(funcname)($(args_in...)) = ccall(($(string(funcname)), libCLFFT),
                                                 cl.CL_int, #clfftStatus
                                                 $arg_types,
                                                 $(args_in...))
    end
end

struct SetupData
    major::cl.CL_uint
    minor::cl.CL_uint
    patch::cl.CL_uint
    debug_flags::cl.CL_ulong
end

const PlanHandle = Csize_t
const Callback =   Ptr{Nothing}
const UserData =   Ptr{Nothing}

const Precision =  Cint
const Dim =        Cint
const Direction =  Cint
const Layout =     Cint
const ResultLocation =   Cint
const ResultTransposed = Cint

clfft_dim(x) = convert(Dim, x)
clfft_direction(x) = convert(Direction, x)

# ERROR CODES
const CLFFT_INVALID_GLOBAL_WORK_SIZE          = cl.CL_INVALID_GLOBAL_WORK_SIZE
const CLFFT_INVALID_MIP_LEVEL                 = cl.CL_INVALID_MIP_LEVEL
const CLFFT_INVALID_BUFFER_SIZE               = cl.CL_INVALID_BUFFER_SIZE
const CLFFT_INVALID_GL_OBJECT                 = cl.CL_INVALID_GL_OBJECT
const CLFFT_INVALID_OPERATION                 = cl.CL_INVALID_OPERATION
const CLFFT_INVALID_EVENT                     = cl.CL_INVALID_EVENT
const CLFFT_INVALID_EVENT_WAIT_LIST           = cl.CL_INVALID_EVENT_WAIT_LIST,
const CLFFT_INVALID_GLOBAL_OFFSET             = cl.CL_INVALID_GLOBAL_OFFSET,
const CLFFT_INVALID_WORK_ITEM_SIZE            = cl.CL_INVALID_WORK_ITEM_SIZE
const CLFFT_INVALID_WORK_GROUP_SIZE           = cl.CL_INVALID_WORK_GROUP_SIZE
const CLFFT_INVALID_WORK_DIMENSION            = cl.CL_INVALID_WORK_DIMENSION
const CLFFT_INVALID_KERNEL_ARGS               = cl.CL_INVALID_KERNEL_ARGS
const CLFFT_INVALID_ARG_SIZE                  = cl.CL_INVALID_ARG_SIZE
const CLFFT_INVALID_ARG_VALUE                 = cl.CL_INVALID_ARG_VALUE
const CLFFT_INVALID_ARG_INDEX                 = cl.CL_INVALID_ARG_INDEX
const CLFFT_INVALID_KERNEL                    = cl.CL_INVALID_KERNEL
const CLFFT_INVALID_KERNEL_DEFINITION         = cl.CL_INVALID_KERNEL_DEFINITION
const CLFFT_INVALID_KERNEL_NAME               = cl.CL_INVALID_KERNEL_NAME
const CLFFT_INVALID_PROGRAM_EXECUTABLE        = cl.CL_INVALID_PROGRAM_EXECUTABLE
const CLFFT_INVALID_PROGRAM                   = cl.CL_INVALID_PROGRAM
const CLFFT_INVALID_BUILD_OPTIONS             = cl.CL_INVALID_BUILD_OPTIONS
const CLFFT_INVALID_BINARY                    = cl.CL_INVALID_BINARY
const CLFFT_INVALID_SAMPLER                   = cl.CL_INVALID_SAMPLER
const CLFFT_INVALID_IMAGE_SIZE                = cl.CL_INVALID_IMAGE_SIZE
const CLFFT_INVALID_IMAGE_FORMAT_DESCRIPTOR   = cl.CL_INVALID_IMAGE_FORMAT_DESCRIPTOR
const CLFFT_INVALID_MEM_OBJECT                = cl.CL_INVALID_MEM_OBJECT
const CLFFT_INVALID_HOST_PTR                  = cl.CL_INVALID_HOST_PTR
const CLFFT_INVALID_COMMAND_QUEUE             = cl.CL_INVALID_COMMAND_QUEUE
const CLFFT_INVALID_QUEUE_PROPERTIES          = cl.CL_INVALID_QUEUE_PROPERTIES
const CLFFT_INVALID_CONTEXT                   = cl.CL_INVALID_CONTEXT
const CLFFT_INVALID_DEVICE                    = cl.CL_INVALID_DEVICE
const CLFFT_INVALID_PLATFORM                  = cl.CL_INVALID_PLATFORM
const CLFFT_INVALID_DEVICE_TYPE               = cl.CL_INVALID_DEVICE_TYPE
const CLFFT_INVALID_VALUE                     = cl.CL_INVALID_VALUE
const CLFFT_MAP_FAILURE                       = cl.CL_MAP_FAILURE
const CLFFT_BUILD_PROGRAM_FAILURE             = cl.CL_BUILD_PROGRAM_FAILURE
const CLFFT_IMAGE_FORMAT_NOT_SUPPORTED        = cl.CL_IMAGE_FORMAT_NOT_SUPPORTED
const CLFFT_IMAGE_FORMAT_MISMATCH             = cl.CL_IMAGE_FORMAT_MISMATCH
const CLFFT_MEM_COPY_OVERLAP                  = cl.CL_MEM_COPY_OVERLAP
const CLFFT_PROFILING_INFO_NOT_AVAILABLE      = cl.CL_PROFILING_INFO_NOT_AVAILABLE
const CLFFT_OUT_OF_HOST_MEMORY                = cl.CL_OUT_OF_HOST_MEMORY
const CLFFT_OUT_OF_RESOURCES                  = cl.CL_OUT_OF_RESOURCES
const CLFFT_MEM_OBJECT_ALLOCATION_FAILURE     = cl.CL_MEM_OBJECT_ALLOCATION_FAILURE
const CLFFT_COMPILER_NOT_AVAILABLE            = cl.CL_COMPILER_NOT_AVAILABLE
const CLFFT_DEVICE_NOT_AVAILABLE              = cl.CL_DEVICE_NOT_AVAILABLE
const CLFFT_DEVICE_NOT_FOUND                  = cl.CL_DEVICE_NOT_FOUND
const CLFFT_SUCCESS                           = cl.CL_SUCCESS

const CLFFT_BUGCHECK                  = Int32(4*1024)
const CLFFT_NOTIMPLEMENTED            = Int32(4*1024+1)  # Functionality is not implemented yet.
const CLFFT_TRANSPOSED_NOTIMPLEMENTED = Int32(4*1024+2)  # Transposed functionality is not implemented for this transformation.
const CLFFT_FILE_NOT_FOUND            = Int32(4*1024+3)  # Tried to open an existing file on the host system, but failed.
const CLFFT_FILE_CREATE_FAILURE       = Int32(4*1024+4)  # Tried to create a file on the host system, but failed.
const CLFFT_VERSION_MISMATCH          = Int32(4*1024+5)  # Version conflict between client and library.
const CLFFT_INVALID_PLAN              = Int32(4*1024+6)  # Requested plan could not be found.
const CLFFT_DEVICE_NO_DOUBLE          = Int32(4*1024+7)  # Double precision not supported on this device.

const CLFFT_1D = Int32(1) # 1 Dimensional FFT transform (default).
const CLFFT_2D = Int32(2) # 2 Dimensional FFT transform.
const CLFFT_3D = Int32(3) # 3 Dimensional FFT transform.

const CLFFT_COMPLEX_INTERLEAVED   = Int32(1) # An array of complex numbers, with real and imaginary components together (default).
const CLFFT_COMPLEX_PLANAR        = Int32(2) # Arrays of real componets and arrays of imaginary components that have been seperated out.
const CLFFT_HERMITIAN_INTERLEAVED = Int32(3) # Compressed form of complex numbers; complex-conjugates not stored, real and imaginary components in same array.
const CLFFT_HERMITIAN_PLANAR      = Int32(4) # Compressed form of complex numbers; complex-conjugates not stored, real and imaginary components in separate arrays.
const CLFFT_REAL                  = Int32(5) # An array of real numbers, with no corresponding imaginary components.

const CLFFT_SINGLE       = Int32(1) # An array of complex numbers, with real and imaginary components as floats (default).
const CLFFT_DOUBLE       = Int32(2) # An array of complex numbers, with real and imaginary components as doubles.
const CLFFT_SINGLE_FAST  = Int32(3) # Faster implementation preferred.
const CLFFT_DOUBLE_FAST  = Int32(4) # Faster implementation preferred.

const CLFFT_FORWARD      = Int32(-1) # FFT transform from the time to the frequency domain.
const CLFFT_BACKWARD     = Int32(1)  # FFT transform from the frequency to the time domain
const CLFFT_MINUS        = Int32(-1) # Alias for the forward transform.
const CLFFT_PLUS         = Int32(1)  # Alias for the backward transform.

const CLFFT_INPLACE      = Int32(1) # The input and output buffers are the same (default).
const CLFFT_OUTOFPLACE   = Int32(2) # Seperate input and output buffers.

const CLFFT_NOTRANSPOSE   = Int32(1) # The results are returned in the original preserved order (default)
const CLFFT_TRANSPOSED    = Int32(2) # The result is transposed where transpose kernel is supported (possibly faster)

# @brief Initialize an clfftSetupData struct for the client
# @details clfftSetupData is passed to clfftSetup to control behavior of the FFT runtime
# @param[out] setupData Data structure is cleared,
#             initialized with version information and default values
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(InitSetupData, (Ptr{SetupData},))

# @brief Initialize internal FFT resources.
# @details AMD's FFT implementation caches kernels, programs and buffers for its internal use.
# @param[in] setupData Data structure that can be passed into the setup routine to control FFT generation behavior
#            and debug functionality
# @return    Enum describing error condition; superset of OpenCL error codes
@clfft(Setup, (Ptr{SetupData},))

# @brief Release all internal resources.
# @details Call when client is done with this FFT library, allowing the library to destroy all resources it has cached
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(Teardown, ())

# @brief Query the FFT library for version information
# @details Return the major, minor and patch version numbers associated with this FFT library
# @param[out] major Major functionality change
# @param[out] minor Minor functionality change
# @param[out] patch Bug fixes, documentation changes, no new features introduced
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetVersion, (Ptr{cl.CL_uint}, Ptr{cl.CL_uint}, Ptr{cl.CL_uint}))

# @brief Create a plan object initialized entirely with default values.
# @details A plan is a repository of state for calculating FFT's.  Allows the runtime to pre-calculate kernels, programs
# and buffers and associate them with buffers of specified dimensions.
# @param[out] plHandle Handle to the newly created plan
# @param[in] context Client is responsible for providing an OpenCL context for the plan
# @param[in] dim The dimensionality of the FFT transform; describes how many elements are in the array
# @param[in] clLengths An array of lengths, of size 'dim'.  Each value describes the length of additional dimensions
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(CreateDefaultPlan, (Ptr{PlanHandle}, cl.CL_context, Dim, Ptr{Csize_t}))

# @brief Create a copy of an existing plan.
# @details This API allows a client to create a new plan based upon an existing plan.  This is a convenience function
# provided for quickly creating plans that are similar, but may differ slightly.
# @param[out] out_plHandle Handle to the newly created plan that is based on in_plHandle
# @param[in] new_context Client is responsible for providing a new context for the new plan
# @param[in] in_plHandle Handle to a plan to be copied, previously created
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(CopyPlan, (Ptr{PlanHandle}, cl.CL_context, Dim, Ptr{Csize_t}))

# @brief Prepare the plan for execution.
# @details After all plan parameters are set, the client has the option of 'baking' the plan, which tells the runtime that
# no more changes to the plan's parameters are expected, and the OpenCL kernels should be compiled.
# This optional function allows the client application to perform this function when the
# application is being initialized instead of on the first execution.
# At this point, the clfft runtime will apply all implimented optimizations,
# possibly including running kernel experiments on the devices in the plan context.
# Users should assume that this function will take a long time to execute.
# If a plan is not baked before being executed, users should assume that the first call to clfftEnqueueTransform will take a long time to execute.
# If any significant parameter of a plan is changed after the plan is baked (by a subsequent call to one of
# the clfftSetPlan____ functions), that will not be considered an error.
# Instead, the plan will revert back to the unbaked state, discarding the benefits of the baking operation.

# @param[in] plHandle Handle to a plan previously created
# @param[in] numQueues Number of command queues in commQueueFFT; 0 is a valid value, in which case client does not want
#            the runtime to run load experiments and only pre-calculate state information
# @param[in] commQueueFFT An array of cl_command_queues created by the client; the command queues must be a proper subset of
#            the devices included in the plan context
# @param[in] pfn_notify A function pointer to a notification routine. The notification routine is a callback function that
#            an application can register and which will be called when the program executable has been built (successfully or unsuccessfully).
#            Currently, this parameter MUST be NULL or nullptr.
# @param[in] user_data Passed as an argument when pfn_notify is called.
#            Currently, this parameter MUST be NULL or nullptr.
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(BakePlan, (PlanHandle, cl.CL_uint, Ptr{cl.CL_command_queue}, Callback, UserData))

# @brief Release the resources of a plan.
# @details A plan may include kernels, programs and buffers associated with it that consume memory. When a plan  is not needed anymore, the client should release the plan.
# @param[in,out] plHandle Handle to a plan previously created
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(DestroyPlan, (Ptr{PlanHandle},))

# @brief Retrieve the OpenCL context of a previously created plan.
# @details User should pass a reference to an cl_context variable, which will be changed to point to a
# context set in the specified plan.
# @param[in] plHandle Handle to a plan previously created
# @param[out] context Reference to user allocated cl_context, which will point to context set in plan
 # @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanContext, (PlanHandle, Ptr{cl.CL_context}))

# @brief Retrieve the floating point precision of the FFT data
# @details User should pass a reference to an clfftPrecision variable, which will be set to the
# precision of the FFT complex data in the plan.
# @param[in] plHandle Handle to a plan previously created
# @param[out] precision Reference to user clfftPrecision enum
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanPrecision, (PlanHandle, Ptr{Precision}))

# @brief Set the floating point precision of the FFT data
# @details Set the plan property which will be the precision of the FFT complex data in the plan.
# @param[in] plHandle Handle to a plan previously created
# @param[in] precision Reference to user clfftPrecision enum
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(SetPlanPrecision, (PlanHandle, Precision))

# @brief Retrieve the scaling factor that should be applied to the FFT data
# @details User should pass a reference to an cl_float variable, which will be set to the
# floating point scaling factor that will be multiplied across the FFT data.
# @param[in] plHandle Handle to a plan previously created
# @param[in] dir Which direction does the scaling factor apply to
# @param[out] scale Reference to user cl_float variable
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanScale, (PlanHandle, Direction, Ptr{cl.CL_float}))

# @brief Set the scaling factor that should be applied to the FFT data
# @details Set the plan property which will be the floating point scaling factor that will be
# multiplied across the FFT data.
# @param[in] plHandle Handle to a plan previously created
# @param[in] dir Which direction does the scaling factor apply to
# @param[in] scale Reference to user cl_float variable
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(SetPlanScale, (PlanHandle, Direction, cl.CL_float))

# @brief Retrieve the number of discrete arrays that this plan can handle concurrently
# @details User should pass a reference to an cl_uint variable, which will be set to the
# number of discrete arrays (1D or 2D) that will be batched together for this plan
# @param[in] plHandle Handle to a plan previously created
# @param[out] batchSize How many discrete number of FFT's are to be performed
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanBatchSize, (PlanHandle, Ptr{Csize_t}))

# @brief Set the number of discrete arrays that this plan can handle concurrently
# @details Set the plan property which will be set to the number of discrete arrays (1D or 2D)
# that will be batched together for this plan
# @param[in] plHandle Handle to a plan previously created
# @param[in] batchSize How many discrete number of FFT's are to be performed
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(SetPlanBatchSize, (PlanHandle, Csize_t))

# @brief Retrieve the dimensionality of FFT's to be transformed in the plan
# @details Queries a plan object and retrieves the dimensionality that the plan is set for.  A size is returned to
# help the client allocate the proper storage to hold the dimensions in a further call to clfftGetPlanLength
# @param[in] plHandle Handle to a plan previously created
# @param[out] dim The dimensionality of the FFT's to be transformed
# @param[out] size Value used to allocate an array to hold the FFT dimensions.
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanDim, (PlanHandle, Ptr{Dim}, Ptr{cl.CL_uint}))

# @brief Set the dimensionality of FFT's to be transformed by the plan
# @details Set the dimensionality of FFT's to be transformed by the plan
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimensionality of the FFT's to be transformed
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(SetPlanDim, (PlanHandle, Dim))

# @brief Retrieve the length of each dimension of the FFT
# @details User should pass a reference to a size_t array, which will be set to the
# length of each discrete dimension of the FFT
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the length parameters; describes how many elements are in the array
# @param[out] clLengths An array of lengths, of size 'dim'.  Each array value describes the length of each dimension
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(GetPlanLength, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Set the length of each dimension of the FFT
# @details Set the plan property which will be the length of each discrete dimension of the FFT
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the length parameters; describes how many elements are in the array
# @param[in] clLengths An array of lengths, of size 'dim'.  Each value describes the length of additional dimensions
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(SetPlanLength, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Retrieve the distance between consecutive elements for input buffers in a dimension.
# @details Depending on how the dimension is set in the plan (for 2D or 3D FFT's), strideY or strideZ can be safely ignored
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the stride parameters; describes how many elements are in the array
# @param[out] clStrides An array of strides, of size 'dim'.
@clfft(GetPlanInStride, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Set the distance between consecutive elements for input buffers in a dimension.
# @details Set the plan properties which will be the distance between elements in a given dimension (units are in terms of clfftPrecision)
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the stride parameters; describes how many elements are in the array
# @param[in] clStrides An array of strides, of size 'dim'. Usually strideX=1 so that successive elements in the first dimension are stored contiguously.
# Typically strideY=LenX, strideZ=LenX*LenY such that successive elements in the second and third dimensions are stored in packed format.
@clfft(SetPlanInStride, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Retrieve the distance between consecutive elements for output buffers in a dimension.
# @details Depending on how the dimension is set in the plan (for 2D or 3D FFT's), strideY or strideZ can be safely ignored
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the stride parameters; describes how many elements are in the array
# @param[out] clStrides An array of strides, of size 'dim'.
@clfft(GetPlanOutStride, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Set the distance between consecutive elements for output buffers in a dimension.
# @details Set the plan properties which will be the distance between elements in a given dimension
# (units are in terms of clfftPrecision)
# @param[in] plHandle Handle to a plan previously created
# @param[in] dim The dimension of the stride parameters; describes how many elements are in the array
# @param[in] clStrides An array of strides, of size 'dim'.  Usually strideX=1 so that successive elements in the first dimension are stored contiguously.
# Typically strideY=LenX, strideZ=LenX*LenY such that successive elements in the second and third dimensions are stored in packed format.
@clfft(SetPlanOutStride, (PlanHandle, Dim, Ptr{Csize_t}))

# @brief Retrieve the distance between Array objects
# @details Pitch is the distance between each discrete array object in an FFT array. This is only used
# for 'array' dimensions in clfftDim; see clfftSetPlanDimension (units are in terms of clfftPrecision)
# @param[in]  plHandle Handle to a plan previously created
# @param[out] iDist The distance between the beginning elements of the discrete array objects in memory on input.
#             For contiguous arrays in memory, iDist=(strideX*strideY*strideZ)
# @param[out] oDist The distance between the beginning elements of the discrete array objects in memory on output.
#             For contiguous arrays in memory, oDist=(strideX*strideY*strideZ)
@clfft(GetPlanDistance, (PlanHandle, Ptr{Csize_t}, Ptr{Csize_t}))

# @brief Set the distance between Array objects
# @details Pitch is the distance between each discrete array object in an FFT array. This is only used
# for 'array' dimensions in clfftDim; see clfftSetPlanDimension (units are in terms of clfftPrecision)
# @param[in] plHandle Handle to a plan previously created
# @param[out] iDist The distance between the beginning elements of the discrete array objects in memory on input.
#             For contiguous arrays in memory, iDist=(strideX*strideY*strideZ)
# @param[out] oDist The distance between the beginning elements of the discrete array objects in memory on output.
#             For contiguous arrays in memory, oDist=(strideX*strideY*strideZ)
@clfft(SetPlanDistance, (PlanHandle, Csize_t, Csize_t))

# @brief Retrieve the expected layout of the input and output buffers
# @details Output buffers can be filled with either hermitian or complex numbers.  Complex numbers can be stored
# in various layouts; this informs the FFT engine what layout to produce on output
# @param[in] plHandle Handle to a plan previously created
# @param[out] iLayout Indicates how the input buffers are laid out in memory
# @param[out] oLayout Indicates how the output buffers are laid out in memory
@clfft(GetLayout, (PlanHandle, Ptr{Layout}, Ptr{Layout}))

# @brief Set the expected layout of the input and output buffers
# @details Output buffers can be filled with either hermitian or complex numbers.  Complex numbers can be stored
# in various layouts; this informs the FFT engine what layout to produce on output
# @param[in] plHandle Handle to a plan previously created
# @param[in] iLayout Indicates how the input buffers are laid out in memory
# @param[in] oLayout Indicates how the output buffers are laid out in memory
@clfft(SetLayout, (PlanHandle, Layout, Layout))

# @brief Retrieve whether the input buffers are going to be overwritten with results
# @details If the setting is to do an in-place transform, the input buffers are overwritten with the results of the
# transform.  If the setting is for out-of-place transforms, the engine knows to look for separate output buffers on the Enqueue call.
# @param[in] plHandle Handle to a plan previously created
# @param[out] placeness Tells the FFT engine to clobber the input buffers or to expect output buffers for results
@clfft(GetResultLocation, (PlanHandle, Ptr{ResultLocation}))

# @brief Set whether the input buffers are going to be overwritten with results
# @details If the setting is to do an in-place transform, the input buffers are overwritten with the results of the
# transform.  If the setting is for out-of-place transforms, the engine knows to look for separate output buffers on the Enqueue call.
# @param[in] plHandle Handle to a plan previously created
# @param[in] placeness Tells the FFT engine to clobber the input buffers or to expect output buffers for results
@clfft(SetResultLocation, (PlanHandle, ResultLocation))

# @brief Retrieve the final transpose setting of a muti-dimensional FFT
# @details A multi-dimensional FFT typically transposes the data several times during calculation.  If the client
# does not care about the final transpose to put data back in proper dimension, the final transpose can be skipped
# for possible speed improvements
# @param[in] plHandle Handle to a plan previously created
# @param[out] transposed Parameter specifies whether the final transpose can be skipped
@clfft(GetPlanTransposeResult, (PlanHandle, Ptr{ResultTransposed}))

# @brief Set the final transpose setting of a muti-dimensional FFT
# @details A multi-dimensional FFT typically transposes the data several times during calculation.  If the client
# does not care about the final transpose to put data back in proper dimension, the final transpose can be skipped
# for possible speed improvements
# @param[in] plHandle Handle to a plan previously created
# @param[in] transposed Parameter specifies whether the final transpose can be skipped
@clfft(SetPlanTransposeResult, (PlanHandle, ResultTransposed))

# @brief Get buffer size (in bytes), which may be needed internally for an intermediate buffer
# @details Very large FFT transforms may need multiple passes, and the operation would need a temporary buffer to hold
# intermediate results. This function is only valid after the plan is baked, otherwise an invalid operation error
# is returned. If buffersize returns as 0, the runtime needs no temporary buffer.
# @param[in]  plHandle Handle to a plan previously created
# @param[out] buffersize Size in bytes for intermediate buffer
@clfft(GetTmpBufSize, (PlanHandle, Ptr{Csize_t}))

# @brief Enqueue an FFT transform operation, and return immediately (non-blocking)
# @details This transform API is the function that actually computes the FFT transfrom. It is non-blocking as it
#          only enqueues the OpenCL kernels for execution. The synchronization step has to be managed by the user.
# @param[in]  plHandle Handle to a plan previously created
# @param[in]  dir Forwards or backwards transform
# @param[in]  numQueuesAndEvents Number of command queues in commQueues; number of expected events to be returned in outEvents
# @param[in]  commQueues An array of cl_command_queues created by the client; the command queues must be a proper subset of
#             the devices included in the plan context
# @param[in]  numWaitEvents Specify the number of elements in the eventWaitList array
# @param[in]  waitEvents Events that this transform should wait to complete before executing on the device
# @param[out] outEvents The runtime fills this array with events corresponding 1 to 1 with the input command queues passed
#             in commQueues.  This parameter can be NULL or nullptr, in which case client is not interested in receiving notifications
#             when transforms are finished, otherwise if not NULL the client is responsible for allocating this array, with at least
#             as many elements as specified in numQueuesAndEvents.
# @param[in]  inputBuffers An array of cl_mem objects that contain data for processing by the FFT runtime.  If the transform
#             is in place, the FFT results will overwrite the input buffers
# @param[out] outputBuffers An array of cl_mem objects that will store the results of out of place transforms.  If the transform
#             is in place, this parameter may be NULL or nullptr.  It is completely ignored
# @param[in]  tmpBuffer A cl_mem object that is reserved as a temporary buffer for FFT processing. If clTmpBuffers is NULL or nullptr,
#             and the runtime needs temporary storage, an internal temporary buffer will be created on the fly managed by the runtime.
# @return Enum describing error condition; superset of OpenCL error codes
@clfft(EnqueueTransform, (PlanHandle,
                          Direction,
                          cl.CL_uint,
                          Ptr{cl.CL_command_queue},
                          cl.CL_uint,
                          Ptr{cl.CL_event},
                          Ptr{cl.CL_event},
                          Ptr{cl.CL_mem},
                          Ptr{cl.CL_mem},
                          cl.CL_mem))
end # end module


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

plot(x=rand(100), y=rand(100),
    shape=rand([Shape.circle,Shape.square,Shape.utriangle], 100),
    Geom.point)


text:
# SPOT example from file doc/examples/sdsos_example.m
# of https://github.com/anirudhamajumdar/spotless/tree/spotless_isos
# See Section 4.1 of [AM17].
#
# [AM17] Section 4.1 of
# A. A. Ahmadi, and A. Majumdar
# DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization
# 2017

using Test
using SumOfSquares
using DynamicPolynomials

function quartic_comparison_test(
    optimizer, config::MOIT.TestConfig,
    cone::SumOfSquares.PolyJuMP.PolynomialSet, expected_objective_value)
    atol = config.atol
    rtol = config.rtol

    @polyvar x[1:3]
    vx = monomials(x, 4) # Degree 4 homogeneous
    # Coefficient of polynomial
    cp = 15:-1:1
    p = polynomial(cp, vx)

    model = _model(optimizer)
    @variable(model, γ)
    @constraint(model, p - γ * sum(x .* x)^2 in cone)
    @objective(model, Max, γ)
    JuMP.optimize!(model)
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    @test JuMP.objective_value(model) ≈ expected_objective_value atol=atol rtol=rtol
end

sos_quartic_comparison_test(optimizer, config)   = quartic_comparison_test(optimizer, config, SOSCone(), -0.184667)
sd_tests["sos_quartic_comparison"]      = sos_quartic_comparison_test
sdsos_quartic_comparison_test(optimizer, config) = quartic_comparison_test(optimizer, config, SDSOSCone(), -3.172412)
soc_tests["sdsos_quartic_comparison"]   = sdsos_quartic_comparison_test
dsos_quartic_comparison_test(optimizer, config)  = quartic_comparison_test(optimizer, config, DSOSCone(), -11/3)
linear_tests["dsos_quartic_comparison"] = dsos_quartic_comparison_test


text:
@testset "LOAD" begin

  mktempdir() do path
    @testset "array" begin
      file = joinpath(path, "arr.bin")
      x0 = libtorch.randn(())
      x1 = libtorch.randn((5,))
      x2 = libtorch.randn((3,3))
      x3 = libtorch.randn((4,3,2))
      x4 = libtorch.randn((5,4,3,2))
      libtorch.save(Dict("x0"=>x0, "x1"=>x1,
               "x2"=>x2, "x3"=>x3, "x4"=>x4), file;
               _use_new_zipfile_serialization=false)
      load_arr = Torch.THload(file)
      @test load_arr["x0"] ≈ x0.numpy()
      @test load_arr["x1"] ≈ x1.numpy()
      @test load_arr["x2"] ≈ x2.numpy()
      @test load_arr["x3"] ≈ x3.numpy()
      @test load_arr["x4"] ≈ x4.numpy()
    end

    x = libtorch.randn(10,10)
    @testset "slice" begin
      file = joinpath(path, "slice.bin")
      x_slice = random_slice(x)
      libtorch.save([x, x_slice], file;
                    _use_new_zipfile_serialization=false)
      test_load_slice = Torch.THload(file)
      load_x, load_slice = test_load_slice
      @test load_x ≈ x.numpy()
      @test load_slice ≈ x_slice.numpy()
      @test pointer(load_x.parent) == pointer(load_slice.parent)
    end

    @testset "stride" begin
      file = joinpath(path, "stride.bin")
      x_stride = random_stride(x)
      libtorch.save([x, x_stride], file;
                    _use_new_zipfile_serialization=false)
      test_load_stride = Torch.THload(file)
      load_x, load_stride = test_load_stride
      @test load_x ≈ x.numpy()
      @test load_stride ≈ x_stride.numpy()
      @test pointer(load_x.parent) == pointer(load_stride.parent)
    end

    @testset "reshape" begin
      file = joinpath(path, "reshape.bin")
      x_reshape = x.reshape(2,5,5,2)
      libtorch.save([x, x_reshape], file;
                    _use_new_zipfile_serialization=false)
      test_load_reshape = Torch.THload(file)
      load_x, load_reshape = test_load_reshape
      @test load_x ≈ x.numpy()
      @test load_reshape ≈ x_reshape.numpy()
      @test pointer(load_x.parent) == pointer(load_reshape.parent)
    end

    @testset "offset" begin
      file = joinpath(path, "offset.bin")
      x_offset = pyslice(x.reshape(-1), 5).reshape(5, -1)
      libtorch.save([x, x_offset], file;
                    _use_new_zipfile_serialization=false)
      test_load_offset = Torch.THload(file)
      load_x, load_offset = test_load_offset
      @test load_x ≈ x.numpy()
      @test load_offset ≈ x_offset.numpy()
      @test pointer(load_x.parent) == pointer(load_offset.parent)
    end

    @testset "mutation" begin
      file = joinpath(path, "mutate.bin")
      x1 = x.clone()
      x2 = pyslice(x1, 5)
      libtorch.save([x1, x2], file;
                    _use_new_zipfile_serialization=false)
      test_load_mutate = Torch.THload(file)
      load_x1, load_x2 = test_load_mutate
      set!(x2, 0, 0)
      load_x2[1, :] .= 0
      @test all(iszero, load_x1[6, :])
      @test load_x1 ≈ x1.numpy()
      @test load_x2 ≈ x2.numpy()
    end

    @testset "transpose" begin
      file = joinpath(path, "transpose.bin")
      x_transpose = x.transpose(1, 0)
      libtorch.save([x, x_transpose], file;
                    _use_new_zipfile_serialization=false)
      test_load_transpose = Torch.THload(file)
      load_x, load_transpose = test_load_transpose
      @test load_x ≈ x.numpy()
      @test load_transpose ≈ x_transpose.numpy()
      @test pointer(load_x.parent) == pointer(load_transpose)
    end

    @testset "all" begin
      file = joinpath(path, "all.bin")
      x_offset = pyslice(x.reshape(-1), 10).reshape(6, 15)
      x_stride = random_stride(x_offset)
      x_slice = random_slice(x_stride)
      x_transpose = x_slice.transpose(0,1)
      set!(x_transpose, (0,0), 0)
      libtorch.save([x, x_offset, x_stride, x_slice, x_transpose], file;
                    _use_new_zipfile_serialization=false)
      test_load_all = Torch.THload(file)
      load_x, load_offset, load_stride, load_slice, load_transpose = test_load_all
      @test load_x ≈ x.numpy()
      @test load_offset ≈ x_offset.numpy()
      @test load_stride ≈ x_stride.numpy()
      @test load_slice ≈ x_slice.numpy()
      @test load_transpose ≈ x_transpose.numpy()
      @test pointer(load_x.parent) == pointer(load_offset.parent)
      @test pointer(load_x.parent) == pointer(load_stride.parent)
      @test pointer(load_x.parent) == pointer(load_slice.parent)
      @test pointer(load_x.parent) == pointer(load_transpose.parent)
    end

  end

end


text:
module VarinnerTests

using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField, MockBasis
using InplaceArrays.TensorValues
using FillArrays

p1 = Point(2,2)
p2 = Point(4,2)
p3 = Point(1,3)
p4 = Point(5,2)
x = [p1,p2,p3,p4]
np = length(x)

va = 3.0
vb = 3.3
d = 2
fa = MockField{d}(va)
fb = MockField{d}(vb)
fax = evaluate(fa,x)
fbx = evaluate(fb,x)

wa = 4.5
wb = 4.5
ndofa = 8
ndofb = 5
ba = MockBasis{d}(wa,ndofa)
bb = MockBasis{d}(wb,ndofb)
bax = evaluate(ba,x)
bbx = evaluate(bb,x)

l = 10
ax = fill(x,l)
afa = Fill(fa,l)
afb = Fill(fb,l)
aba = Fill(ba,l)
abb = Fill(bb,l)

g = varinner(fa,fb)
gx = [inner(ai,bj) for (ai,bj) in zip(fax,fbx)]
test_field(g,x,gx)

ag = varinner(afa,afb)
agx = fill(gx,l)
test_array_of_fields(ag,ax,agx)

g = varinner(ba,fb)
gx = zeros(size(bax))
for p in 1:np
  for i in 1:ndofa
    gx[p,i] = inner(bax[p,i],fbx[p])
  end
end
test_field(g,x,gx)

ag = varinner(aba,afb)
agx = fill(gx,l)
test_array_of_fields(ag,ax,agx)

g = varinner(ba,bb)
gx = zeros(np,ndofa,ndofb)
for p in 1:np
  for i in 1:ndofa
    for j in 1:ndofb
    gx[p,i,j] = inner(bax[p,i],bbx[p,j])
    end
  end
end
test_field(g,x,gx)

ag = varinner(aba,abb)
agx = fill(gx,l)
test_array_of_fields(ag,ax,agx)

end # module


text:
push!(LOAD_PATH,"../src/")

using Documenter, LinearAlgebraicRepresentation

makedocs(
	format = :html,
	sitename = "LinearAlgebraicRepresentation.jl",
	assets = ["assets/lar.css", "assets/logo.png"],
	pages = [
		"Home" => "index.md",
		"L.A.R. Intro" => "lar.md",
		"Interface" => "interface.md",
		"Arrangement" => "arrangement.md",
		"Parametric primitives" => [
			"Mapper" => "mapper.md",
			"Assemblies" => "struct.md"
		],
		"Grid generation" => [
			"Cuboidal grids" => "largrid.md",
			"Simplicial grids" => "simplexn.md"
		],
		"Domain integration" => "integr.md",
	]
)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 53c7d9b0-6b89-11eb-398b-675411f493dc
begin
	using Plots
	# Use Plotly plotting backend
	plotlyjs()
end

# ╔═╡ 15c25e30-6bbe-11eb-2e55-cfe5abc44d2b
md"""
Using Gauss' Theorem

$\oint_S \overline{E} \cdot d\overline{s} = \frac{Q_{in}}{\epsilon}$

"""

# ╔═╡ 951df740-6b8a-11eb-2caa-b5a1291edc31
md"""
## E-field in spherical geometry
This claculation describes the strength of an electric field with chrage density $\rho_v$. It is of a spehrical form, with vaccum on the inside (radius less than a) and the outside (radius bigger than b).
"""

# ╔═╡ f87b9c1e-6d6f-11eb-0805-0923965aba3b
ϵ = 8.8541878128e-12

# ╔═╡ 5a5fb850-6b8f-11eb-2d2c-8dcedf10efd9
ρ = 5*10^(-3)

# ╔═╡ 1ed77c60-6bc5-11eb-1a68-e9e8d51f7d87
@bind a html"<input value='1' type='range' max='1' step='' min='0'>"

# ╔═╡ 8fc3f460-6bc7-11eb-26b8-61d35bfa4fef
@bind b html"<input value='2' type='range' step='0.1' max='3' min='2' label='$(b)'>"

# ╔═╡ f0717252-6bc8-11eb-3bfd-d5d0d0fd299c
(a,b)

# ╔═╡ 660d2082-6b89-11eb-100d-8b577e9ae037
function E(r)
	if(r<=a)
		return 0
	elseif(a<r<b)
		return ρ*(r^3-a^3)/(ϵ*(3*r^2))
	else
		return ρ*(b^3-a^3)/(ϵ*(3*r^2))
	end
end

# ╔═╡ 6a2b5b50-6b89-11eb-38e5-33a649d074ae
r = 0:0.001:10

# ╔═╡ 7418b540-6b89-11eb-1f53-89165d55d825
Y = E.(r)

# ╔═╡ 78ff2210-6b89-11eb-29df-c5714c003c7d
plot(r, Y, label="|E|", xlabel="Radius in m (r)", ylabel="Field strength |E|")

# ╔═╡ 6d19078e-6b98-11eb-265b-c913b44ff56a
md"""
## Voltage everywhere around a spherical shell
The shell has a voltage of  $V_o$, determine the voltage as a function of the radius. It is located at r=a.

1. Decide E-field
2. Integrate for the voltage
3. Solve for boundary values

By using symmetry when creating a Gauss surface for r>a

$E_r(r)\cdot 4\pi r^2=\frac{Q_{in}}{\epsilon}$

Which gives us the electric field

$E_r(r)=\frac{Q_{in}}{ 4\pi\epsilon r^2}$

Next we use

$\overline{E} = - \nabla V$
Because of symmetry, it can be described as

$E_r(r)=-\frac{\partial V}{\partial r}$

Solve for V:

$V(r) = - \int E_r(r)dr = -\int \frac{Q_{in}}{ 4\pi\epsilon r^2}dr = \frac{Q_{in}}{ 4\pi\epsilon r} + C$

Use boundary values

$V(\infty) = C = 0$

$V(a) = V_0 = \frac{Q_{in}}{ 4\pi\epsilon a} \Rightarrow Q_{in} = 4V_0\pi\epsilon a$

For r>a V is:

$V(r) = \frac{V_0 a}{r}$

For r<a V is:

$V(r) = V_0$

"""

# ╔═╡ eaee7830-6b98-11eb-0a38-3b3060407cc1
begin
	#Initilize variables
	V_0 = 10 #Voltage on spehrical shell
	a_r = 2 #Radius where the shell is located
end

# ╔═╡ 0de35890-6bc3-11eb-2e29-0f5ebe681620
function V(r)
	if r<=a_r
		return V_0
	end
	return (V_0*a_r)/r;
end

# ╔═╡ 4e92a440-6bc3-11eb-108c-e1cb845eb432
r_vec = 0:0.01:10

# ╔═╡ 7cf30c32-6bc3-11eb-1e1d-6d3a80c175f6
y = V.(r_vec)

# ╔═╡ 995855b0-6bc3-11eb-013d-29ec57548faf
plot(r_vec, y, xlabel="Radius(r)", ylabel="Volatge(V)", ylim=(0, V_0+1), label="V")

# ╔═╡ Cell order:
# ╠═53c7d9b0-6b89-11eb-398b-675411f493dc
# ╠═15c25e30-6bbe-11eb-2e55-cfe5abc44d2b
# ╟─951df740-6b8a-11eb-2caa-b5a1291edc31
# ╠═f87b9c1e-6d6f-11eb-0805-0923965aba3b
# ╠═5a5fb850-6b8f-11eb-2d2c-8dcedf10efd9
# ╟─1ed77c60-6bc5-11eb-1a68-e9e8d51f7d87
# ╟─8fc3f460-6bc7-11eb-26b8-61d35bfa4fef
# ╠═f0717252-6bc8-11eb-3bfd-d5d0d0fd299c
# ╠═660d2082-6b89-11eb-100d-8b577e9ae037
# ╠═6a2b5b50-6b89-11eb-38e5-33a649d074ae
# ╠═7418b540-6b89-11eb-1f53-89165d55d825
# ╟─78ff2210-6b89-11eb-29df-c5714c003c7d
# ╟─6d19078e-6b98-11eb-265b-c913b44ff56a
# ╠═eaee7830-6b98-11eb-0a38-3b3060407cc1
# ╠═0de35890-6bc3-11eb-2e29-0f5ebe681620
# ╠═4e92a440-6bc3-11eb-108c-e1cb845eb432
# ╠═7cf30c32-6bc3-11eb-1e1d-6d3a80c175f6
# ╟─995855b0-6bc3-11eb-013d-29ec57548faf

max_stars_repo_name: tskarhed/julia-notebooks

max_issues_repo_name: tskarhed/julia-notebooks

max_forks_repo_name: tskarhed/julia-notebooks


text:
import StatisticalNonlocality: ou_transition_matrix, uniform_phase
import StatisticalNonlocality: discrete_laplacian

@testset "Transition Matrix Column Sum Zero" begin
    for n in 1:10
        M = ou_transition_matrix(n)
        @test all(abs.(sum(M, dims = 1)) .≤ eps(n / 2))
        if n > 1
            M = uniform_phase(n)
            @test all(abs.(sum(M, dims = 1)) .≤ eps(n / 2))
            Δ = discrete_laplacian(n)
            @test all(abs.(sum(Δ, dims = 1)) .≤ eps(n^2 / 2))
        end
    end
end

@testset "OU Transition Matrix Correctness Test for n in 1:2" begin
    n = 1
    M = ou_transition_matrix(n)
    Mh = [
        -1/2 1/2
        1/2 -1/2
    ]
    @test all(abs.(Mh - M) .≤ eps(n / 2))

    n = 2
    M = ou_transition_matrix(n)
    Mh = [
        -1 1/2 0
        1 -1 1
        0 1/2 -1
    ]
    @test all(abs.(Mh - M) .≤ eps(n / 2))
end

@testset "Uniform Phase Transition Matrix Correctness Test for n = 4" begin
    n = 4
    M = uniform_phase(n)
    Mh = [-1 1 0 0; 0 -1 1 0; 0 0 -1 1; 1 0 0 -1]
    @test all(abs.(Mh - M) .≤ eps(n / 2))

end

@testset "Eigenvalue Test for n in 1:10" begin
    for n in 1:10
        λ_exact = collect((-n):0)
        λ_approx = eigvals(ou_transition_matrix(n))
        @test all(abs.(λ_exact - λ_approx) .≤ eps(20.0 * n))
    end
end

@testset "Symmetry Test for Discrete Laplacian" begin
    for n in 4:10
        Δ = discrete_laplacian(n)
        @test norm(Δ - Δ') ≤ eps(norm(Δ))
    end
end

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
#######################################################################
#
# An example of an Excel chart with patterns using XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_pattern.xlsx")
	ws = add_worksheet!(wb)
	# Formats used in the workbook.
	bold = add_format!(wb, Dict("bold"=> 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Shingle", "Brick"]
	data = [
		[105 150 130 90 ];
		[50  120 100 110];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])

	# Create a new Chart object.
	chart = add_chart!(wb, Dict("type"=> "column"))

	# Configure the charts. Add two series with patterns. The gap is used to make
	# the patterns more visible.
	add_series!(chart, Dict(
		"name"=>   "=Sheet1!\$A\$1",
		"values"=> "=Sheet1!\$A\$2:\$A\$5",
		"pattern"=> Dict(
			"pattern"=>  "shingle",
			"fg_color"=> "#804000",
			"bg_color"=> "#c68c53"
		),
		"border"=>  Dict("color"=> "#804000"),
		"gap"=>     70,
	))

	add_series!(chart, Dict(
		"name"=>   "=Sheet1!\$B\$1",
		"values"=> "=Sheet1!\$B\$2:\$B\$5",
		"pattern"=> Dict(
			"pattern"=>  "horizontal_brick",
			"fg_color"=> "#b30000",
			"bg_color"=> "#ff6666"
		),
		"border"=>  Dict("color"=> "#b30000"),
	))

	# Add a chart title and some axis labels.
	set_title!(chart, Dict("name"=> "Cladding types"))
	set_x_axis!(chart, Dict("name"=> "Region"))
	set_y_axis!(chart, Dict("name"=> "Number of houses"))

	# Insert the chart into the worksheet.
	insert_chart!(ws, "D2", chart)

	close(wb)
	isfile("chart_pattern.xlsx")
end

test()


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to fill data in the table that will be printed.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Compute the vectors that will be used to fill the rows of the matrix that will
# be printed. This is required due to the different crop modes available.
function _compute_row_fill_vectors(id_rows::Vector{Int},
                                   num_printed_rows::Int,
                                   vcrop_mode::Symbol)

    # Compute the array separation if the vertical crop mode is `:middle`.
    num_rows = length(id_rows)
    Δ₀ = cld(num_printed_rows, 2)

    # Create the vector that will be used to fill the rows given the crop mode.
    jvec  = Vector{Int}(undef, num_printed_rows)
    jrvec = Vector{Int}(undef, num_printed_rows)

    if vcrop_mode == :bottom
        for i = 1:num_printed_rows
            jvec[i]  = i
            jrvec[i] = id_rows[i]
        end
    elseif vcrop_mode == :middle
        i  = 1
        i₀ = 1
        i₁ = 1

        len₀ = Δ₀
        len₁ = num_printed_rows - Δ₀

        while (i₀ ≤ len₀) || (i₁ ≤ len₁)
            if i₀ ≤ len₀
                j        = i₀
                jvec[i]  = j
                jrvec[i] = id_rows[j]

                i  += 1
                i₀ += 1
            end

            if i₁ ≤ len₁
                j = num_printed_rows - i₁ + 1
                jvec[i]  = j
                jrvec[i] = id_rows[num_rows - i₁ + 1]

                i  += 1
                i₁ += 1
            end
        end
    else
        error("Unknown vertical crop mode.")
    end

    return jvec, jrvec
end

# Fill the header and matrix data.
function _fill_matrix_data!(header_str::Matrix{String},
                            header_len::Matrix{Int},
                            data_str::Matrix{Vector{String}},
                            data_len::Matrix{Vector{Int}},
                            cols_width::Vector{Int},
                            id_cols::Vector{Int},
                            id_rows::Vector{Int},
                            num_lines_in_row::Vector{Int},
                            jvec::Vector{Int},
                            jrvec::Vector{Int},
                            Δc::Int,
                            columns_width::Vector{Int},
                            data::Any,
                            header::Any,
                            formatters::Tuple,
                            display::Display,
                            # Configuration options.
                            autowrap::Bool,
                            cell_first_line_only::Bool,
                            compact_printing::Bool,
                            crop_subheader::Bool,
                            fixed_col_width::Vector{Bool},
                            linebreaks::Bool,
                            maximum_columns_width::Vector{Int},
                            noheader::Bool,
                            renderer::Union{Val{:print}, Val{:show}},
                            vcrop_mode::Symbol)

    num_printed_rows, num_printed_cols = size(data_str)
    header_num_rows, ~ = size(header_len)

    # This variable stores the predicted table width. If the user wants
    # horizontal cropping, then it can be use to avoid unnecessary processing of
    # columns that will not be displayed.
    pred_tab_width = 0

    for i = (1+Δc):num_printed_cols
        # Here we store the number of processed rows. This is used to save
        # processing if the user wants to crop the output and has cells with
        # multiple lines.
        num_processed_rows = 0

        # Index of the i-th printed column in `data`.
        ic = id_cols[i-Δc]

        if !noheader
            for j = 1:header_num_rows
                id = (ic-1)*header_num_rows + j
                header_ij = isassigned(header,id) ? header[id] : undef

                # NOTE: For headers, we always use `print` instead of `show` to
                # avoid quotes.
                hstr, hlstr, cell_width =
                    _parse_cell_text(header_ij;
                                     autowrap = false,
                                     cell_first_line_only = false,
                                     column_width = -1,
                                     compact_printing = compact_printing,
                                     has_color = display.has_color,
                                     linebreaks = false,
                                     renderer = Val(:print))

                header_str[j,i] = first(hstr)
                header_len[j,i] = first(hlstr)

                num_processed_rows += 1

                # If the user does not want a fixed column width, then we must
                # store the information to automatically compute the field size.
                if !fixed_col_width[ic]
                    # Check if we should consider the sub-header size when
                    # computing the column width.
                    if (j == 1) || (!crop_subheader)
                        # Check if we need to increase the columns size.
                        cols_width[i] < cell_width && (cols_width[i] = cell_width)

                        # Make sure that the maximum column width is respected.
                        if maximum_columns_width[ic] > 0 &&
                            maximum_columns_width[ic] < cols_width[i]
                            cols_width[i] = maximum_columns_width[ic]
                        end
                    end
                end
            end
        end

        for k = 1:num_printed_rows
            j  = jvec[k]
            jr = jrvec[k]

            # Apply the formatters.
            data_ij = isassigned(data,jr,ic) ? data[jr,ic] : undef

            data_ij_type = typeof(data_ij)

            for f in formatters
                data_ij = f(data_ij, jr, ic)
            end

            # Parse the cell.
            data_str[j,i], data_len[j,i], cell_width =
                _parse_cell_text(data_ij;
                                 autowrap = autowrap && fixed_col_width[ic],
                                 cell_data_type = data_ij_type,
                                 cell_first_line_only = cell_first_line_only,
                                 column_width = columns_width[ic],
                                 compact_printing = compact_printing,
                                 has_color = display.has_color,
                                 linebreaks = linebreaks,
                                 renderer = renderer)

            # Check if we must update the number of lines in this row.
            num_lines_ij = length(data_str[j,i])
            num_lines_in_row[j] < num_lines_ij && (num_lines_in_row[j] = num_lines_ij)

            num_processed_rows += num_lines_ij

            # If the user does not want a fixed column width, then we must store
            # the information to automatically compute the field size.
            if !fixed_col_width[ic]
                # Check if we need to increase the columns size.
                cols_width[i] < cell_width && (cols_width[i] = cell_width)

                # Make sure that the maximum column width is respected.
                if maximum_columns_width[ic] > 0 && maximum_columns_width[ic] < cols_width[i]
                    cols_width[i] = maximum_columns_width[ic]
                end
            end

            # If the crop mode if `:middle`, then we need to always process a
            # row in the top and in another in the bottom before stopping due to
            # display size. This is required to avoid printing from a cell that
            # is undefined. Notice that due to the printing order in `jvec` we
            # just need to check if `k` is even.
            if ( (vcrop_mode == :bottom) ||
                 ( (vcrop_mode == :middle) && (k % 2 == 0) ) ) &&
                (display.size[1] > 0) && (num_processed_rows ≥ display.size[1])
                break
            end
        end

        # If the user horizontal cropping, then check if we need to process
        # another column.
        #
        # TODO: Should we take into account the dividers?
        if display.size[2] > 0
            pred_tab_width += cols_width[i]

            if pred_tab_width > display.size[2]
                num_printed_cols = i
                break
            end
        end
    end

    return num_printed_cols, num_printed_rows
end

# Fill the information related to the row number column.
function _fill_row_number_column!(header_str::Matrix{String},
                                  header_len::Matrix{Int},
                                  data_str::Matrix{Vector{String}},
                                  data_len::Matrix{Vector{Int}},
                                  cols_width::Vector{Int},
                                  id_rows::Vector{Int},
                                  jvec::Vector{Int},
                                  jrvec::Vector{Int},
                                  noheader::Bool,
                                  num_rows::Int,
                                  row_number_column_title::String)

    num_printed_rows = size(data_str)[1]

    num_printed_rows == 0 && return nothing

    # Set the header of the row column.
    header_str[1,1]      = row_number_column_title
    header_str[2:end,1] .= ""
    header_len[1,1]      = textwidth(row_number_column_title)

    # Set the data of the row column.
    for i = 1:num_printed_rows
        j  = jvec[i]
        jr = jrvec[i]

        data_str[j,1] = [string(jr)]

        # Here we can use `length` because there will be no UTF-8 character
        # in this cell.
        data_len[j,1] = [length(data_str[j,1][1])]
    end

    # The row number width depends on how many digits the total number of
    # rows has and the length of the header. Notice that if `noheader` is
    # set to `true`, then we should not take the word "Row" into account.
    cols_width[1] = max(noheader ? 0 : header_len[1,1],
                        floor(Int, log10(num_rows)) + 1)

    return nothing
end

# Fill the information related to the row name column.
function _fill_row_name_column!(header_str::Matrix{String},
                                header_len::Matrix{Int},
                                data_str::Matrix{Vector{String}},
                                data_len::Matrix{Vector{Int}},
                                cols_width::Vector{Int},
                                row_names::AbstractVector,
                                jvec::Vector{Int},
                                jrvec::Vector{Int},
                                Δc::Int,
                                compact_printing::Bool,
                                renderer::Union{Val{:print}, Val{:show}},
                                row_name_column_title::String)

    num_printed_rows = size(data_str)[1]

    num_printed_rows == 0 && return nothing

    # Escape the row name column title.
    header_str[1,Δc]      = first(_render_text(Val(:print),
                                               row_name_column_title,
                                               compact_printing = compact_printing,
                                               linebreaks = false))
    header_str[2:end,Δc] .= ""

    # Compute the length of the row name column title.
    str_len = textwidth(header_str[1,Δc])
    header_len[1,Δc] = str_len

    # Convert the row names to string.
    max_size = 0
    for i = 1:num_printed_rows
        j  = jvec[i]
        jr = jrvec[i]

        row_names_j = isassigned(row_names,jr) ? row_names[jr] : undef
        row_name_str, row_name_lstr, cell_width =
            _parse_cell_text(row_names_j;
                             autowrap = false,
                             cell_first_line_only = false,
                             column_width = -1,
                             compact_printing = compact_printing,
                             linebreaks = false,
                             renderer = Val(:print))

        data_str[j,Δc] = row_name_str
        data_len[j,Δc] = row_name_lstr

        cell_width > max_size && (max_size = cell_width)
    end

    # Obtain the size of the row name column.
    cols_width[Δc] = max(str_len, max_size)

    return nothing
end


text:
using Hello;
Hello.command_main();


text:
"""
process padding in shared memory - where we have results that may affect other blocks
"""
module ProcessPadding
using CUDA, ..CUDAGpuUtils, Logging

"""
Analizing all padding planes - so we need properly pass mainly x,y,z indexes ...

x,y,z - this are coordinates of upper left corner of our data block
shmem- shared memory boolean array that will be used for reduction as we need to establish weather neighbouring block needs to be activated or not
    - we are reusing main shared memory array
nextBlockXchange,nextBlockYchange,nextBlockZchange - indicates wheather we need to add or subtracct from current block ids  in order to get to the  block that is neighbouring from this padding side that we are analyzing
currBlockX, Y, Z - points to position of the current block
sourceArray,referenceArray - global memory arrays with primary data  needed to establish weather we need to write data there
resArr - 3 dimensional array where we put results
metaData - basic data about data blocks - stored in 3 dimensional array
metadataDims - dimensions of metaData array
isPassGold - true if we are dilatating currently gold standard mask false if we are dilatating other mask
resArraysCounter - counter needed to add keep track on how many results we have in our result array
                    - we have separate counter for gold pass and other pass
"""
function processAllPaddingPlanes(    x::UInt16
                                    ,y::UInt16
                                    ,z::UInt16
                                    ,shmem
                                    ,currBlockX::UInt16
                                    ,currBlockY::UInt16
                                    ,currBlockZ::UInt16
                                    ,sourceArray
                                    ,referenceArray
                                    ,resArr
                                    ,metaData
                                    ,metadataDims
                                    ,isPassGold::Bool
                                    ,locArr
                                    ,resArraysCounter
                                    ,iterationNumber)
   #neded to clear memory for establishing weather we should process given pading or not
    #- basically wheather we care about block next to this - so is it full or is it on edge ...
    locArr = Int32(0)

   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,false,false,false,false,false,true,1,shmem,metadataDims,isPassGold)
   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,false,false,true,false,false,false,3 ,shmem,metadataDims,isPassGold)

   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,false,true,false,false,false,false,5,shmem,metadataDims,isPassGold)
   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,false,false,false,false,true,false,7 ,shmem,metadataDims,isPassGold)

   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,false,false,false,true,false,false,9,shmem,metadataDims,isPassGold)
   isNextBlockOfIntrest(currBlockX,currBlockY,currBlockZ,true,false,false,false,false,false,11,shmem,metadataDims,isPassGold)

    for i in [UInt8(1),UInt8(3),UInt8(5),UInt8(7),UInt8(9),UInt8(11)]
        #we need to check weather we are in the same warp as we used in isNextBlockOfIntrest
        if(threadIdxY()==i )
            shmem[i,i,13]=reduce_warp_and(locArr[i], UInt8(32))
        end#if
    end#for
    sync_threads()#now we should have all required booleans and we can reduce them

   # so we have 6 planes to analyze all around our data cube and we already collected data wheather are blocks adjacent to those planes are of our intrest

   processsPaddingTOP(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
   processsPaddingBOTTOM(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)

   processsPaddingFront(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
   processsPaddingBack(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)

   processsPaddingLeft(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
   processsPaddingRight(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)

 end#processAllPaddingPlanes

"""
this is a little bit tricky to get it correctly as we are looking at paddings from all sides
    - this is also crucial from perspective how we should access indexes of target although we are looking for paddings from all directions
we need to pass data to source array about new "trues"
    paddingVal - value of shared memory associated with this lane
    correctedX ,y,z - coordinates of the source arrays - important we need to set those smartly as we are looking at planes all around
    shmem- shared memory boolean array that will be used for reduction as we need to establish weather neighbouring block needs to be activated or not
        - we are reusing main shared memory array
    nextBlockXchange,nextBlockYchange,nextBlockZchange - indicates wheather we need to add or subtracct from current block ids  in order to get to the  block that is neighbouring from this padding side that we are analyzing
    currBlockX, Y, Z - points to position of the current block
    sourceArray,referenceArray - global memory arrays with primary data  needed to establish weather we need to write data there
    resArr - 3 dimensional array where we put results
    sliceNumbManual - controlls what slice of shared memory we will use
    metaData - basic data about data blocks - stored in 3 dimensional array
    metadataDims - dimensions of metaData array
    isPassGold - true if we are dilatating currently gold standard mask false if we are dilatating other mask
"""
function processPaddingPlane(paddingVal::Bool
                            ,isOfIntrest::Bool
                            ,correctedX::UInt16
                            ,correctedY::UInt16
                            ,correctedZ::UInt16
                            ,nextBlockXIncrease::Bool
                            ,nextBlockYIncrease::Bool
                            ,nextBlockZIncrease::Bool
                            ,nextBlockXDecrease::Bool
                            ,nextBlockYDecrease::Bool
                            ,nextBlockZDecrease::Bool
                            ,sliceNumbManual::UInt8
                            ,shmem
                            ,currBlockX::UInt16
                            ,currBlockY::UInt16
                            ,currBlockZ::UInt16
                            ,sourceArray
                            ,referenceArray
                            ,resArr
                            ,metaData
                            ,metadataDims
                            ,isPassGold::Bool
                            ,mainQuesCounter
                            ,mainWorkQueue
                            ,resArraysCounter
                            ,iterationNumber)
    #we can ignore all of this if this is edge block - and we can ignore paddings in this situations from those directions, the same for blocks that are full
    newBlockX = currBlockX+nextBlockXIncrease-nextBlockXDecrease
    newBlockY= currBlockY+nextBlockYIncrease-nextBlockYDecrease
    newBlockZ= currBlockZ+nextBlockZIncrease-nextBlockZDecrease

        #before we will start processing padding we need to check weather it makes sense - weather we are not at the edge  or next block is full
        #we need to check is it there  sth of our intrest there - we check shared memory value established in isNextBlockOfIntrestHelper function
        if( shmem[sliceNumbManual,sliceNumbManual,13]  )

            #we are intrested in futher processing only if we have some true here
            if(paddingVal)
                setGlobalsFromPadding(correctedX,correctedY,correctedZ,resArr,sourceArray,resArraysCounter, iterationNumber)
            end#if
            #we need to reduce now  the values  of padding vals to establish weather there is any true there if yes we put the neighbour block to be active
            reducePaddingPlane(shmem,paddingVal,sliceNumbManual )
            #we have true in shmem[1,1,sliceNumbManual+1] if there is at least one true in this plane and false otherwise
            if(shmem[1,1,sliceNumbManual+1])
                activateNextBlock(newBlockX,newBlockY,newBlockZ, metaData,metadataDims,isPassGold,mainQuesCounter,mainWorkQueue)
            end#if

        end#if

end#processPaddingPlane

"""
    activate neighbour  block - invoked when in padding plane we have at least one true
    metaData - basic data about data blocks - stored in 3 dimensional array
    metadataDims - dimensions of metaData array
    isPassGold - true if we are dilatating currently gold standard mask false if we are dilatating other mask
    mainQuesCounter - counter that we will update atomically and will be usefull to populate the work queue
mainWorkQueue - the list of the indicies of  data blocks in metadata with additional information is it referencing the goldpass or second one

"""
function activateNextBlock(newBlockX,newBlockY,newBlockZ,metadataDims,metaData,isPassGold,mainQuesCounter,mainWorkQueue)
    if(threadIdxY()==1&& threadIdxX()==1)
    metaData[newBlockX,newBlockY,newBlockZ,isPassGold+1]= true
    end
   #so in case it not empty and not full we need to put it into the work queue and increment appropriate counter
    if(threadIdxY()==3&& threadIdxX()==3)
        mainWorkQueue[CUDA.atomic_inc!(pointer(mainQuesCounter), UInt16(1))+1,:]=[newBlockX,newBlockY,newBlockZ,UInt8(isPassGold)] #x,y,z dim of block in metadata
    end#if

end#activateNextBlock

"""
we need to reduce now  the values  of padding vals to establish weather there is any true there if yes we put the neighbour block to be active
    if all goes well it should write true to shmem[1,1,sliceNumbManual+1] if there is at leas one true in this plane and false otherwise
"""
function reducePaddingPlane(shmem,paddingVal,sliceNumbManual )::Bool
    #to be experimented is it better to reduce to varyiong y or x
    @inbounds shmem[1,threadIdxY(),sliceNumbManual]=  reduce_warp_or(paddingVal, UInt8(32))
    #so now we have 32 booleans in shared memory so we need to reduce it one more time using single warp
    #TODO() check weather warps are column or row wise this below also I am not sure is it good
    if(threadIdxY()==1  && threadIdxX()==1  )
        @inbounds  shmem[1,1,sliceNumbManual+1]= reduce_warp_or(shmem[1,threadIdxY(),sliceNumbManual], UInt8(32))
    end
end

"""
accesses source arrays and modifies it if needed - invoked from procesing of padding
    sourceArrValue - value of source array in x,y,z position
    generally function will be invoked only if value in padding was set to true

resArraysCounter - counter needed to add keep track on how many results we have in our result array
                    - we have separate counter for gold pass and other pass

"""
function setGlobalsFromPadding(correctedX::UInt16,correctedY::UInt16,correctedZ::UInt16,resArr,sourceArray,resArraysCounter, iterationNumber)
       processResArrBasedOnoldVal(CUDA.atomic_inc!(pointer(resArr[correctedX,correctedY,correctedZ]), UInt16(iterationNumber)),resArraysCounter )
       sourceArray[correctedX,correctedY,correctedZ]= true
end#setGlobalsFromPadding

"""
when we set new result from padding we need to take into account possibility that neighbour block already did it
hence when we set atomic we check the old value if old value was 0 all is good if not we  do not increse  the rescounter
"""
function processResArrBasedOnoldVal(oldVal,resArraysCounter  )
    if(oldVal==0)
        CUDA.atomic_inc!(pointer(resArraysCounter), UInt16(1))
    end
end

"""
process Top Padding
"""
function processsPaddingTOP(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[threadIdxX(), threadIdxY(),1],x,y,z-1 # to specialize
       ,false,false,false,false,false,true # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,1 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold ,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)# this arguments will be the same for all invocations
end

"""
process Bottom Padding
"""
function processsPaddingBOTTOM(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[threadIdxX(), threadIdxY(),34],x,y,z+34
       ,false,false,true,false,false,false # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,3 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,mainQuesCounter,mainWorkQueue ,resArraysCounter,iterationNumber)# this arguments will be the same for all invocations
end

"""
process  Padding Front
"""
function processsPaddingFront(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[threadIdxX(), 34 ,threadIdxY()],x,y+34,z
       ,false,true,false,false,false,false # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,5 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,mainQuesCounter,mainWorkQueue ,resArraysCounter,iterationNumber)# this arguments will be the same for all invocations
end

"""
process  Padding back
"""
function processsPaddingBack(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[threadIdxX(), 1 ,threadIdxY()],x,y-1,z
       ,false,false,false,false,true,false # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,7 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold ,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)# this arguments will be the same for all invocations
end

"""
process  Padding left
"""
function processsPaddingLeft(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[1,threadIdxX(),threadIdxY()],x-1,y,z
       ,false,false,false,true,false,false # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,9 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,mainQuesCounter,mainWorkQueue ,resArraysCounter,iterationNumber)# this arguments will be the same for all invocations
end

"""
process  Padding right
"""
function processsPaddingRight(shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,x,y,z,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber)
    processPaddingPlane(
        sourceArray[34,threadIdxX(),threadIdxY()],x+34,y,z
       ,true,false,false,false,false,false # nextBlockXIncrease,nextBlockYIncrease,nextBlockZIncrease,nextBlockXDecrease,nextBlockYDecrease,nextBlockZDecrease
       ,11 # to specialize
      ,shmem,currBlockX,currBlockY,currBlockZ,sourceArray,referenceArray,resArr ,metaData,metadataDims,isPassGold,mainQuesCounter,mainWorkQueue,resArraysCounter,iterationNumber )# this arguments will be the same for all invocations
end

end#ProcessPadding

# """
#     before we will start processing padding we need to check weather it makes sense - weather we are not at the edge  or next block is full
#     we will aslo experiment with looking in given direction if all block in this direction are full or empty in both masks - sth to experiment upon
#     In order to parallelize the work we will  use separate warps to check the conditions for all
#     we will write results to appropriate spots in shared memory
# """
# function isNextBlockOfIntrest(currBlockX::UInt16
#                             ,currBlockY::UInt16
#                             ,currBlockZ::UInt16
#                             ,nextBlockXIncrease::Bool
#                             ,nextBlockYIncrease::Bool
#                             ,nextBlockZIncrease::Bool
#                             ,nextBlockXDecrease::Bool
#                             ,nextBlockYDecrease::Bool
#                             ,nextBlockZDecrease::Bool
#                             ,sliceNumbManual::UInt8
#                             ,shmem
#                             ,metadataDims
#                             ,isPassGold)
#        nextBlockX = currBlockX+nextBlockXIncrease-nextBlockXDecrease
#        nextBlockY = currBlockY+nextBlockYIncrease-nextBlockYDecrease
#        nextBlockZ = currBlockZ+nextBlockZIncrease-nextBlockZDecrease

# #we are selecting single warp
# if(threadIdxY()==sliceNumbManual )
#     #we need to check is it there at all - it will not if we are in border case
#     if(threadIdxX()==1)
#         locArr[sliceNumbManual]= (nextBlockX>0)
#     elseif(threadIdxX()==2)
#         locArr[sliceNumbManual]= (nextBlockX<=metadataDims[1])
#     elseif(threadIdxX()==3)
#         locArr[sliceNumbManual]= (nextBlockY>0)
#     elseif(threadIdxX()==4)
#         locArr[sliceNumbManual]= (nextBlockY<=metadataDims[2])
#     elseif(threadIdxX()==5)
#         locArr[sliceNumbManual]= (nextBlockZ>0)
#     elseif(threadIdxX()==6)
#         locArr[sliceNumbManual]= (nextBlockZ<=metadataDims[3])
#     elseif(threadIdxX()==7)
#         locArr[sliceNumbManual] = metaData[newBlockX,newBlockY,newBlockZ,isPassGold+3] #then we need to check weather mask is already full - in this case we can not activate it
#     else
#         locArr[sliceNumbManual] = true # so when all are true we know it
#     end#inner if

# end#outer if

# #we will aslo experiment with looking in given direction if all block in this direction are full or empty in both masks we can also


text:
using EffectiveWaves

using JLD

include("optimise_wavenumber.jl")

medium = Acoustic(2; ρ=1.0,c=1.0)
ωs = 0.01:0.01:1.0
MaxTime = 10000.0

function f_slow1(sps::Species, medium::PhysicalMedium,ωs,kTs)
  N = Int(round(length(ωs)/2))
  mean(imag(kTs[1:N])) + mean(ωs[1:N]./real(kTs[1:N]))
end
function f_slow2(sps::Species, medium::PhysicalMedium,ωs,kTs)
  N = Int(round(length(ωs)/2))
  mean(imag(kTs[N:end])) + mean(ωs[N:end]./real(kTs[N:end]))
end
function f_fast1(sps::Species, medium::PhysicalMedium,ωs,kTs)
  N = Int(round(length(ωs)/2))
  mean(imag(kTs[1:N])) - mean(ωs[1:N]./real(kTs[1:N]))
end
function f_fast2(sps::Species, medium::PhysicalMedium,ωs,kTs)
  N = Int(round(length(ωs)/2))
  mean(imag(kTs[N:end])) - mean(ωs[N:end]./real(kTs[N:end]))
end

fs = [f_slow1,f_slow2,f_fast1,f_fast2]

for f in fs
  species = optimal_species(f,medium,ωs; fix_fields = [(:c,[1.0+0.0im,1.0+0.0im])],
              opt_fields = [(:r,(0.01,2.0)),(:volfrac,(0.,0.14)),(:ρ,(0.0,100.))],
              num_species=2, MaxTime=MaxTime, method=:xnes)
  save("data/$f.jld","objective", string(f), "species", species, "ωs", ωs, "medium", medium)
end

# using Plots
# plotly()
# strfs = ["f_slow1","f_slow2","f_fast1","f_fast2"];
# species = load("data/$f.jld")["species"]
#
# kTs_arr = [
#     [ sqrt(wavenumber_low_volumefraction(ω,medium, sps)) for ω in ωs]
# for sps in [[species[1]],[species[2]],species]];
#
# speed_arr = [ ωs./real(kTs) for kTs in kTs_arr]
# atten_arr = imag(kTs_arr)
#
# labs = ["void" "stone" "mix"];
# p1 = plot(ωs./medium.c,speed_arr, xlabel="k", ylabel="wave speed", labels=labs,);
# p2 = plot(ωs./medium.c,atten_arr, xlabel="k", ylabel="attenuation", labels=labs,);
# plot(p1,p2)

max_stars_repo_name: arturgower/EffectiveWaves


text:
@testset "ast" begin
    @testset "Base.show overload" begin
        ast = JSExpr.@crawl foo = bar
        @test sprint(
            io -> Base.show(io, MIME("text/plain"), ast),
        ) == """JSAST(:(=), JSTerminal(foo), JSTerminal(bar))"""
    end
end


text:
logfk(m::OverlappingDiscs, p, xbefore) = log(fk(m, p, xbefore))
function fk(m::OverlappingDiscs, p, xbefore)
  m.theta(countoverlaps(m, p, xbefore), length(xbefore))
end
function countoverlaps(m::OverlappingDiscs, p, xbefore)
  k = 0
  x, y = p
  for (i, (x1, y1)) in enumerate(xbefore)
    if (x1-x)^2 + (y1-y)^2 <= m.R(i)^2
      k += 1
    end
  end
  k
end
function add1tocircle!(image, ks, x, y, R)
  add1tocircle!(image, ks, round(Int, x), round(Int, y), R)
end
function add1tocircle!(image, ks, x::Int, y::Int, R)
  Ri = ceil(Int, R) + 1
  xi = round(Int, x)
  yi = round(Int, y)
  checkbounds(image, xi, yi)
  xr = intersect(xi .+ (-Ri:Ri), axes(image, 1)) .- xi
  yr = intersect(yi .+ (-Ri:Ri), axes(image, 2)) .- yi
  #println((x, y, xr, yr, R, Ri, axes(image, 1), axes(image, 2)))
  for dx in xr, dy in yr
    if dx^2 + dy^2 <= R^2
      k = image[x+dx, y+dy]
      ks[k+1] -= 1
      k += 1
      ks[k+1] += 1
      image[x+dx, y+dy] = k
    end
  end
end

function compute_integral(m::OverlappingDiscs, data, scale, nx, ny)
  cellsize = 1/scale^2
  image = zeros(Int, nx, ny)
  ks = zeros(Int, length(data) + 1)
  # Remember that arrays have 1 based indexing whereas counts start from 0
  ks[1] = length(image)
  I = Vector{Float64}(undef, length(data))
  for (i, (x, y)) in enumerate(data)
    add1tocircle!(image, ks, 0.5 + x*scale, 0.5 + y*scale, m.R(i)*scale)
    Ii = 0.0
    for (ki, k) in enumerate(ks)
      Ii += k * m.theta(ki-1, i)
    end
    I[i] = Ii * cellsize
  end
  I
end

function discretization(window, nx::Int)
  (x1,x2), (y1,y2) = window.x, window.y
  ny = round(Int, nx*(y2-y1)/(x2-x1))
  scale = nx/(x2-x1)
  scale, x1, y1, nx, ny
end

function normconstants(m::OverlappingDiscs, xs, window, int)
  int.parallel === false || throw(ArgumentError("Model doesn't support parallization"))
  nx = int.nx
  compute_integral(m, xs, discretization(window, nx)...)
end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses