text:
using ChanceConstraintExtensions
import Ipopt
using JuMP
using Test

@testset "Γ" begin
    @test ChanceConstraintExtensions.Γ(-1.1, 1.0) == 1.0
    @test ChanceConstraintExtensions.Γ(-0.11, 0.1) == 1.0
    @test ChanceConstraintExtensions.Γ(0.11, 0.1) == 0.0
    @test ChanceConstraintExtensions.Γ(1.1, 1.0) == 0.0
    @test ChanceConstraintExtensions.Γ(0.0, 0.1) ≈ 0.5
    @test ChanceConstraintExtensions.Γ(0.0, 0.2) ≈ 0.5
    @test ChanceConstraintExtensions.Γ(0.05, 0.1) ≈ 53 / 512
    @test ChanceConstraintExtensions.Γ(-0.05, 0.1) ≈ 1 - 53 / 512
end

@testset "Example" begin
    Ω = [-1.0, -0.2, 0.2, 1.0]
    K = length(Ω)
    model = Model(with_optimizer(Ipopt.Optimizer))
    set_parameter(model, "print_level", 0)
    @variable(model, -2 <= x <= 2, start = 0.4)
    chance_constraint(
        model,
        @NLexpression(model, [k = 1:K], x^2 - 2 + Ω[k]),
        probability = 0.3, smoothing = 0.01, bisection_atol=1e-4
    )
    @objective(model, Max, x)
    optimize!(model)
    @test value(x) ≈ 1.34035 atol = 1e-4
end


text:
using NEAT
using Cambrian
using Test

function test_identical(ind::NEATInd, ind2::NEATInd)
    @test ind.n_in == ind2.n_in
    @test ind.n_out == ind2.n_out
    @test ind.n_hidden == ind2.n_hidden

    for i in eachindex(ind.connections)
        @test ind.connections[i].in_node == ind2.connections[i].in_node
        @test ind.connections[i].out_node == ind2.connections[i].out_node
        @test ind.connections[i].weight == ind2.connections[i].weight
        @test ind.connections[i].enabled == ind2.connections[i].enabled
        @test ind.connections[i].innovation == ind2.connections[i].innovation
    end

    @test all(ind.fitness .== ind2.fitness)
end

include("individual.jl")
include("process.jl")
include("mutation.jl")

ext: jl


text:
# Tests for the Decision stump classifier and regressor
function t_stump()

Ac = j4pr.DataGenerator.fish(20)  # classification dataset
Ar = j4pr.DataGenerator.fishr(20) # regression dataset

tol  = 1e-6; # tolerance when comparing results
Wclass = [j4pr.stump(vartypes=v, nthresh=n, split=s, count=c, prop=p, crit=cr ) for
	  						v in [:nominal, :real, Dict(1=>:nominal)],
	  						n in [2,5,100],
							s in [:linear,:density],
							c in [0,10],
							p in [0.0],
							cr in [:gini,:entropy,:misclassification]
]

# Test classification
for w in Wclass
	Test.@test try
		wt = Ac |> w
		result = Ac |> wt
		result2 = +Ac |> wt
		sum(abs.(result.x - result2) .>= tol) > 0 ? false : true
	catch
		false
	end
end

# Test classification on a 1-D dataset
A1d = j4pr.datacell([rand(100);2*rand(100)],[zeros(100);ones(100)])
Test.@test try
	wt1 = A1d |> Wclass[1]
	result = A1d |> wt1
	buf = IOBuffer()
	Base.show(buf,wt1.x.data)
	true
	catch
		false
	end

# Test Prin# model=:linear fit is not tested because it is unstable in certain situations
Wreg = [j4pr.stumpr(model=m, errcrit=ec, vartypes=v, nthresh=n, split=s, count=c, prop=p) for
							m in [:mean, :median],
						        ec in [(x,y)->mean((x-y).^2), (x,y)->mean(abs.(x-y))],
							v in [:nominal, :real, Dict(1=>:nominal)],
	  						n in [2,5],
							s in [:linear,:density],
							c in [0,10],
							p in [0.0]
]

# Test regression
for w in Wreg
	Test.@test try
		wt = Ar |> w

		result = Ar |> wt
		result2 = +Ar |> wt

		sum(abs.(result.x - result2) .>= tol) > 0 ? false : true
	catch
		false
	end
end

# Test regression on a 1-D dataset
A1d = j4pr.datacell([rand(100);2*rand(100)],rand(200))
Test.@test try
	wt1 = A1d |> Wreg[1]
	result = A1d |> wt1
	buf = IOBuffer()
	Base.show(buf,wt1.x.data)
	true
	catch
		false
	end

end


text:
"""
Google API URL roots.
"""
module root

export API_ROOT, SCOPE_ROOT, AUD_ROOT, METADATA_ROOT, isurl

"""
    isurl(path)

Return true if `path` is a URL and false a path fragment.
"""
isurl(path::AbstractString) = occursin(r"^https?://", path)

const API_ROOT = "https://www.googleapis.com"
const SCOPE_ROOT = "$API_ROOT/auth"
const AUD_ROOT = "$API_ROOT/oauth2/v4/token"
const METADATA_ROOT = "http://metadata.google.internal/computeMetadata/v1"

end


text:
k = Kosmo()

mutable struct Health <: Component
    v::Int
end

ply = Entity!(k, Health(100))

system!(k) do x::Health
    x.v = x.v-3
    println(x.v)
end

t = Kosmos.run_systems(k)


text: wait_for_key(prompt) = (print(stdout, prompt); readline(stdin); nothing)


text:
# Single choice variable
dipar = Dict(:nK => 40, :nz => 5, :γ => 2., :F=>0.)
p_1 = createmodel(:CapitalAdjustModel; dipar..., intdim = :Separable)
sol_1 = solve(p_1; monotonicity = false, concavity = false)
sol_1_mon = solve(p_1; monotonicity = true, concavity = false)
sol_1_conc = solve(p_1; monotonicity = false, concavity = true)

# two choice variables
dipar = Dict(:nK => 10, :nN => 7, :nz => 3, :F=>0.)
p_2 = createmodel(:CapitalAdjustModel2; dipar..., intdim = :Separable)
optdict = Dict(:rewardcall=>:jit)
sol_2 = solve(p_2; optdict..., monotonicity=[false,false], concavity=[false,false])
sol_2_mon = solve(p_2; optdict..., monotonicity=[true,true], concavity=[false,false])
sol_2_conc = solve(p_2; optdict..., monotonicity=[false,false], concavity=[true,true])

@testset "monotonicity concavity" begin

    @test isapprox(sol_1, sol_1_mon,  rtol=mytol)
    @test isapprox(sol_1, sol_1_conc, rtol=mytol)
    @test isapprox(sol_2, sol_2_mon, rtol=mytol)
    @test isapprox(sol_2, sol_2_conc, rtol=mytol)

end


text:
using Psychometrics
using Test

include("reliabilitytests.jl")

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


text:
"""
    AbstractBridge

Subtype of [`MathOptInterface.Bridges.AbstractBridge`](@ref) for variable
bridges.
"""
abstract type AbstractBridge <: MOIB.AbstractBridge end

"""
    IndexInVector

Index of variable in vector of variables.
"""
struct IndexInVector
    value::Int
end

"""
    bridge_constrained_variable(BT::Type{<:AbstractBridge}, model::MOI.ModelLike,
                                set::MOI.AbstractSet)

Bridge the constrained variable in `set` using bridge `BT` to `model` and returns
a bridge object of type `BT`. The bridge type `BT` should be a concrete type,
that is, all the type parameters of the bridge should be set. Use
[`concrete_bridge_type`](@ref) to obtain a concrete type for given set types.
"""
function bridge_constrained_variable end

"""
    function MOI.get(model::MOI.ModelLike, attr::MOI.AbstractVariableAttribute,
                     bridge::AbstractBridge)

Return the value of the attribute `attr` of the model `model` for the
variable bridged by `bridge`.
"""
function MOI.get(::MOI.ModelLike, attr::MOI.AbstractVariableAttribute,
                 bridge::AbstractBridge)
    throw(ArgumentError("Variable bridge of type `$(typeof(bridge))` does not support accessing the attribute `$attr`."))
end

"""
    function MOI.get(model::MOI.ModelLike, attr::MOI.AbstractVariableAttribute,
                     bridge::AbstractBridge, i::IndexInVector)

Return the value of the attribute `attr` of the model `model` for the
variable at index `i` in the vector of variables bridged by `bridge`.
"""
function MOI.get(::MOI.ModelLike, attr::MOI.AbstractVariableAttribute,
                 bridge::AbstractBridge, ::IndexInVector)
    throw(ArgumentError("Variable bridge of type `$(typeof(bridge))` does not support accessing the attribute `$attr`."))
end

"""
    supports_constrained_variable(::Type{<:AbstractBridge},
                                   ::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging
constrained variables in `S`.
"""
function supports_constrained_variable(::Type{<:AbstractBridge},
                                        ::Type{<:MOI.AbstractSet})
    return false
end

"""
    added_constrained_variable_types(BT::Type{<:MOI.Bridges.Variable.AbstractBridge},
                                     S::Type{<:MOI.AbstractSet})

Return a list of the types of constrained variables that bridges of type `BT`
add for bridging constrained variabled in `S`. This falls back to
`added_constrained_variable_types(concrete_bridge_type(BT, S))`
so bridges should not implement this.

## Examples

As a variable in [`MathOptInterface.GreaterThan`](@ref) is bridged into
variables in [`MathOptInterface.Nonnegatives`](@ref) by the
[`VectorizeBridge`](@ref):
```jldoctest
BT = MOI.Bridges.Variable.VectorizeBridge{Float64}
S = MOI.GreaterThan{Float64}
MOI.Bridges.added_constrained_variable_types(BT, S)

# output

[(MOI.Nonnegatives,)]
```
"""
function MOIB.added_constrained_variable_types(
    BT::Type{<:AbstractBridge}, S::Type{<:MOI.AbstractSet})
    MOIB.added_constrained_variable_types(concrete_bridge_type(BT, S))
end

"""
    added_constraint_types(BT::Type{<:MOI.Bridges.Variable.AbstractBridge},
                           S::Type{<:MOI.AbstractSet})

Return a list of the types of constraints that bridges of type `BT` add for
for bridging constrained variabled in `S`. This falls back to
`added_constraint_types(concrete_bridge_type(BT, S))`
so bridges should not implement this method.

## Examples

In addition to creating variables in
[`MathOptInterface.PositiveSemidefiniteConeTriangle`](@ref), the
[`RSOCtoPSDBridge`](@ref) also creates
[`MathOptInterface.SingleVariable`](@ref)-in-[`MathOptInterface.EqualTo`](@ref) and
[`MathOptInterface.ScalarAffineFunction`](@ref)-in-[`MathOptInterface.EqualTo`](@ref)
constraints:
```jldoctest
BT = MOI.Bridges.Variable.RSOCtoPSDBridge{Float64}
S = MOI.RotatedSecondOrderCone{
MOI.Bridges.added_constraint_types(BT, S)

# output

[(MOI.SingleVariable, MOI.EqualTo{Float64}), (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})]
```
"""
function MOIB.added_constraint_types(
    BT::Type{<:AbstractBridge}, S::Type{<:MOI.AbstractSet})
    MOIB.added_constraint_types(concrete_bridge_type(BT, S))
end

"""
    concrete_bridge_type(BT::Type{<:AbstractBridge},
                         S::Type{<:MOI.AbstractSet})::DataType

Return the concrete type of the bridge supporting variables in `S` constraints.
This function can only be called if `MOI.supports_constrained_variable(BT, S)`
is `true`.

## Examples

As a variable in [`MathOptInterface.GreaterThan`](@ref) is bridged into
variables in [`MathOptInterface.Nonnegatives`](@ref) by the
[`VectorizeBridge`](@ref),
```jldoctest
BT = MOI.Bridges.Variable.VectorizeBridge{Float64}
S = MOI.GreaterThan{Float64}
MOI.Bridges.Variable.concrete_bridge_type(BT, S)

# output

MOI.Bridges.Variable.VectorizeBridge{Float64,MOI.Nonnegatives}
```
"""
function concrete_bridge_type(bridge_type::DataType,
                              ::Type{<:MOI.AbstractSet})
    return bridge_type
end

function concrete_bridge_type(b::MOIB.AbstractBridgeOptimizer,
                              S::Type{<:MOI.AbstractSet})
    return concrete_bridge_type(MOIB.bridge_type(b, S), S)
end

"""
   unbridged_map(bridge::MOI.Bridges.Variable.AbstractBridge,
                 vi::MOI.VariableIndex)

For a bridged variable in a scalar set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the
bridged variable `vi`.

   unbridged_map(bridge::MOI.Bridges.Variable.AbstractBridge,
                 vis::Vector{MOI.VariableIndex})

For a bridged variable in a vector set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the bridged
variable `vis`. If this method is not implemented, it falls back to calling
the following method for every variable of `vis`.

    unbridged_map(bridge::MOI.Bridges.Variable.AbstractBridge,
                  vi::MOI.VariableIndex, i::IndexInVector)

For a bridged variable in a vector set, return a tuple of pairs mapping the
variables created by the bridge to an affine expression in terms of the bridged
variable `vi` corresponding to the `i`th variable of the vector.

If there is no way to recover the expression in terms of the bridged variable(s)
`vi(s)`, return `nothing`. See [`ZerosBridge`](@ref) for an example of bridge
returning `nothing`.
"""
function unbridged_map end

function unbridged_map(bridge::AbstractBridge, vis::Vector{MOI.VariableIndex})
    mappings = Pair{MOI.VariableIndex, MOI.AbstractScalarFunction}[]
    for (i, vi) in enumerate(vis)
        vi_mappings = unbridged_map(bridge, vi, IndexInVector(i))
        if vi_mappings === nothing
            return nothing
        end
        for mapping in vi_mappings
            push!(mappings, mapping)
        end
    end
    return mappings
end


text:
module ImageGather

    using JUDI, DSP, PyCall

    const dbr = PyNULL()

    function __init__()
        pushfirst!(PyVector(pyimport("sys")."path"),dirname(pathof(ImageGather)))
        copy!(dbr, pyimport("double_rtm"))
    end
    # Utility functions
    include("utils.jl")
    # Surface offset gathers
    include("surface_gather.jl")

end # module


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

using Logging
using MPI
using JLD2
using BenchmarkTools

using Oceananigans
using Oceananigans.Distributed
using Benchmarks

Logging.global_logger(OceananigansLogger())

MPI.Init()

      comm = MPI.COMM_WORLD
local_rank = MPI.Comm_rank(comm)
         R = MPI.Comm_size(comm)

 #assigns one GPU per rank, could increase efficiency but must have enough GPUs
 #CUDA.device!(local_rank)

Nx = parse(Int, ARGS[1])
Ny = parse(Int, ARGS[2])
Rx = parse(Int, ARGS[3])
Ry = parse(Int, ARGS[4])

@assert Rx * Ry == R

@info "Setting up distributed shallow water model with N=($Nx, $Ny) grid points and ranks=($Rx, $Ry) on rank $local_rank..."

topo = (Periodic, Periodic, Flat)
distributed_grid = RegularRectilinearGrid(topology=topo, size=(Nx, Ny), extent=(1, 1))
arch = MultiCPU(grid=distributed_grid, ranks=(Rx, Ry, 1))
model = DistributedShallowWaterModel(architecture=arch, grid=distributed_grid, gravitational_acceleration=1.0)
set!(model, h=1.0)

@info "Warming up distributed shallow water model on rank $local_rank..."

time_step!(model, 1) # warmup

@info "Benchmarking distributed shallow water model on rank $local_rank..."

MPI.Barrier(comm)

trial = @benchmark begin
    @sync_gpu time_step!($model, 1)
end samples=10 evals=1

MPI.Barrier(comm)

t_median = BenchmarkTools.prettytime(median(trial).time)
@info "Done benchmarking on rank $(local_rank). Median time: $t_median"

jldopen("distributed_shallow_water_model_$(R)ranks_$local_rank.jld2", "w") do file
    file["trial"] = trial
end


text:
module DistanceFields

function euclidean_distance_transform!(out::AbstractVector, f::AbstractVector)
    n = length(f)
    k = 1
    v = Array{Int}(undef, n)
    v[1] = 1
    z = Array{Float64}(undef, n + 1)
    z[1] = -Inf
    z[2] = Inf

    for q in 2:n
        @label compute_s
        s = ((f[q] + q ^ 2) - (f[v[k]] + v[k] ^ 2)) / (2 * (q - v[k]))
        if s <= z[k]
            k = k - 1
            @goto compute_s
        else
            k = k + 1
            v[k] = q
            z[k] = s
            z[k + 1] = Inf
        end
    end

    k = 1

    for q in 1:n
        while z[k + 1] < q
            k = k + 1
        end

        out[q] = (q - v[k]) ^ 2 + f[v[k]]
    end

    return nothing
end

end


text:
module TestTMap

using Test
using ThreadingTools
using ThreadingTools: Batches
using BenchmarkTools
using OffsetArrays

@testset "Batches" begin
    for _ in 1:10
        start = rand(-1000:1000)
        stop = start + rand(1:1000)
        inds = start:stop
        batch_size = rand(1:1000)
        b = Batches(inds, batch_size)
        @test vcat([b[i] for i in eachindex(b)]...) == inds
    end

    inds = Base.OneTo(rand(1:1000))
    batch_size = rand(1:1000)
    b = Batches(inds, batch_size)
    @test vcat([b[i] for i in eachindex(b)]...) == inds
end

@testset "test tmap, tmap!" begin

    src = OffsetArray(rand(3), -2:0)
    dst = similar(src)
    @test tmap!(sqrt, dst, src) == Base.map(sqrt, src)
    @test Base.map(sqrt, src) == tmap(sqrt, src, batch_size=2)

    src = (1:2, [1.0, 2.0], [false, true])
    @test Base.map(*, src...) ==  tmap(*, src...)

    for (f, src) in [
              (tuple, tuple([randn(2) for _ in 1:3]...)),
              (sin  , (randn(3),)                      ),
              (+    , (randn(3), randn(3),)            ),
              (+    , (randn(6), randn(2,3),)          ),
              (x->x , (randn(100,100),)                ),
              (atan , (randn(10^4), randn(10^4))       ),
              (*    , (1:2, randn(2), [true, false])   ),
              (tuple, tuple([randn(1) for _ in 1:3]...)),
              (sqrt , (OffsetArray(rand(3), -2:0),)   ),
             ]
        src1 = first(src)
        y = f(map(first, src)...)
        T = typeof(y)
        dst_map = similar(src1, T)
        dst_TTmap = similar(src1, T)

        res_map!  = @inferred map!(f, dst_map, src...)
        res_TTmap! = @inferred tmap!(f, dst_TTmap, src...)
        res_map   = @inferred map(f, src...)
        res_TTmap  = @inferred tmap(f, src...)
        @test typeof(res_map) == typeof(res_TTmap)
        @test res_map == res_TTmap
        @test res_map! == res_TTmap!
    end

    # tmap empty
    @test tmap(+, Int[], Float64[]) == Float64[]

    @test_throws DimensionMismatch tmap(+, randn(5), randn(4))
    @test_throws DimensionMismatch tmap!(+, randn(4), randn(5), randn(4))
    @test_throws DimensionMismatch tmap!(+, randn(5), randn(5), randn(4))

    # allocations
    for n in 1:4
        srcs = [randn(10^4) for _ in 1:n]
        # circumvent bug in julia < 1.2:
        # https://discourse.julialang.org/t/debugging-strange-allocations/25488/5?u=jw3126
        tmap(tuple, srcs...)
        b = @benchmark tmap($tuple, $(srcs...)) evals=1 samples=1
        @test b.allocs < 200
    end
end

end#module


text:
using CUDArt
using CUBLAS
using CudaNNUtils
using Quaternions
using QuaternionArgs
using Base.Test

function test_equiv(n,m)
    for typ in [Float32, Float64]
        A = rand(Quaternion{typ}, n,m)

        d_A = CudaArray(A)
        d_B = CudaArray(Complex{typ}, 2n,2m)

        CudaNNUtils.equiv!(d_B, d_A)
        B = Quaternions.equiv(A)
        h_B = to_host(d_B)
        @test h_B == B

        CudaNNUtils.equiv!(d_A, d_B)
        h_A = to_host(d_A)
        @test h_A == A

        CUDArt.free(d_A)
        CUDArt.free(d_B)
    end
end

function test_equiv_mul(n)
    # parameters
    typ = Float64
    alpha = one(Complex{typ})
    beta = one(Complex{typ})

    # generate quaternion matrices
    A = rand(Quaternion{typ},n,n)
    B = rand(Quaternion{typ},n)
    C = zeros(Quaternion{typ},n)

    # generate equivalent complex matrices
    Ac = Quaternions.equiv(A)
    Bc = Quaternions.equiv(B)
    Cc = Quaternions.equiv(C)

    # allocate device memory
    d_A = CudaArray(Ac)
    d_B = CudaArray(Bc)
    d_C = CudaArray(Cc)

    # calc
    CUBLAS.gemm!('N','N',alpha,d_A,d_B,beta,d_C)
    C = (alpha*A)*B + beta*C

    # compare
    h_C = to_host(d_C)
    @test_approx_eq Quaternions.equiv(C) h_C

    for ptr in [d_A, d_B, d_C]
        CUDArt.free(ptr)
    end
end

function test_fill(n,m)
    for typ in [Float32, Float64]
        A = rand(typ, n,m)

        d_A = CudaArray(A)

        CudaNNUtils.fill!(d_A, zero(typ))
        s = sum(to_host(d_A))
        @test s == zero(typ)

        CUDArt.free(d_A)
    end
end

function test_transpose(n,m)
    for typ in [Quaternion{Float32}, Quaternion{Float64}]
        A = rand(typ, n,m)

        d_A = CudaArray(A)
        d_tA = CudaArray(eltype(A), size(A,2), size(A,1))

        CudaNNUtils.transpose!(d_tA, d_A)
        h_tA = to_host(d_tA)
        @test h_tA == transpose(A)

        CUDArt.free(d_A)
        CUDArt.free(d_tA)
    end
end

function test_ctranspose(n,m)
    for typ in [Quaternion{Float32}, Quaternion{Float64}]
        A = rand(typ, n,m)

        d_A = CudaArray(A)
        d_tA = CudaArray(eltype(A), size(A,2), size(A,1))

        CudaNNUtils.ctranspose!(d_tA, d_A)
        h_tA = to_host(d_tA)
        @test h_tA == A'

        CUDArt.free(d_A)
        CUDArt.free(d_tA)
    end
end

devlist = devices(dev->capability(dev)[1]>=1.3, nmax=1)
if length(devlist) != 1
    warn("NO DEVICE FOR TESTING")
    exit()
end

device_synchronize()
CudaNNUtils.init(devlist)

N = 100
M = 200
test_equiv(N, M)
test_equiv_mul(N)
test_fill(N, M)
test_transpose(N, M)
test_ctranspose(N, M)

CudaNNUtils.close()


text:
function mapper(arr)
    for i = 1:18
        for j = i+1:19
            arr[j] = arr[j]%arr[i] == 0 ? arr[j]/arr[i] : arr[j]
        end
    end
    arr
end

println(reduce(*, mapper(collect(2:20))))

max_stars_repo_path: 005 - Smallest Multiple/005 - Smallest_Multiple.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
using FourierCore.Grid

@testset "Constructing the grid" begin
    grid = FourierGrid(4, Circle())
    k⃗ = grid.wavenumbers[]
    @test k⃗[1] == 0.0
    @test k⃗[2] == 1.0
    @test k⃗[3] == 0.0
    @test k⃗[4] == -1.0
    @test length(k⃗) == 4
    @test grid[1][1] == 0.0

    grid = FourierGrid(4, Circle(1π))
    k⃗ = grid.wavenumbers[]
    @test k⃗[1] == 0.0
    @test k⃗[2] == 2.0
    @test k⃗[3] == 0.0
    @test k⃗[4] == -2.0
    @test length(k⃗) == 4

    grid = FourierGrid((8, 8), Circle()^2)
    (; nodes, wavenumbers) = grid
    @test length(nodes) == 2
    @test length(wavenumbers) == 2
    @test all(wavenumbers[1][:] .== wavenumbers[2][:])
    @test all(LinearAlgebra.size(wavenumbers[1]) .== (8, 1))
    @test all(LinearAlgebra.size(wavenumbers[2]) .== (1, 8))

    grid = FourierGrid((3, 5, 7), Circle()^3)
    (; nodes, wavenumbers) = grid
    @test length(nodes) == 3
    @test length(wavenumbers) == 3
    @test all(LinearAlgebra.size(wavenumbers[1]) .== (3, 1, 1))
    @test all(LinearAlgebra.size(wavenumbers[2]) .== (1, 5, 1))
    @test all(LinearAlgebra.size(wavenumbers[3]) .== (1, 1, 7))

end

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
#***************#
#   CART POLE   #
#***************#
using Test
model, obj0 = Dynamics.cartpole_analytical
n,m = model.n, model.m

Q = copy(obj0.cost.Q)
R = copy(obj0.cost.R)
Qf = copy(obj0.cost.Qf)

x0 = [0;0;0;0.]
xf = [0.5;pi;0;0]
tf = 2.
u_bnd = 50
x_bnd = [0.6,Inf,Inf,Inf]
obj = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd, x_min=-x_bnd, x_max=x_bnd)
obj_con = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd)
cost_mintime = LQRCost(Q*0,R*0.001,Qf,xf)
obj_mintime = update_objective(obj_con,cost=cost_mintime,tf=:min)
dt = 0.05

# Check Jacobians
method = :hermite_simpson
solver = Solver(model,ConstrainedObjective(obj),N=101,integration=:rk3)
N = solver.N
NN = (n+m)N
nG, = get_nG(solver,method)
U0 = ones(1,N)*1
X0 = line_trajectory(obj.x0, obj.xf, N)
Z0 = packZ(X0,U0)

eval_f, eval_g, eval_grad_f, eval_jac_g = gen_usrfun_ipopt(solver,method)
x_L, x_U, g_L, g_U = get_bounds(solver,method)
P = length(g_L)  # Total number of constraints

g = zeros(P)
grad_f = zeros(NN)
rows = zeros(nG)
cols = zeros(nG)
vals = zeros(nG)

reset_evals(solver)
eval_f(Z0)
evals(solver,:f)

reset_evals(solver)
eval_g(Z0,g)
evals(solver,:f)
evals(solver,:c)

reset_evals(solver)
eval_grad_f(Z0,grad_f)
evals(solver,:f)
evals(solver,:Fc)

reset_evals(solver)
eval_jac_g(Z0,:Structure,rows,cols,vals)
evals(solver,:f)
evals(solver,:Fc)
evals(solver,:c)
eval_jac_g(Z0,:vals,rows,cols,vals)
evals(solver,:f)
evals(solver,:Fc)
evals(solver,:c)

grad_f == grad_f2
grad_f2 = copy(grad_f)

solver.opts.verbose = false
sol,stats = solve_dircol(solver,X0,U0, method=method)
@test norm(sol.X[:,end] - obj.xf) < 1e-6
plot(sol.X')
plot(sol.U')

include("../development/dircol_old.jl")
method = :hermite_simpson
function check_grads(solver,method)
    n,m = get_sizes(solver)
    N,N_ = get_N(solver, method)
    U0 = ones(1,N)*1
    X0 = line_trajectory(obj.x0, obj.xf, N)

    res = DircolResults(n,m,solver.N,method)
    res.X .= X0
    res.U .= U0
    Z = res.Z
    Z[1:15] = 1:15

    weights = get_weights(method,N)*dt
    TrajectoryOptimization.update_derivatives!(solver,res,method)
    TrajectoryOptimization.get_traj_points!(solver,res,method)
    TrajectoryOptimization.get_traj_points_derivatives!(solver,res,method)
    TrajectoryOptimization.update_jacobians!(solver,res,method)

    function eval_ceq(Z)
        X,U = TrajectoryOptimization.unpackZ(Z,(n,m,N))
        TrajectoryOptimization.collocation_constraints(X,U,method,solver.dt,solver.model.f)
    end

    function eval_f(Z)
        X,U =TrajectoryOptimization.unpackZ(Z,(n,m,N))
        cost(solver,X,U,weights,method)
    end

    # Check constraints
    g = eval_ceq(Z)
    g_colloc = collocation_constraints(solver, res, method)
    @test g_colloc ≈ g

    Xm = res.X_
    X = res.X

    X1 = zero(Xm)
    X1[:,end] = Xm[:,end]
    for k = N-1:-1:1
        X1[:,k] = 2Xm[:,k] - X1[:,k+1]
    end
    2Xm[:,N-1]-Xm[:,N]

    # Check constraint jacobian
    jacob_g = constraint_jacobian(solver,res,method)
    jacob_g_auto = ForwardDiff.jacobian(eval_ceq,Z)
    @test jacob_g_auto ≈ jacob_g

    # Check cost
    @test cost(solver,res) ≈ eval_f(Z)

    N,N_ = get_N(solver,method)
    # Check cost gradient
    jacob_f_auto = ForwardDiff.gradient(eval_f,Z)
    jacob_f = cost_gradient(solver,res,method)
    @test jacob_f_auto ≈ jacob_f
end

check_grads(solver,:midpoint)
check_grads(solver,:trapezoid)
check_grads(solver,:hermite_simpson)
check_grads(solver,:hermite_simpson_separated)

# Minimum time
method = :hermite_simpsonj
# Set up problem
solver = Solver(model,obj_mintime,N=51,integration=:rk3_foh)
N,N_ = TrajectoryOptimization.get_N(solver,method)
m̄, = get_num_controls(solver)
NN = N*(n+m̄)
nG,Gpart = TrajectoryOptimization.get_nG(solver,method)
# nG = Gpart.collocation
U0 = rand(1,N)*1
if is_min_time(solver)
    U0 = [U0; ones(1,N)*0.1]
end
# X0 = line_trajectory(obj.x0, obj.xf, N)
X0 = rand(n,N)
Z = TrajectoryOptimization.packZ(X0,U0)

# Init vars
g = zeros((N-1)n + N-2)
grad_f = zeros(NN)
rows = zeros(nG)
cols = zeros(nG)
vals = zeros(nG)
weights = get_weights(method,N_)

# Get functions and evaluate
eval_f, eval_g, eval_grad_f, eval_jac_g = gen_usrfun_ipopt(solver,method)

# Test Cost
eval_f(Z)
J = cost(solver::Solver,X0,U0,weights::Vector,method::Symbol)
@test eval_f(Z) ≈ J

# Cost Gradient
eval_grad_f(Z,grad_f)
grad_f
grad_f_check = ForwardDiff.gradient(eval_f,Z)
@test grad_f_check ≈ grad_f

# Constraint Jacobian
eval_g(Z,g)
eval_jac_g(Z, :Structure, rows, cols, vals)
eval_jac_g(Z, :vals, rows, cols, vals)
function eval_g_wrapper(Z)
    gin = zeros(eltype(Z),size(g))
    eval_g(Z,gin)
    return gin
end
eval_g_wrapper(Z)
jac_g_check = ForwardDiff.jacobian(eval_g_wrapper,Z)
jac_g = Array(sparse(rows,cols,vals))
@test jac_g ≈ jac_g_check

# Solver integration scheme should set the dircol scheme
U0 = ones(1,N)
solver = Solver(model,obj,dt=dt,integration=:midpoint)
sol,stats = solve_dircol(solver,X0,U0)
@test norm(sol.X[:,end]-obj.xf) < 1e-5
@test stats["info"] == :Solve_Succeeded

solver = Solver(model,obj,dt=dt)
sol2,stats = solve_dircol(solver,X0,U0,method=:midpoint)
@test norm(sol2.X[:,end]-obj.xf) < 1e-5
@test sol.X == sol2.X
@test sol.U == sol2.U
@test stats["info"] == :Solve_Succeeded

solver = Solver(model,obj,dt=dt,integration=:rk3_foh)
sol,stats = solve_dircol(solver,X0,U0)
@test norm(sol.X[:,end]-obj.xf) < 1e-5
@test stats["info"] == :Solve_Succeeded

solver = Solver(model,obj,dt=dt)
sol2,stats = solve_dircol(solver,X0,U0,method=:hermite_simpson)
@test norm(sol2.X[:,end]-obj.xf) < 1e-5
@test sol.X == sol2.X
@test sol.U == sol2.U
@test stats["info"] == :Solve_Succeeded

# No initial guess
mesh = [0.2]
t1 = @elapsed sol,stats = solve_dircol(solver)
@test norm(sol.X[:,end]-obj.xf) < 1e-5

# Test dircol constraint stuff
n,m = 3,2
cE!(c,x,u) = begin c[1:2] = 2x[1:2]+u; c[3] = x'x + 5 end
pE = 3
cE!(c,x) = begin c[1] = cos(x[1]) + x[2]*x[3]; c[2] = x[1]*x[2]^2 end
pE_N = 2
cI!(c,x,u) = begin c[1] = x[3]-x[2]; c[2] = u[1]*x[1] end
pI = 2
pI_N = 0

x = rand(n); u = rand(m);
c = zeros(3)
cI!(c,x,u)
is_inplace_function(cE!,x,u)
model, obj = Dynamics.dubinscar
obj_con = ConstrainedObjective(obj,cE=cE!,cI=cI!,cE_N=cE!,use_xf_equality_constraint=false,tf=3)
@test obj_con.p == pE + pI
@test obj_con.pI == pI
@test obj_con.pI_N == pI_N
@test obj_con.p_N == pE_N + pI_N

method = :trapezoid
solver = Solver(model,obj_con,dt=1.)
N, = TrajectoryOptimization.get_N(solver,method)
NN = (n+m)*N
X = [1. 2 3 4; 1 2 3 4; 1 2 3 4]
U = [0. 1 0 0; -1 0 -1 0]
# U = ones(m,N)
# X = line_trajectory(obj.x0,obj.xf,N)
x,u = X[:,1],U[:,1]
c = zeros(3)
@test is_inplace_function(cE!,x,u)
@test 3 == count_inplace_output(cE!,x,u)

# Constraint function
pI_obj, pE_obj = TrajectoryOptimization.count_constraints(obj_con)
@test pI_obj == (pI, pI, 0, 0)
@test pE_obj == (pE, pE, pE_N, pE_N)
p_total = (pE + pI)*(N-1) + pE_N + pI_N
p_colloc = (N-1)n

c_fun!, jac_c = TrajectoryOptimization.gen_custom_constraint_fun(solver, method)
C = zeros(p_total)
c_fun!(C,X,U)
Z = packZ(X,U)
cE(x,u) = begin c=zeros(3); cE!(c,x,u); return c end
cE(x) = begin c=zeros(2); cE!(c,x); return c end
cI(x,u) = begin c=zeros(2); cI!(c,x,u); return c end
C_expected = [cE(X[:,1],U[:,1]);
              cE(X[:,2],U[:,2]);
              cE(X[:,3],U[:,3]);
              cE(X[:,4]);
              cI(X[:,1],U[:,1]);
              cI(X[:,2],U[:,2]);
              cI(X[:,3],U[:,3])]
@test C_expected == C

x_L,x_U, g_L,g_U = get_bounds(solver,method)

# Jacobian
c_fun!, jac_c!, jac_sparsity = TrajectoryOptimization.gen_custom_constraint_fun(solver, method)
nG,Gpart = get_nG(solver,method)
nP = Gpart.custom

J_struct = jac_sparsity()
rows,cols,inds = findnz(J_struct)
v = sortperm(inds)
rows = rows[v]
cols = cols[v]

J = zeros(nP)
jac_c!(J,X,U)
J1 = Array(sparse(rows,cols,J,p_total,NN))
Z = packZ(X,U)

function c_funZ!(C,Z)
    X,U = unpackZ(Z,(n,m,N))
    c_fun!(C,X,U)
end
C2 = zeros(p_total)
c_funZ!(C2,Z)

J2 = zeros(p_total,NN)
jac_cZ!(J,Z) = ForwardDiff.jacobian!(J,c_funZ!,C2,Z)
jac_cZ!(J2,Z)
@test J2 == J1

vals = zeros(nG)
jac_c(vals,X,U)
vals

eval_f, eval_g, eval_grad_f, eval_jac_g = gen_usrfun_ipopt(solver,method)
g = zeros(p_total + p_colloc)
eval_g(Z,g)
@test g[p_colloc+1:end] == C

collocation_constraint_jacobian_sparsity(solver,method)
vals = zeros(nG)
rows = zeros(nG)
cols = zeros(nG)
eval_jac_g(Z,:Structure,rows,cols,vals)
eval_jac_g(Z,:vals,rows,cols,vals)
J_all = Array(sparse(rows,cols,vals))
J3 = J_all[p_colloc.+1:end,:]
@test J3 == J1

function eval_g_wrapper(Z)
    gin = zeros(eltype(Z),size(g))
    eval_g(Z,gin)
    return gin
end

auto_diff(Z) = ForwardDiff.jacobian(eval_g_wrapper,Z)
jac_g_check = auto_diff(Z)
@test jac_g_check == sparse(rows,cols,vals)

using BenchmarkTools
using Juno
using Profile
Profile.init(delay=1e-6)
Profile.clear()
@profile eval_jac_g(Z,:vals,rows,cols,vals)
Juno.profiler()
@btime eval_jac_g(Z,:vals,rows,cols,vals)
@btime auto_diff(Z)

@time jac_c(J,X,U)
@time jac_cZ!(J2,Z)

# DircolResults
n,m,N = (4,2,51)
res = DircolResults(n,m,N,:midpoint)
res.vars
res.Z[1] = 10.
@test res.X[1] == 10
@test res.U[1] === res.Z[5]
@test res.U[1] === res.U_[1]
@test size(res.U_) == (m,N)
@test res.vars.X == res.X
@test res.vars.Z == res.Z

res = DircolResults(n,m,N,:hermite_simpson_separated)
@test size(res.X) == (n,2N-1)
@test size(res.X_) == (n,2N-1)
@test res.X[1] === res.Z[1]
@test size(res.U) == (m,2N-1)
@test res.U === res.U_
@test res.X === res.X_

res = DircolResults(n,m,N,:hermite_simpson)
@test size(res.X) == (n,N)
@test size(res.X_) == (n,2N-1)
@test res.X[1] === res.Z[1]

res = DircolResults(n,m,N,:trapezoid)
@test res.X === res.X_

# Dircol Vars
X0 = rand(n,N)
U0 = rand(m,N)
Z0 = packZ(X0,U0)
vars = DircolVars(Z0,n,m,N)
@test X0 == vars.X
@test U0 == vars.U
@test Z0 === vars.Z
vars = DircolVars(X0,U0)
@test X0 == vars.X
@test U0 == vars.U
@test Z0 == vars.Z


text:
@testset "Solves" begin
# Test max iterations
solver = ALTROSolver(Problems.Pendulum()...)
set_options!(solver,
    iterations=100,
    cost_tolerance_intermediate=1e-10,
    gradient_tolerance_intermediate=1e-10,
    constraint_tolerance=1e-10,
    show_summary=false
)
solve!(solver)
status(solver)
@test iterations(solver) == solver.opts.iterations
@test status(solver) == Altro.MAX_ITERATIONS

## test max cost
model = RobotZoo.Quadrotor()
N = 101
tf = 5.0
x0,u0 = zeros(model)
xf = MArray(x0); xf[1] = 100
Q = Diagonal(@SVector ones(13))*1e6
R = Diagonal(@SVector ones(4))
obj = LQRObjective(Q,R,Q*N,xf,N,uf=u0)
prob = Problem(model, obj, xf, tf, x0=x0)
U0 = [u0*1e1 for k = 1:N-1]
initial_controls!(prob, U0)

solver = ALTROSolver(prob, show_summary=false)
Altro.is_constrained(solver)
solve!(solver)
@test iterations(solver) == 1
@test status(solver) == Altro.MAXIMUM_COST

## test max state value
model = RobotZoo.Quadrotor()
N = 101
tf = 5.0
x0,u0 = zeros(model)
xf = MArray(x0); xf[1] = -1e3
Q = Diagonal(SA[10,10,10, 1,1,1,1, 10,10,10, 10,10,10.])*1e-1
R = Diagonal(@SVector ones(4))
obj = LQRObjective(Q,R,Q*N,xf,N,uf=u0)
prob = Problem(model, obj, xf, tf, x0=x0)
U0 = [u0 for k = 1:N-1]

solver = ALTROSolver(prob, show_summary=false)
set_options!(solver, max_state_value=100)
@test Altro.is_constrained(solver) == false
solve!(solver)
@test status(solver) == Altro.STATE_LIMIT

end


text:
using FileIO
using DepthRenderer
using Geometry: CameraIntrinsics, I4, TriangleMesh, SceneGraphNode
using Test
using Printf: @sprintf

width = 20
height = 20
near = 0.001
far = 5.0
cam = CameraIntrinsics(
    width=width, height=height,
    fx=width, fy=height, cx=div(width, 2), cy=div(height, 2), near=near, far=far)
renderer = Renderer(cam)

# triangle 1
a = Float32[-0.25, -0.25, -1]
b = Float32[0.25, -0.25, -1]
c = Float32[0.0, 0.25, -1]
vertices = hcat(a, b, c)
indices = hcat(UInt32[0, 1, 2])
triangle1 = TriangleMesh(vertices, indices)
register_mesh!(renderer, triangle1)

# triangle 2
a = Float32[-1, -1, -2]
b = Float32[1, -1, -2]
c = Float32[0, 1, -2]
vertices = hcat(a, b, c)
indices = hcat(UInt32[0, 1, 2])
triangle2 = TriangleMesh(vertices, indices)
register_mesh!(renderer, triangle2)

# triangle 3
a = Float32[1.5, 2, -4]
b = Float32[2, 2, -4]
c = Float32[2, 1.5, -4]
vertices = hcat(a, b, c)
indices = hcat(UInt32[0, 1, 2])
triangle3 = TriangleMesh(vertices, indices)
register_mesh!(renderer, triangle3)

# scene graph
root = SceneGraphNode(
    children=[
        SceneGraphNode(mesh=triangle1),
        SceneGraphNode(mesh=triangle2),
        SceneGraphNode(mesh=triangle3)])

model = I4(Float32)
view = I4(Float32)

draw!(renderer, root, model, view)

depth_image = get_depth_image!(renderer; show_in_window=true)
@test size(depth_image) == (width, height)

destroy!(renderer)

# print it out
using Printf: @sprintf
for col=height:-1:1
    for row=1:width
        str = @sprintf("%d ", round(depth_image[row, col]))
        print(str)
    end
    println()
end

#
#  depth_image
#
#     row=1                              row=width
#
#  ^
#  |   5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 4   col=height
#  |   5 5 5 5 5 5 5 5 5 2 2 5 5 5 5 5 5 5 5 4
#  |   5 5 5 5 5 5 5 5 5 2 2 5 5 5 5 5 5 5 5 5
#  |   5 5 5 5 5 5 5 5 2 2 2 2 5 5 5 5 5 5 5 5
#  |   5 5 5 5 5 5 5 5 2 2 2 2 5 5 5 5 5 5 5 5
#  |   5 5 5 5 5 5 5 2 2 2 2 2 2 5 5 5 5 5 5 5
#  |   5 5 5 5 5 5 5 2 2 1 1 2 2 5 5 5 5 5 5 5
#  |   5 5 5 5 5 5 2 2 2 1 1 2 2 2 5 5 5 5 5 5
#  Y   5 5 5 5 5 5 2 2 1 1 1 1 2 2 5 5 5 5 5 5
#  |   5 5 5 5 5 2 2 2 1 1 1 1 2 2 2 5 5 5 5 5
#  |   5 5 5 5 5 2 2 1 1 1 1 1 1 2 2 5 5 5 5 5
#  |   5 5 5 5 2 2 2 1 1 1 1 1 1 2 2 2 5 5 5 5
#  |   5 5 5 5 2 2 1 1 1 1 1 1 1 1 2 2 5 5 5 5
#  |   5 5 5 2 2 2 1 1 1 1 1 1 1 1 2 2 2 5 5 5
#  |   5 5 5 2 2 1 1 1 1 1 1 1 1 1 1 2 2 5 5 5
#  |   5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 5
#  |   5 5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5 5
#  |   5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5
#  |   5 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5
#  |   2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2   col=1
#  |
#  *----------------X------------------------->
#
# NOTE Z-axis is point towards camera (objects in front of camera have negative Z-coordinate)

# TODO figure out why it's not exact, why there is < 1e-3 error
@test isapprox(depth_image[1, height], far, atol=1e-3) # background (upper left corner)
@test isapprox(depth_image[1, 1], 2.0, atol=1e-3) # triangle 2 (lower left corner)
@test isapprox(depth_image[width, 1], 2.0, atol=1e-3) # triangle 2 (lower right corner)
@test isapprox(depth_image[div(width, 2), height-1], 2.0, atol=1e-3) # triangle 2 (top middle)
@test isapprox(depth_image[div(width, 2), div(height, 2)], 1.0, atol=1e-3) # triangle 1 (center)
@test isapprox(depth_image[width, height], 4.0, atol=1e-3) # triangle 3 (upper right corner)

###################################
# test round trip image save/load #
###################################

depth_meas = rand(UInt16, 10, 5)
fname = "$(tempname()).png"
save_depth_image(depth_meas, fname)
loaded_meas = load_depth_image(fname)
rm(fname)
@test depth_meas == loaded_meas

max_stars_repo_name: probcomp/DepthRenderer

max_issues_repo_name: probcomp/DepthRenderer

max_forks_repo_name: probcomp/DepthRenderer


text:
using Test
using LinearAlgebra
using SummationByPartsOperators

D_test_list = (MattssonNordström2004(), MattssonSvärdNordström2004(),
                MattssonSvärdShoeybi2008(), Mattsson2014(),
                MattssonAlmquistCarpenter2014Extended())
Di_test_list = (MattssonSvärdNordström2004(),)

# Nonperiodic: Test symmetry and eigenvalues
for source_D in D_test_list, source_Di in Di_test_list, acc_order in 2:2:8, T in (Float32,Float64)
    xmin = -one(T)
    xmax = 2*one(T)
    N = 101
    der_order = 1

    D = try
        derivative_operator(source_D, der_order, acc_order, xmin, xmax, N)
    catch err
        !isa(err, ArgumentError) && throw(err)
        nothing
    end
    D == nothing && continue

    @inferred mass_matrix(D)
    H = mass_matrix(D)

    for order in 2:2:8
        Di = try
            dissipation_operator(source_Di, D, order=order)
        catch err
            !isa(err, ArgumentError) && throw(err)
            nothing
        end
        Di == nothing && continue

        println(devnull, Di)
        println(devnull, Di.coefficients)
        HDi = H*Matrix(Di)
        @test norm(HDi - HDi') < 10*eps(T)
        @test minimum(real, eigvals(HDi)) < 10*eps(T)
    end
end

# Nonperiodic: Compare mul! with β=0 and mul! without β.
for T in (Float32, Float64), order in (2,4,6,8)
    xmin = zero(T)
    xmax = 5*one(T)
    N = 51
    D = derivative_operator(MattssonSvärdNordström2004(), 1, order, xmin, xmax, N)
    x = grid(D)
    u = x.^5
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    Di_serial = dissipation_operator(D, order=order, mode=FastMode())
    Di_threads= dissipation_operator(D, order=order, mode=ThreadedMode())
    Di_full   = Matrix(Di_serial)

    mul!(dest1, Di_serial, u, one(T), zero(T))
    mul!(dest2, Di_serial, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest1, Di_threads, u, one(T), zero(T))
    mul!(dest2, Di_threads, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    dest3 = Di_serial*u
    @test all(i->dest1[i] ≈ dest3[i], eachindex(u))
    dest3 = Di_full*u
    @test all(i->isapprox(dest1[i], dest3[i], atol=800000*eps(T)), eachindex(u))
end

# Periodic: Test symmetry and eigenvalues
for acc_order in 2:2:8, T in (Float32,Float64)
    xmin = -one(T)
    xmax = 2*one(T)
    N = 100
    der_order = 1

    D = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N)
    for order in 2:2:8
        Di = dissipation_operator(D, order=order)
        println(devnull, Di)
        HDi = Matrix(Di)
        @test norm(HDi - HDi') < 10*eps(T)
        @test minimum(real, eigvals(HDi)) < 10*eps(T)
    end
end

# Periodic: Compare mul! with β=0 and mul! without β.
for T in (Float32, Float64), order in (2,4,6,8)
    xmin = zero(T)
    xmax = 5*one(T)
    N = 50
    D = periodic_derivative_operator(1, order, xmin, xmax, N)
    x = grid(D)
    u = x.^5
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))

    Di_serial = dissipation_operator(D, order=order, mode=FastMode())
    Di_threads= dissipation_operator(D, order=order, mode=ThreadedMode())
    Di_safe   = dissipation_operator(D, order=order, mode=SafeMode())
    Di_full   = Matrix(Di_serial)

    mul!(dest1, Di_serial, u, one(T), zero(T))
    mul!(dest2, Di_serial, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    mul!(dest1, Di_threads, u, one(T), zero(T))
    mul!(dest2, Di_threads, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    dest3 = Di_serial*u
    @test maximum(abs, dest1 - dest3) < u[end]*N*eps(T)
    dest3 = Di_full*u
    @test maximum(abs, dest1 - dest3) < 5*u[end]*N*eps(T)
    dest3 = Di_safe*u
    @test maximum(abs, dest1 - dest3) < 5*u[end]*N*eps(T)
end

# Compare periodic and nonperiodic.
for T in (Float32, Float64), order in (2,4,6,8)
    xmin = zero(T)
    xmax = 5*one(T)
    N = 50
    D = derivative_operator(MattssonSvärdNordström2004(), 1, order, xmin, xmax, N+1)
    Dp = periodic_derivative_operator(1, order, xmin, xmax, N)

    Di = dissipation_operator(D)
    Dip = dissipation_operator(Dp)

    @test norm(Matrix(Di)[15:35,15:35] - Matrix(Dip)[15:35,15:35]) < eps(T)
end


text:
#Reineke SDI calculation
function sdi(s::Stand)
    return s.trees_per_area * ((s.qmd / 10.0) ^ 1.605)
end

function qmd(s::Stand)
    sum_d2 = 0.0
    @simd for i in eachindex(s.treelist)
    @inbounds sum_d2 += (s.treelist[i].diameter ^ 2 * s.treelist[i].expansion)
    end
    qmd = ismissing(s.trees_per_area) ? sqrt(sum_d2 / s.trees_per_area) : sqrt(sum_d2 / trees_per_area(s))
    return qmd
end

function qmd!(s::Stand)
    s.qmd = qmd(s)
end

function basal_area(s::Stand)
    ba = 0.0
    @simd for i in eachindex(s.treelist)
    @inbounds ba += basal_area(s.treelist[i]) * s.treelist[i].expansion
    end
    return ba
end

function basal_area!(s::Stand)
    s.basal_area = basal_area(s)
end

function basal_area!(s::Stand,m::Metric)
    s.basal_area = basal_area(s, metric)
end

function basal_area(t::Tree)
    return ((t.diameter ^ 2) * K)
end

function basal_area(t::Tree,m::Metric)
    return ((t.diameter ^ 2) * KMETRIC)
end

function trees_per_area(s::Stand)
    tpa = 0.0
    @simd for i in eachindex(s.treelist)
    @inbounds tpa += s.treelist[i].expansion
    end
    return tpa
end

function trees_per_area!(s::Stand)
    s.trees_per_area = trees_per_area(s)
end


text:
# using TFIMSampler

using StatsBase
using LinearAlgebra
# using BitBasis
# using ProgressMeter

# using Test
function ED(h)
    Z = [0 1.0 ; 1 0]
    X = [1 0.0 ; 0 -1]
    J = 1
    H_XX =  kron(X, X, I(2), I(2)) + kron(I(2), X, X, I(2)) + kron(I(2), I(2), X, X) # + kron(X, I(2), I(2), X)
    H_Z = kron(Z, I(2), I(2), I(2)) + kron(I(2), Z, I(2), I(2)) + kron(I(2), I(2), Z, I(2)) + kron(I(2), I(2), I(2), Z)
    H = -J * H_XX - h * H_Z |> Matrix
    e, psis = eigen(H)
    e = e[1]
    gs_exact = psis[:,1]

    EJ = gs_exact' * J * H_XX * gs_exact# / norm(gs_exact)
    Eh = gs_exact' * h * H_Z * gs_exact # / norm(gs_exact)

    M = 0.0
    ##### Calculate the groundstate magnetization <m^2> in the Z direction
    magnetization = 0
    Dim = 2^N
    for Ket = 0:Dim - 1  # Loop over Hilbert Space
        SumSz = 0.
        for SpinIndex = 0:N - 1  # Loop over spin index (base zero, stop one spin before the end of the chain)
            Spin1 = 2 * ((Ket >> SpinIndex) & 1) - 1
            SumSz += Spin1 # spin is +1 or -1
            # print(Spin1," ")
        end
        # println(SumSz," ",GroundState[Ket+1])
        magnetization += SumSz * SumSz * gs_exact[Ket + 1]^2  # Don't forgot to square the coefficients...
    end
    M = magnetization / (N * N)
    return Eh + EJ, M
end
    # display(gs_exact)
# function test_ed(h)
#     L = 4
#     @info "h = $h"
#     H_XX =  kron(X, X, I(2), I(2)) + kron(I(2), X, X, I(2)) + kron(I(2), I(2), X, X) + kron(X, I(2), I(2), X)
#     H_Z = kron(Z, I(2), I(2), I(2)) + kron(I(2), Z, I(2), I(2)) + kron(I(2), I(2), Z, I(2)) + kron(I(2), I(2), I(2), Z)
#     H = -H_XX + h * H_Z |> Matrix
#     e, psis = eigen(H)
#     e = e[1]
#     gs_exact = normalize(psis[:,1].^2, 1)

#     M = TFIMSampler.pairing(L, h)
#     unnormalizedamps, amps = TFIMSampler.get_all_probs(M)

#     bits = [bitarray(i, L) for i in 0:2^L - 1]
#     display(gs_exact)
#     display(amps)
#     # display(hcat(bits, amps, gs_exact))
#     @test isapprox(gs_exact, amps; atol = 1e-5)
# end

# # @testset "ED test" begin
# #     test_ed(1)
# #     test_ed(2)
# #     test_ed(0.21)
# # end

# # configs = TFIMSampler.sample()

# function test_mc(L, h)
#     N = 10^4
#     bits = [bitarray(i, L) for i in 0:2^L - 1]
#     configs = TFIMSampler.sample(L = L, h = h, N = N)# |> sort
#     amps = countmap(configs)
#     amps = filter(x->40506903398877536 > x > -40506031398877536, amps.vals) |> sort
#     amps = amps / (2L)

#     M = TFIMSampler.pairing(L, h)
#     amps_exact = TFIMSampler.get_all_probs(M)[2] |> sort
#     amps_exact *= N
#     amps_exact = filter(x->x > 2, amps_exact) |> sort
#     amps_exact = round.(amps_exact, digits = 1)

#     Z = sum(amps)
#     Z2 = sum(amps_exact)
#     @test isapprox(amps/Z, amps_exact/Z2; atol = 1e-2)
# end

# @testset "MCMC test" begin
#     test_mc(4, 1.0)
#     test_mc(4, 2.0)
#     test_mc(4, 0.21)
#     test_mc(16, 1.0)
# end

# test_mc(32, 1.0)


text:
@cdef :PySequence_Length Py_ssize_t (PyPtr,)
@cdef :PySequence_GetItem PyPtr (PyPtr, Py_ssize_t)
@cdef :PySequence_SetItem Cint (PyPtr, Py_ssize_t, PyPtr)
@cdef :PySequence_Contains Cint (PyPtr, PyPtr)


text:
using AVLTrees, BenchmarkTools
using Random
using Plots
using DataFrames

function batch_insert!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        insert!(t, i, i)
    end
end

function batch_delete!(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        delete!(t, i)
    end
end

function batch_find(t::AVLTree{K,D}, v::Vector{K}) where {K,D}
    for i in v
        i in t
    end
end

insertion_vec = []
deletion_vec = []
search_vec = []

d = DataFrame((op=[], time=[], n=[]))
x = [1_000, 10_000, 100_000, 1_000_000, 10_000_000]

function prepare_t(t)
    _t = deepcopy(t)
    for i in nums_test
        insert!(_t, i, i)
    end
    _t
end

for attempt in 1:1
    for N in x
        global t = AVLTree{Int64,Int64}()
        rng = MersenneTwister(1111)
        global nums_fill = rand(rng, Int64, N)
        global nums_test = rand(rng, Int64, 10_000)

        for i in nums_fill
            insert!(t, i, i)
        end

        insertion = @benchmark batch_insert!(_t, nums_test) setup=(_t =deepcopy(t))
        search = @benchmark batch_find(t, nums_test) setup=(_t = prepare_t(t))
        deletion = @benchmark batch_delete!(t, nums_test) setup=(_t = prepare_t(t))

        push!(d, ("insert", minimum(insertion).time, N))
        push!(d, ("delete", minimum(deletion).time,N))
        push!(d, ("search", minimum(search).time,N))
        println("done $N")
    end
end

c = combine(groupby(d, [:op,:n]), :time => minimum)

# plot(x, insertion_vec/1000, xscale=:log10, ylabel="us")
# plot(x, deletion_vec/1000, xscale=:log10, ylabel="us")
# plot(x, search_vec/1000, xscale=:log10, ylabel="us")

plot(
    x,
    [c[(c.op.=="insert"),:].time_minimum,c[(c.op.=="delete"),:].time_minimum, c[(c.op.=="search"),:].time_minimum],
    xscale = :log10,
    ylabel = "operation time [us]",
    xlabel = "N",
    xticks = [1e3, 1e4, 1e5, 1e6, 1e7],
    markershape =[:diamond :utriangle :dtriangle],
    labels= ["insert" "delete" "lookup"],
    legend=:topleft,
)

savefig("branch_results_new2.svg")
savefig("result_new2.png")
using CSV
CSV.write("branch_results_new2.csv", c)


text:
# This file contains functions/methods related to images/image generation.

"""
    contraststretch!(image::Matrix{Float64};
                     minval::Float64 = 0.0, maxval::Float64 = 1.0)

Stretch the histogram of `image` to fill the range `minval` to `maxval`.

# Inputs
-`image`: Image to be contrast stretched.
-`minval`: Minimum value of pixels in `image` once stretched.
-`maxval`: Maximum value of pixels in `image` once stretched.
"""
function contraststretch!(image::Matrix{Float64};
    minval::Float64 = 0.0, maxval::Float64 = 1.0)
    # Stretch the image to the range [0.0, 1.0].
    image .-= StatsBase.minimum(image)
    image ./= StatsBase.maximum(image)
    image = (maxval - minval) * image .+ minval

    return image
end

"""
    outimage = contraststretch(image::Matrix{Float64};
                            minval::Float64 = 0.0, maxval::Float64 = 1.0)

Stretch the histogram of `image` to fill the range `minval` to `maxval`.

# Inputs
-`image`: Image to be contrast stretched.
-`minval`: Minimum value of pixels in `image` once stretched.
-`maxval`: Maximum value of pixels in `image` once stretched.

# Outputs
-`outimage`: Contrast stretched copy of input `image`.
"""
function contraststretch(image::Matrix{Float64};
    minval::Float64 = 0.0, maxval::Float64 = 1.0)
    return SMLMData.contraststretch!(deepcopy(image);
        minval = minval, maxval = maxval)
end

"""
    image = makecircleim(coords::Matrix{Float64},
                         σ::Vector{Float64},
                         datasize::Vector{Int},
                         mag::Float64 = 20.0)

Make a circle image of the localizations in `coords`.

# Description
This function creates an image of the localizations in `coords` by adding a
circle centered at the locations `coords` with radii `σ`.

# Inputs
-`coords`: Localization coordinates. ([y x])
-`σ`: Standard error of localizations in `coords`. (nlocx1)
-`datasize`: Size of the data image. ([ysize xsize])
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} histogram image of localizations.
"""
function makecircleim(coords::Matrix{Float64},
    σ::Vector{Float64},
    datasize::Vector{Int},
    mag::Float64 = 20.0)
    # Rescale the coordinates based on `mag`.
    coords = mag * (coords .- 0.5) .+ 0.5
    σ *= mag

    # Loop through localizations and add them to our output image.
    imagesize = Int.(round.(datasize * mag))
    image = zeros(Float64, imagesize[1], imagesize[2])
    for nn = 1:size(coords, 1)
        # If σ[nn] isn't positive, skip this localization.
        if !(σ[nn] > 0.0)
            continue
        end

        # Define the pixel locations that fall along the circle.
        # NOTE: The extra factor of 4 improves circle appearance.
        θ = range(0, 2 * pi, length = max(4, Int(ceil(4 * (2 * pi * σ[nn])))))
        rows = Int.(round.(coords[nn, 1] .+ σ[nn] * sin.(θ)))
        cols = Int.(round.(coords[nn, 2] .+ σ[nn] * cos.(θ)))
        validind = findall((rows .>= 1) .* (rows .< imagesize[1]) .*
                           (cols .>= 1) .* (cols .< imagesize[1]))

        # Set the pixels of the output image to 1.0 wherever met by the circle.
        for ii in validind
            image[rows[ii], cols[ii]] = 1.0
        end
    end

    return image
end

"""
    image = makecircleim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)

Make a circle image of the localizations in `smld`.

# Description
This function creates an image of the localizations in `smld` by adding a
circle for each localization.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} circle image of localizations.
"""
function makecircleim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)
    coords = [smld.y smld.x]
    σ = vec(mean([smld.σ_y smld.σ_x], dims = 2))
    return SMLMData.makecircleim(coords, σ, smld.datasize, mag)
end

"""
    image = circleim(smld::SMLMData.SMLD2D, pxsize::Float64;
                     pxsize_out::Float64 = 0.005)

Generate and save a circle image.

# Description
This method is a wrapper for makecircleim() with modified inputs.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`pxsize`: Pixel size of localizations in `smld`. (micrometers)
-`pxsize_out`: Desired output pixel size. (micrometers)
"""
function circleim(smld::SMLMData.SMLD2D, pxsize::Float64;
    pxsize_out::Float64 = 0.005)
    # Determine the requested magnfication factor.
    mag = pxsize / pxsize_out

    # Generate the circle image.
    image = SMLMData.makecircleim(smld, mag)

    # Perform a contrast stretch so that the pixels fills the range [0, 1].
    SMLMData.contraststretch!(image)

    return image
end

"""
    image = circleim(smld::SMLMData.SMLD2D, pxsize::Float64, filename::String;
                     pxsize_out::Float64 = 0.005)

Generate and save a circle image.

# Description
This method is a wrapper for makecircleim() which uses different inputs
arguments and saves the image.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`pxsize`: Pixel size of localizations in `smld`. (micrometers)
-`pxsize_out`: Desired output pixel size. (micrometers)
"""
function circleim(smld::SMLMData.SMLD2D, pxsize::Float64, filename::String;
    pxsize_out::Float64 = 0.005)
    # Generate the circle image.
    image = SMLMData.circleim(smld, pxsize; pxsize_out = pxsize_out)

    # Save the image.
    save(filename, Images.colorview(Gray, image))

    return image
end

"""
    image = makegaussim(μ::Matrix{Float64},
                        σ_μ::Matrix{Float64},
                        datasize::Vector{Int};
                        mag::Float64 = 20.0,
                        nsigma::Float64 = 5.0)

Make a Gaussian image of the localizations in `μ`.

# Description
This function creates an image of the localizations defined by `μ` and `σ_μ`
in which Gaussians with standard deviations `σ_μ` are added at positions `μ`.

# Inputs
-`μ`: Positions of the Gaussians. (pixels)([y x])
-`σ_μ`: Standard errors of `μ` estimates. (pixels)([y x])
-`datasize`: Size of the region of data collection. (pixels)([y; x])
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)
-`nsigma`: Number of standard deviations from the localization coordinate at
           which we truncate the Gaussian. (Default = 5.0)

# Outputs
-`image`: Matrix{Float64} Gaussian image in which each localization in `smld`
          is plotted as a Gaussian.
"""
function makegaussim(μ::Matrix{Float64},
    σ_μ::Matrix{Float64},
    datasize::Vector{Int};
    mag::Float64 = 20.0,
    nsigma::Float64 = 5.0)
    # Loop through emitters and add them to our output Gaussian image.
    imagesize = Int.(round.(datasize * mag))
    image = zeros(Float64, imagesize[1], imagesize[2])
    for nn = 1:size(μ, 1)
        # Prepare a normal distribution for this emitter.
        if !all(σ_μ[nn, :] .> 0.0)
            continue
        end
        Σ = [σ_μ[nn, 1]^2 0.0; 0.0 σ_μ[nn, 2]^2]
        distrib = Distributions.MvNormal(μ[nn, :], Σ)

        # Loop through pixels of the image and add this emitter.
        ystart = max(1,
            Int(round(mag * (μ[nn, 1] - nsigma * σ_μ[nn, 2] - 0.5))))
        yend = min(imagesize[1],
            Int(round(mag * (μ[nn, 1] + nsigma * σ_μ[nn, 2]))))
        xstart = max(1,
            Int(round(mag * (μ[nn, 2] - nsigma * σ_μ[nn, 1] - 0.5))))
        xend = min(imagesize[2],
            Int(round(mag * (μ[nn, 2] + nsigma * σ_μ[nn, 1]))))
        for ii = ystart:yend, jj = xstart:xend
            image[ii, jj] +=
                Distributions.pdf(distrib, ([ii; jj] .- 0.5) / mag .+ 0.5)
        end
    end

    # Normalize the image to sum to 1.0.  If any image values are NaN, set them
    # to 0.0.
    nanpixels = isnan.(image)
    if any(nanpixels)
        image[nanpixles] .= 0.0
    end
    image = image ./ sum(image)
    if !isapprox(sum(image), 1.0)
        @warn "Image is non-normalizable!  Returning flat image."
        image = ones(Float64, imagesize[1], imagesize[2]) ./ prod(imagesize)
    end

    return image
end

"""
    image = makegaussim(smld::SMLMData.SMLD2D;
                        mag::Float64 = 20.0,
                        nsigma::Float64 = 5.0)

Make a Gaussian image of the localizations in `smld`.

# Description
This function creates an image of the localizations in `smld` by placing a
Gaussian truncated to `nsigma` at the localization coordinates, where the
standard deviation is given by `smld.σ_x` and `smld.σ_y`.  The image is then
normalized such that it sums to 1.0.  The background signal is not accounted
for in this method.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)
-`nsigma`: Number of standard deviations from the localization coordinate at
           which we truncate the Gaussian. (Default = 5.0)

# Outputs
-`image`: Matrix{Float64} Gaussian image in which each localization in `smld`
          is plotted as a Gaussian.
"""
function makegaussim(smld::SMLMData.SMLD2D;
    mag::Float64 = 20.0,
    nsigma::Float64 = 5.0)
    return SMLMData.makegaussim([smld.y smld.x], [smld.σ_y smld.σ_x],
        smld.datasize; mag = mag, nsigma = nsigma)
end

"""
    image = gaussim(smld::SMLMData.SMLD2D, pxsize::Float64;
                    pxsize_out::Float64 = 0.005,
                    prctileceiling::Float64 = 99.5,
                    nsigma::Float64 = 5.0)

Generate a user-friendly (i.e., scaled and thresholded) Gaussian image.

# Description
This method is a wrapper for makegaussim() which generates a more user-friendly
Gaussian image.  That is, this method calls makegaussim(), applies a percentile
ceiling to the results, and then performs a contrast stretch so its pixels
occupy the range [0.0, 1.0].

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`pxsize`: Pixel size of localizations in `smld`. (micrometers)
-`pxsize_out`: Desired output pixel size. (micrometers)
-`prctileceiling`: Upper percentile used to threshold the pixel values of
                   particularly bright pixels.
-`nsigma`: Number of standard deviations from the localization coordinate at
           which we truncate the Gaussian. (Default = 5.0)
"""
function gaussim(smld::SMLMData.SMLD2D, pxsize::Float64;
    pxsize_out::Float64 = 0.005,
    prctileceiling::Float64 = 99.5,
    nsigma::Float64 = 5.0)
    # Determine the requested magnfication factor.
    mag = pxsize / pxsize_out

    # Generate a "raw" Gaussian image (i.e., no scaling or thresholding).
    image = SMLMData.makegaussim(smld; mag = mag, nsigma = nsigma)

    # Apply a percentile ceiling to improve contrast.
    upperbound = StatsBase.percentile(image[:], prctileceiling)
    image[image.>upperbound] .= upperbound

    # Perform a contrast stretch so that the pixels fills the range [0, 1].
    SMLMData.contraststretch!(image)

    return image
end

"""
    image = gaussim(smld::SMLMData.SMLD2D, pxsize::Float64, filename::String;
                    pxsize_out::Float64 = 0.005,
                    prctileceiling::Float64 = 99.5,
                    nsigma::Float64 = 5.0)

Generate and save a user-friendly (i.e., scaled and thresholded) Gaussian image.

# Description
This method is a wrapper for makegaussim() which generates a more user-friendly
Gaussian image.  That is, this method calls makegaussim(), applies a percentile
ceiling to the results, and then performs a contrast stretch so its pixels
occupy the range [0.0, 1.0].  This method of gaussim() saves the image in the
user-specified `filename`.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`pxsize`: Pixel size of localizations in `smld`. (micrometers)
-`filename`: String specifying the location of the output image.
-`pxsize_out`: Desired output pixel size. (micrometers)
-`prctileceiling`: Upper percentile used to threshold the pixel values of
                   particularly bright pixels.
-`nsigma`: Number of standard deviations from the localization coordinate at
           which we truncate the Gaussian. (Default = 5.0)
"""
function gaussim(smld::SMLMData.SMLD2D, pxsize::Float64, filename::String;
    pxsize_out::Float64 = 0.005,
    prctileceiling::Float64 = 99.5,
    nsigma::Float64 = 5.0)
    # Prepare the image.
    image = SMLMData.gaussim(smld, pxsize;
        pxsize_out = pxsize_out,
        prctileceiling = prctileceiling,
        nsigma = nsigma)

    # Create a color image.
    redchannel(x) = min(1.0, 10.0 * x)
    greenchannel(x) = min(1.0, max(10.0 * (x - 0.4), 0.0))
    bluechannel(x) = min(1.0, max(5.0 * (x - 0.7), 0.0))
    imageRGB = Colors.RGB.(redchannel.(image),
        greenchannel.(image),
        bluechannel.(image))

    # Save the image.
    Images.save(filename, imageRGB)

    return image
end

"""
    image = makebinim(coords::Matrix{Float64},
                      datasize::Vector{Float64},
                      mag::Float64 = 20.0)

Make a binary image of the localizations in `coords`.

# Description
This function creates an image of the localizations in `coords` by placing a
hot pixel at the coordinates of each localization.

# Inputs
-`coords`: Localization coordinates. ([y x])
-`datasize`: Size of the data image. ([ysize xsize])
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} binary image of localizations.
"""
function makebinim(coords::Matrix{Float64},
    datasize::Vector{Int},
    mag::Float64 = 20.0)
    # Loop through localizations and add them to our output binary image.
    imagesize = Int.(round.(datasize * mag))
    image = zeros(Float64, imagesize[1], imagesize[2])
    inds = max.(1.0, (coords .- 0.5) * mag)
    inds[:, 1] = min.(imagesize[1], inds[:, 1])
    inds[:, 2] = min.(imagesize[2], inds[:, 2])
    inds = Int.(round.(inds))
    for nn = 1:size(coords, 1)
        image[inds[nn, 1], inds[nn, 2]] = 1.0
    end

    # Normalize the image to sum to 1.0.
    image = image ./ max(sum(image), 1.0)
    if !isapprox(sum(image), 1.0)
        @warn "Image is non-normalizable!  Returning flat image."
        image = ones(Float64, imagesize[1], imagesize[2]) ./ prod(imagesize)
    end

    return image
end

"""
    image = makebinim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)

Make a binary image of the localizations in `smld`.

# Description
This function creates an image of the localizations in `smld` by placing a
hot pixel at the coordinates of each localization.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} binary image of localizations.
"""
function makebinim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)
    return makebinim([smld.y smld.x], smld.datasize, mag)
end

"""
    image = makehistim(coords::Matrix{Float64},
                       datasize::Vector{Int},
                       mag::Float64 = 20.0)

Make a histogram image of the localizations in `coords`.

# Description
This function creates an image of the localizations in `coords` by adding 1.0
to a pixel for each localization present within that pixel.  The final image
is then scaled so that it sums to 1.0.

# Inputs
-`coords`: Localization coordinates. ([y x])
-`datasize`: Size of the data image. ([ysize xsize])
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} histogram image of localizations.
"""
function makehistim(coords::Matrix{Float64},
    datasize::Vector{Int},
    mag::Float64 = 20.0)
    # Loop through localizations and add them to our output image.
    imagesize = Int.(round.(datasize * mag))
    image = zeros(Float64, imagesize[1], imagesize[2])
    inds = max.(1.0, (coords .- 0.5) * mag)
    inds[:, 1] = min.(imagesize[1], inds[:, 1])
    inds[:, 2] = min.(imagesize[2], inds[:, 2])
    inds = Int.(round.(inds))
    for nn = 1:size(coords, 1)
        image[inds[nn, 1], inds[nn, 2]] += 1.0
    end

    # Normalize the image to sum to 1.0.
    image = image ./ max(sum(image), 1.0)
    if !isapprox(sum(image), 1.0)
        @warn "Image is non-normalizable!  Returning flat image."
        image = ones(Float64, imagesize[1], imagesize[2]) ./ prod(imagesize)
    end

    return image
end

"""
    image = makehistim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)

Make a histogram image of the localizations in `smld`.

# Description
This function creates an image of the localizations in `smld` by adding 1.0
to a pixel for each localization present within that pixel.  The final image
is then scaled so that it sums to 1.0.

# Inputs
-`smld`: SMLMData.SMLD2D data structure containing localizations.
-`mag`: Approximate magnfication from data coordinates to SR coordinates.
        (Default = 20.0)

# Outputs
-`image`: Matrix{Float64} histogram image of localizations.
"""
function makehistim(smld::SMLMData.SMLD2D, mag::Float64 = 20.0)
    return makehistim([smld.y smld.x], smld.datasize, mag)
end


text:
using CSV
using DataFrames

OUTPUT_DIR = joinpath(@__DIR__, "output")
MAX_TIME = 300.0

# Read output
output_ooqp = readdir(OUTPUT_DIR)
res_ooqp_ = [f[1:end-9] for f in output_ooqp]
sort!(res_ooqp_)

instance = copy(res_ooqp_)
time_OOQP = MAX_TIME .* ones(length(instance))
success_OOQP = zeros(Int, length(instance))

for (i, f) in enumerate(res_ooqp_)

    # read output file
    s = read(OUTPUT_DIR*"/$(f).mps.ooqp", String)

    # Check output
    if !occursin("*** SUCCESSFUL TERMINATION ***", s)
        # solver experienced issues
        time_OOQP[i] = MAX_TIME
        success_OOQP[i] = 0
    else
        # Parse ouput to get solving time
        s_ = split(s, "QP solved in ")[2]
        t = parse(Float64, split(s_, " seconds.\n")[1])
        if t <= MAX_TIME
            time_OOQP[i] = t
            success_OOQP[i] = 1
        else
            time_OOQP[i] = MAX_TIME
            success_OOQP[i] = 0
        end

    end
end

# Create dataframe of results
df = DataFrame()
df.instance = instance
df.time_OOQP = time_OOQP
df.success_OOQP = success_OOQP

CSV.write("./res_ooqp.csv", df)


text:
# f(x)=A'*X
function FAtransposeX(A,X)
    if size(A)==(1,1)
        return (A[1].*X,nothing)
    elseif size(X)==(1,1)
        return (A'.*X[1],nothing)
    else
        return (A'*X,nothing)
    end
end

function FAtransposeX_inplace(value,auxvalue,A,X)
    if size(A)==(1,1)
        copy!(value,A[1]*X)
    elseif size(X)==(1,1)
        copy!(value,A'*X[1])
    else
        BLAS.gemm!('T','N',1.0,A,X,0.0,value)
    end
end

function DAtransposeX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X)
    if derivativeIDX==1
        if size(A)==(1,1)
            axpy!(1.0,[sum(X.*grad_c)],grad_n)
        elseif size(X)==(1,1)
            axpy!(X[1],grad_c',grad_n)
        else
            BLAS.gemm!('N','T',1.0,X,grad_c,1.0,grad_n) # X*grad_c'
        end
    elseif derivativeIDX==2
        if size(A)==(1,1)
            axpy!(A[1],grad_c,grad_n)
        elseif size(X)==(1,1)
            axpy!(1.0,[sum(A'.*grad_c)],grad_n)
        else
            BLAS.gemm!('N','N',1.0,A,grad_c,1.0,grad_n) # A*grad_c
        end
    end
end

if PROC=="GPU"
#    function FAtransposeX(A::CudaArray,X::CudaArray)
#        if size(A)==(1,1)
#            value=CudaArray(Float64,size(X));
#            copy!(value,X); scale!(A,value)
#            return (value,[])
#        elseif size(X)==(1,1)
#            value=CudaArray(Float64,size(A));
#            copy!(value,A); scale!(X,value) # nb argument converse of Base.scale!
#            return (value,[])
#        else
#            return (CUBLAS.gemm('T','N',A,X),[])
#        end
#    end

    function FAtransposeX_inplace(value::CudaArray,auxvalue,A::CudaArray,X::CudaArray)
        if size(A)==(1,1)
            copy!(value,X); scale!(A,value) # nb argument converse of Base.scale!
        elseif size(X)==(1,1)
            CUBLAS.geam!('T','N',1.0,A,0.0,value,value) # nb argument converse of Base.scale!
            scale!(X,value) # nb argument converse of Base.scale!
        else
            CUBLAS.gemm!('T','N',1.0,A,X,0.0,value)
        end
    end

    function DAtransposeX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray)
        if derivativeIDX==1
            if size(A)==(1,1)
                tmp=CudaArray(Float64,size(X))
                vmult!(1.0,X,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n); free(tmp)
            elseif size(X)==(1,1)
                alphaaxpy!(1.0,X,grad_c,grad_n)
            else
                CUBLAS.gemm!('N','T',1.0,X,grad_c,1.0,grad_n) # X*grad_c'
            end
        elseif derivativeIDX==2
            if size(A)==(1,1)
                alphaaxpy!(1.0,A,grad_c,grad_n)
            elseif size(X)==(1,1)
                tmp=CudaArray(Float64,size(A))
                vmult!(1.0,A,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n);free(tmp)
            else
                CUBLAS.gemm!('N','N',1.0,A,grad_c,1.0,grad_n) # A*grad_c
            end
        end
    end
end

Derivative[FAtransposeX]=DAtransposeX
export FAtransposeX

Inplace[FAtransposeX]=FAtransposeX_inplace

#import Base.*
#*(A::ADnode,B::ADnode)=ADnode(FAX,[A B])

#*(A::Real,B::ADnode)=ADnode(FAX,[ADconst(A) B])
#*(A::ADnode,B::Real)=ADnode(FAX,[A ADconst(B)])

#@gpu *(A::CudaArray,B::CudaArray)=CUBLAS.gemm('N','N',A,B)
#export *

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
Fast (non-allocating) version of norm(x-y,2)^2
"""
function normdiff2(x::AbstractArray{C}, y::AbstractArray{C}) where {
    R <: Real, C <: Union{R, Complex{R}}
}
    s = R(0)
    for i in eachindex(x)
        s += abs2(x[i]-y[i])
    end
    return s
end

"""
Fast (non-allocating) version of norm(x-y,2)
"""
normdiff(x, y) = sqrt(normdiff2(x, y))


text:
function sparklineit(arr::Vector{<:Integer})
    sparkchars = '\u2581':'\u2588'
    dyn = length(sparkchars)
    lo, hi = extrema(arr)
    b = @. max(ceil(Int, dyn * (arr - lo) / (hi - lo)), 1)
    return join(sparkchars[b])
end

v = rand(0:10, 10)
println("$v → ", sparklineit(v))
v = 10rand(10)
println("$(round.(v, 2)) → ", sparklineit(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#! /usr/bin/julia

# Rosetta Code, Sparkline in unicode

function sparklineit(a)
    const sparkchars = '\u2581':'\u2588'
    const dyn = length(sparkchars)
    (lo, hi) = extrema(a)
    b = max(iceil(dyn*(a-lo)/(hi-lo)), 1)
    return join(sparkchars[b], "")
end

function getnumbers(s)
    a = split(s, r"[,,\s]+")
    a = try
        map(parseint, a)
    catch
        map(parsefloat, a)
    end
end

test = getnumbers("1 2 3 4 5 6 7 8 7 6 5 4 3 2 1")
println(test, " => ", sparklineit(test))

test = getnumbers("1.5, 0.5 3.5, 2.5 5.5, 4.5 7.5, 6.5")
println(test, " => ", sparklineit(test))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#=
Work out the first ten digits of the sum of the following one-hundred 50-digit numbers.
=#

x = [37107287533902102798797998220837590246510135740250,
46376937677490009712648124896970078050417018260538,
74324986199524741059474233309513058123726617309629,
91942213363574161572522430563301811072406154908250,
23067588207539346171171980310421047513778063246676,
89261670696623633820136378418383684178734361726757,
28112879812849979408065481931592621691275889832738,
44274228917432520321923589422876796487670272189318,
47451445736001306439091167216856844588711603153276,
70386486105843025439939619828917593665686757934951,
62176457141856560629502157223196586755079324193331,
64906352462741904929101432445813822663347944758178,
92575867718337217661963751590579239728245598838407,
58203565325359399008402633568948830189458628227828,
80181199384826282014278194139940567587151170094390,
35398664372827112653829987240784473053190104293586,
86515506006295864861532075273371959191420517255829,
71693888707715466499115593487603532921714970056938,
54370070576826684624621495650076471787294438377604,
53282654108756828443191190634694037855217779295145,
36123272525000296071075082563815656710885258350721,
45876576172410976447339110607218265236877223636045,
17423706905851860660448207621209813287860733969412,
81142660418086830619328460811191061556940512689692,
51934325451728388641918047049293215058642563049483,
62467221648435076201727918039944693004732956340691,
15732444386908125794514089057706229429197107928209,
55037687525678773091862540744969844508330393682126,
18336384825330154686196124348767681297534375946515,
80386287592878490201521685554828717201219257766954,
78182833757993103614740356856449095527097864797581,
16726320100436897842553539920931837441497806860984,
48403098129077791799088218795327364475675590848030,
87086987551392711854517078544161852424320693150332,
59959406895756536782107074926966537676326235447210,
69793950679652694742597709739166693763042633987085,
41052684708299085211399427365734116182760315001271,
65378607361501080857009149939512557028198746004375,
35829035317434717326932123578154982629742552737307,
94953759765105305946966067683156574377167401875275,
88902802571733229619176668713819931811048770190271,
25267680276078003013678680992525463401061632866526,
36270218540497705585629946580636237993140746255962,
24074486908231174977792365466257246923322810917141,
91430288197103288597806669760892938638285025333403,
34413065578016127815921815005561868836468420090470,
23053081172816430487623791969842487255036638784583,
11487696932154902810424020138335124462181441773470,
63783299490636259666498587618221225225512486764533,
67720186971698544312419572409913959008952310058822,
95548255300263520781532296796249481641953868218774,
76085327132285723110424803456124867697064507995236,
37774242535411291684276865538926205024910326572967,
23701913275725675285653248258265463092207058596522,
29798860272258331913126375147341994889534765745501,
18495701454879288984856827726077713721403798879715,
38298203783031473527721580348144513491373226651381,
34829543829199918180278916522431027392251122869539,
40957953066405232632538044100059654939159879593635,
29746152185502371307642255121183693803580388584903,
41698116222072977186158236678424689157993532961922,
62467957194401269043877107275048102390895523597457,
23189706772547915061505504953922979530901129967519,
86188088225875314529584099251203829009407770775672,
11306739708304724483816533873502340845647058077308,
82959174767140363198008187129011875491310547126581,
97623331044818386269515456334926366572897563400500,
42846280183517070527831839425882145521227251250327,
55121603546981200581762165212827652751691296897789,
32238195734329339946437501907836945765883352399886,
75506164965184775180738168837861091527357929701337,
62177842752192623401942399639168044983993173312731,
32924185707147349566916674687634660915035914677504,
99518671430235219628894890102423325116913619626622,
73267460800591547471830798392868535206946944540724,
76841822524674417161514036427982273348055556214818,
97142617910342598647204516893989422179826088076852,
87783646182799346313767754307809363333018982642090,
10848802521674670883215120185883543223812876952786,
71329612474782464538636993009049310363619763878039,
62184073572399794223406235393808339651327408011116,
66627891981488087797941876876144230030984490851411,
60661826293682836764744779239180335110989069790714,
85786944089552990653640447425576083659976645795096,
66024396409905389607120198219976047599490197230297,
64913982680032973156037120041377903785566085089252,
16730939319872750275468906903707539413042652315011,
94809377245048795150954100921645863754710598436791,
78639167021187492431995700641917969777599028300699,
15368713711936614952811305876380278410754449733078,
40789923115535562561142322423255033685442488917353,
44889911501440648020369068063960672322193204149535,
41503128880339536053299340368006977710650566631954,
81234880673210146739058568557934581403627822703280,
82616570773948327592232845941706525094512325230608,
22918802058777319719839450180888072429661980811197,
77158542502016545090413245809786882778948721859617,
72107838435069186155435662884062257473692284509516,
20849603980134001723930671666823555245252804609722,
53503534226472524250874054075591789781264330331690]

@time println(string(sum(x))[1:10])

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using WatchJuliaBurn
using Test

@testset "WatchJuliaBurn.jl" begin
    @test_throws ArgumentError c╯°□°ↄ╯(ArgumentError("Great Success"))
end


text:
module HyperbolicCrossApprox

import SmolyakApprox: normalize_node,
                      scale_nodes,
                      combine_nodes

import ChebyshevApprox: chebyshev_polynomial,
                        derivative_of_chebyshev_polynomial,
                        chebyshev_nodes,
                        chebyshev_extrema,
                        chebyshev_extended,
                        vertesi_nodes

using ThreadPools

include("hyperbolic_cross_weights.jl")
include("hyperbolic_cross_evaluate.jl")
include("hyperbolic_cross_grid.jl")
include("generate_multi_index.jl")
include("hyperbolic_cross_derivative.jl")
include("hyperbolic_cross_gradient.jl")
include("determine_grid_size.jl")
include("hyperbolic_cross_inverse_interpolation_matrix.jl")

export chebyshev_nodes,
       chebyshev_extrema,
       chebyshev_extended,
       vertesi_nodes,
       hyperbolic_cross_grid,
       hyperbolic_cross_weights,
       hyperbolic_cross_weights_threaded,
       hyperbolic_cross_inverse_interpolation_matrix,
       hyperbolic_cross_inverse_interpolation_matrix_threaded,
       hyperbolic_cross_polynomial,
       hyperbolic_cross_evaluate,
       hyperbolic_cross_derivative,
       hyperbolic_cross_gradient

end