text:
import Base: convert, getindex, setindex!, similar

## Constructors and converters
## (special methods for AbstractArray{Nullable}, to avoid wrapping Nullable inside Nullable)

NullableCategoricalArray{T, N}(::Type{Nullable{T}}, dims::NTuple{N,Int}; ordered=false) =
    NullableCategoricalArray{T, N}(zeros(DefaultRefType, dims), CategoricalPool(ordered))
NullableCategoricalArray{T}(::Type{Nullable{T}}, dims::Int...; ordered=false) =
    NullableCategoricalArray(T, dims; ordered=ordered)

@compat (::Type{NullableCategoricalArray{Nullable{T}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                       ordered=false) =
    NullableCategoricalArray(zeros(R, dims), CategoricalPool{T, R}(ordered))
@compat (::Type{NullableCategoricalArray{Nullable{T}, N}}){T, N}(dims::NTuple{N,Int};
                                                                 ordered=false) =
    NullableCategoricalArray{T}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m,); ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int, n::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m, n); ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int, n::Int, o::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m, n, o); ordered=ordered)

@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T, R}}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T}}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T, R}}, N}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T}}, N}}){T, N}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N}(dims; ordered=ordered)
# @compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue}, N}}){N}(dims::NTuple{N,Int};
#                                                                               ordered=false) =
#     NullableCategoricalArray{String, N}(dims; ordered=ordered)
# @compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue}}}){N}(dims::NTuple{N,Int};
#                                                                            ordered=false) =
#     NullableCategoricalArray{String, N}(dims; ordered=ordered)

if VERSION >= v"0.5.0-dev"
    NullableCategoricalVector{T}(::Type{Nullable{T}}, m::Integer; ordered=false) =
      NullableCategoricalArray{T}((m,); ordered=ordered)
end
@compat (::Type{NullableCategoricalVector{Nullable{T}}}){T}(m::Int; ordered=false) =
    NullableCategoricalArray{T}((n,); ordered=ordered)

if VERSION >= v"0.5.0-dev"
    NullableCategoricalMatrix{T}(::Type{Nullable{T}}, m::Int, n::Int; ordered=false) =
      NullableCategoricalArray{T}((m, n); ordered=ordered)
end
@compat (::Type{NullableCategoricalMatrix{Nullable{T}}}){T}(m::Int, n::Int; ordered=false) =
    NullableCategoricalArray{T}((m, n); ordered=ordered)

@compat (::Type{NullableCategoricalArray}){T<:Nullable}(A::AbstractArray{T};
                                                        ordered=_ordered(A)) =
    NullableCategoricalArray{eltype(T)}(A, ordered=ordered)
@compat (::Type{NullableCategoricalVector}){T<:Nullable}(A::AbstractVector{T};
                                                         ordered=_ordered(A)) =
    NullableCategoricalVector{eltype(T)}(A, ordered=ordered)
@compat (::Type{NullableCategoricalMatrix}){T<:Nullable}(A::AbstractMatrix{T};
                                                         ordered=_ordered(A)) =
    NullableCategoricalMatrix{eltype(T)}(A, ordered=ordered)

function NullableCategoricalArray{T, N}(A::AbstractArray{T, N},
                                        missing::AbstractArray{Bool};
                                        ordered=false)
    res = NullableCategoricalArray{T, N}(size(A); ordered=ordered)
    @inbounds for (i, x, m) in zip(eachindex(res), A, missing)
        res[i] = ifelse(m, Nullable{T}(), x)
    end

    if method_exists(isless, (T, T))
        levels!(res, sort(levels(res)))
    end

    res
end

if VERSION >= v"0.5.0-dev"
    NullableCategoricalVector{T}(A::AbstractVector{T},
                                 missing::AbstractVector{Bool};
                                 ordered=false) =
        NullableCategoricalArray(A, missing; ordered=ordered)
    NullableCategoricalMatrix{T}(A::AbstractMatrix{T},
                                 missing::AbstractMatrix{Bool};
                                 ordered=false) =
        NullableCategoricalArray(A, missing; ordered=ordered)
end

function getindex(A::NullableCategoricalArray, i::Int)
    j = A.refs[i]
    S = eltype(eltype(A))
    j > 0 ? Nullable{S}(A.pool[j]) : Nullable{S}()
end

function setindex!(A::NullableCategoricalArray, v::Nullable, i::Int)
    if isnull(v)
        A.refs[i] = 0
    else
        A[i] = get(v)
    end
end

levels!(A::NullableCategoricalArray, newlevels::Vector; nullok=false) = _levels!(A, newlevels, nullok=nullok)


text:
module SellmeierFit

using LsqFit
using CSV, DataFrames, Downloads  # used in read.jl
using AbbreviatedTypes
using SimpleConstants

include("read.jl")
include("model.jl")
include("fit.jl")

end # module

ext: jl


text:
module Tabu

    using Tsp.Model

    """
        find_optimal(cities, [kwargs])

    ## Arguments
    * cities - list of cities

    ## Keyword arguments
    * rounds - max number of iterations to perform [default: 100]
    * restart_after - number of rounds without progress before a jump to new solution space region [default: 5]
    * quit_after - max number of jumps before ending the search [default: 15]
    * min_tabu_duration - [default: sqrt(legth(cities))]
    """
    function find_optimal(
        grid::Grid;
        rounds::Int             = 1000,
        f_neighbours            = .2,
        n_neighbours            = round(Int, f_neighbours * length(grid.cities)),
        jump_radius             = div(length(grid.cities), 10),
        restart_after::Int      = 5,
        quit_after::Int         = 15,
        min_tabu_duration::Int  = min(div(rounds, 2), round(Int, sqrt(length(grid.cities)))),
        max_tabu_duration::Int  = 2 * min_tabu_duration
        )
        tic()
        #= println("""
            Initializing search:
                Rounds: $rounds
                Neighbourood size: $n_neighbours
                Neighbourhood radius: $jump_radius
                Restart after: $restart_after rounds
                Restart limit: $quit_after
                Tabu duration: $min_tabu_duration:$max_tabu_duration
                Jump radius: $jump_radius
            """)*/=#

        cities = grid.cities
        best::Route = nearest_neighbour(cities) # best solution so far
        best_round = 0
        origin::Route = best # solution from which the search currently progresses

        rounds_without_progress::Int = 0
        rounds_with_progress::Int = 0
        rounds_without_progress_strict::Int = 0
        jumps::Int64 = 0

        # save initial argument values
        initial_jump_radius = jump_radius
        initial_restart_after = restart_after
        initial_n_neighbours = n_neighbours
        initial_score = origin.score
        n_neighbours = min(length(cities)^2, n_neighbours)

        max_neighbours = min(length(cities) >= 10000 ? 5000 : n_neighbours * 2)

        max_restart_after = 2 * initial_restart_after

        # @debug "Initial solution: %10f (%.4fs)\n", score!(origin), toq()
        last_loop= 0 # last loop's time in seconds
        # The main loop
        time_elapsed = 0.0
        rounds_in_terrible_neighbourhood = 0
        for cur_round in 1:rounds
            tic()
            if rounds_without_progress_strict >= 10
                # @debug("[%5d] Not enough progress, forcing restart\n", cur_round)
                jumps += 1
                restart_after = initial_restart_after
                n_neighbours = initial_n_neighbours
                rounds_without_progress_strict = 0
                rounds_with_progress = 0
                rounds_without_progress = restart_after # force restart
            end
            # ignore terrible solutions faster
            if rounds_in_terrible_neighbourhood == 0 && (origin.score - initial_score)/initial_score * 100 > 100.0
                # @debug("Terrible solution (%+8.5f), not staying for long\n", (origin.score - initial_score)/initial_score * 100 > 100.0)
                rounds_in_terrible_neighbourhood += 1
                restart_after = min(0, 5 - floor(Int, origin.score - initial_score)/initial_score)
            end
            # change origin if there was no progress in the last restart_after rounds
            if rounds_without_progress + rounds_in_terrible_neighbourhood >= restart_after
                # println("[$cur_round] No progress for $rounds_without_progress rounds")
                if jumps >= quit_after # too many jumps - accept current best solution
                    # @debug("[%5d] Restart limit reached - stopping the search\n", cur_round)
                    break
                end
                origin = jump(origin, jump_radius) # get_neighbours(origin, 1, jump_radius*100, grid)[1][1]
                # println("[$cur_round] Restarting.\tNew origin: $(score!(origin)). Restarts: $jumps/$quit_after")
                # jump_radius += 1
                jumps += 1
                rounds_without_progress = min(restart_after, floor(Int, time_elapsed/240.0))
                rounds_with_progress = 0
                rounds_in_terrible_neighbourhood = 0
                restart_after = initial_restart_after
                n_neighbours = initial_n_neighbours
                # @printf("[%5d] Restarting: new origin: %10f (%+8.5f%%)\tRestarts: %d/%d, Time elapsed: %fs (+%d)\n",
                    # cur_round, score!(origin), (origin.score - initial_score)/initial_score * 100,
                    # jumps, quit_after, time_elapsed, rounds_without_progress)
            end

            if rounds_with_progress >= 5
                if restart_after < initial_restart_after
                    restart_after = initial_restart_after
                end
                if restart_after < max_restart_after
                    restart_after += round(Int, .5*restart_after)
                end
                if n_neighbours < max_neighbours
                    n_neighbours += round(Int, .1*n_neighbours)
                end
                rounds_without_progress = 0
                rounds_with_progress = 0
                # @printf("[%5d] Nice neighbourhood, extending search. New params:\n\tNeighbourood size: %d\n\tRestarting after: %d\n",
                    # cur_round, n_neighbours, restart_after)
                # println("[$cur_round] Nice neighbourhood, extending search. New params:
                #         Neighbourood size: $n_neighbours
                #         Restart after $restart_after rounds")
            end

            # tic() # neighbours
            neighbours = Task(()->get_neighbours(origin, n_neighbours, grid)) # "coroutine" to avoid allocating all the neighbours at once
            best_candidate::Route, best_i, best_j = first(neighbours)

            # search through the neighbours
            # tic()
            for (candidate, i, j) in neighbours

                if cities[i].tabu >= cur_round || cities[j].tabu >= cur_round
                    # check aspiration
                    penalty = 0.5 * abs((cities[i].tabu  + cities[j].tabu - cur_round*2))
                    # @debug("tabu: %16f, %d, %d, penalty = %d | best = %16f\n", score!(candidate),i, j, penalty, best.score)
                    if score!(candidate) + penalty < score!(best)
                        # @printf("[%5d] Found a great solution, ignoring tabu [%d,%d]: %10f (%+8.5f%%), swap: (%d, %d)\n",
                        #      cur_round, cities[i].tabu, cities[j].tabu, candidate.score, (candidate.score - initial_score)/initial_score * 100, i, j)
                        best_candidate, best_i, best_j = candidate, i, j
                        best = candidate
                        continue
                    end
                elseif score!(candidate) < score!(best_candidate)
                    best_candidate, best_i, best_j = candidate, i, j
                else
                    make_tabu(cities[i], min(
                        cur_round + div(min_tabu_duration, 2),
                        floor(Int, cur_round + rand(min_tabu_duration:max_tabu_duration) * (origin.score - initial_score)/initial_score * 100))
                    )
                    make_tabu(cities[j], min(
                        cur_round + div(min_tabu_duration, 2),
                        floor(Int, cur_round + rand(min_tabu_duration:max_tabu_duration) * (origin.score - initial_score)/initial_score * 100))
                    )
                end
            end

            # println("[Time][$cur_round] $n_neighbours neighbours checked in $(toq())s")
            # check best neighbour with the origin
            if score!(best_candidate) < score!(origin)
                # @debug("[%5d] New best neighbour: %10f (%+8.5f%%), swap: (%d, %d)\n", cur_round, best_candidate.score, (best_candidate.score - initial_score)/initial_score * 100, best_i, best_j)
                rounds_without_progress = 0
                origin = best_candidate

                if score!(origin) < score!(best)
                    # println("[$cur_round] New best solution: $(origin.score) ($((origin.score - initial_score)/initial_score * 100)%)")
                    # @info @printf("[%5d] New best solution: %10f (%+8.5f%%)\n", cur_round, origin.score, (origin.score - initial_score)/initial_score * 100)
                    best = origin
                    best_round = cur_round
                    rounds_without_progress_strict = 0
                else
                    rounds_without_progress_strict += 1
                end
            else
                rounds_without_progress += 1
            end
            round_time = toq()
            time_elapsed += round_time
            if (round_time > 0.1 || cur_round % 50 == 0)
                # @info("[%5d] %7.5fs\tbest: %10f\tround best: %10f (%+8.5f%%)\n",
                #      cur_round, round_time, best.score, best_candidate.score,
                #      (best_candidate.score - initial_score)/initial_score * 100)
            end
        end
        # @info @printf("Best solution found in round #%d\n", best_round)
        return best
    end #find_optimal

    function get_neighbours(origin::Route, n_neighbours::Int, grid::Grid)
        # neighbours = Array(Tuple{Route, Int, Int}, n_neighbours)
        n = 0
        cities = grid.cities
        while n < n_neighbours
            neighbour = Route(origin)
            i = rand(2:length(origin)-1)
            city = grid.cities[i]
            chunk = grid.chunks[city.chunk_x, city.chunk_y]
            if length(chunk.cities) < 2
                # println("WARNING: requested swap in a single-city chunk, ignoring")
                continue
            end
            j_id = rand(chunk.cities).id # id of city to swap ith city on route with
            # println("[swap: $i] indices: i=$i, j=$j")
            while j_id == city.id || j_id == 1 # ensure there is a swap
                # if canary >= 100
                #     println("WARNING: stuck selecting city to swap for 100 rounds")
                #     break
                # end
                j_id = rand(chunk.cities).id
                # canary += 1
            end
            j = findfirst(city->city.id == j_id, neighbour.cities)
            # middle = i < j ? (i:j) : (i:-1:j)
            # println("old: $neighbour")
            # println("i=$i, j=$j")
            if i < j
                neighbour.score -= distance(neighbour.cities[i], neighbour.cities[i-1]) + distance(neighbour.cities[j], neighbour.cities[j+1])
                neighbour.cities[i:j] = neighbour.cities[reverse(i:j)]
                neighbour.score += distance(neighbour.cities[i], neighbour.cities[i-1]) + distance(neighbour.cities[j], neighbour.cities[j+1])
            else
                neighbour.score -= distance(neighbour.cities[j], neighbour.cities[j-1]) + distance(neighbour.cities[i], neighbour.cities[i+1])
                neighbour.cities[j:i] = neighbour.cities[reverse(j:i)]
                neighbour.score += distance(neighbour.cities[j], neighbour.cities[j-1]) + distance(neighbour.cities[i], neighbour.cities[i+1])
            end
            n += 1
            produce((neighbour, city.id, j_id)) # avoids allocating a shitload of paths for no good reason
        end
    end #get_neighbours

    function jump(origin::Route, jump_radius::Int)
        new_route = Route(origin)
        i = rand(2+jump_radius:length(origin)-1-jump_radius)
        j = rand(2:length(origin)-1)
        # println("$(i+jump_radius), $(length(origin)-1-jump_radius)")
        new_route.cities[i-jump_radius:i+jump_radius] = shuffle(new_route.cities[i-jump_radius:i+jump_radius])
        return new_route
    end #jump
end


text:
# Autogenerated wrapper script for LAZperf_jll for i686-linux-musl-cxx03
export liblazperf

JLLWrappers.@generate_wrapper_header("LAZperf")
JLLWrappers.@declare_library_product(liblazperf, "liblazperf.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liblazperf,
        "lib/liblazperf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module UncertaintyQuantification

using Accessors
using Bootstrap
using DataFrames
using Dates
using Dierckx
using Distributed
using FiniteDifferences
using Formatting
using LinearAlgebra
using Mustache
using Primes
using QuasiMonteCarlo
using Random
using Reexport

@reexport using Distributions

import Base: rand, names, copy, run
import Statistics: mean

abstract type UQType end

abstract type UQInput <: UQType end
abstract type UQModel <: UQType end

abstract type DeterministicUQInput <: UQInput end
abstract type RandomUQInput <: UQInput end

abstract type Copula <: UQType end

abstract type AbstractMonteCarlo end
abstract type AbstractQuasiMonteCarlo <: AbstractMonteCarlo end

# Types
export AbstractMonteCarlo
export AbstractQuasiMonteCarlo
export Copula
export DeterministicUQInput
export RandomUQInput
export UQInput
export UQModel
export UQType

# Structs
export ExternalModel
export Extractor
export GaussianCopula
export HaltonSampling
export JointDistribution
export LatinHypercubeSampling
export LatticeRuleSampling
export LineSampling
export Model
export MonteCarlo
export Parameter
export PolyharmonicSpline
export RandomVariable
export SobolSampling
export Solver
export SubSetSimulation

# Methods
export calc
export count_rvs
export dimensions
export evaluate!
export gradient
export gradient_in_standard_normal_space
export mean
export probability_of_failure
export qmc_samples
export rand
export sample
export sobolindices
export to_copula_space
export to_physical_space!
export to_standard_normal_space
export to_standard_normal_space!

include("inputs/inputs.jl")
include("inputs/parameter.jl")
include("inputs/randomvariable.jl")
include("inputs/jointdistribution.jl")

include("inputs/copulas/gaussian.jl")

include("solvers/solver.jl")
include("solvers/extractor.jl")

include("models/externalmodel.jl")
include("models/model.jl")
include("models/polyharmonicspline.jl")

include("sensitivity/gradient.jl")

include("simulations/linesampling.jl")
include("simulations/montecarlo.jl")
include("simulations/subset.jl")

include("reliability/probabilityoffailure.jl")
include("sensitivity/sobolindices.jl")

end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using PoissonRandom
import Distributions
using Test, Statistics

n_tsamples = 10^5

function test_samples(rand_func,
                      distr::Distributions.DiscreteUnivariateDistribution,
                      n::Int;                                # number of samples to generate
                      q::Float64=1.0e-8,                     # confidence interval, 1 - q as confidence
                      verbose::Bool=false)                   # show intermediate info (for debugging)

    # The basic idea
    # ------------------
    #   Generate n samples, and count the occurences of each value within a reasonable range.
    #   For each distinct value, it computes an confidence interval of the counts
    #   and checks whether the count is within this interval.
    #
    #   If the distribution has a bounded range, it also checks whether
    #   the samples are all within this range.
    #
    #   By setting a small q, we ensure that failure of the tests rarely
    #   happen in practice.
    #
    λ = distr.λ
    n > 1 || error("The number of samples must be greater than 1.")
    0.0 < q < 0.1 || error("The value of q must be within the open interval (0.0, 0.1).")

    # determine the range of values to examine
    vmin = minimum(distr)
    vmax = maximum(distr)

    rmin = floor(Int,quantile(distr, 0.00001))::Int
    rmax = floor(Int,quantile(distr, 0.99999))::Int
    m = rmax - rmin + 1  # length of the range
    p0 = Distributions.pdf.((distr,), rmin:rmax)  # reference probability masses
    @assert length(p0) == m

    # determine confidence intervals for counts:
    # with probability q, the count will be out of this interval.
    #
    clb = Vector{Int}(undef,m)
    cub = Vector{Int}(undef,m)
    for i = 1:m
        bp = Distributions.Binomial(n, p0[i])
        clb[i] = floor(Int,quantile(bp, q/2))
        cub[i] = ceil(Int,Distributions.cquantile(bp, q/2))
        @assert cub[i] >= clb[i]
    end

    # generate samples
    samples = [rand_func(λ) for i in 1:n]
    @assert length(samples) == n

    # scan samples and get counts
    cnts = zeros(Int, m)
    for i = 1:n
        @inbounds si = samples[i]
        if rmin <= si <= rmax
            cnts[si - rmin + 1] += 1
        else
            vmin <= si <= vmax ||
                error("Sample value out of valid range.")
        end
    end

    # check the counts
    for i = 1:m
        verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts[i])")
        clb[i] <= cnts[i] <= cub[i] ||
            error("The counts are out of the confidence interval.")
    end
    return samples
end

println("testing count random sampler")
for λ in [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0]
  test_samples(PoissonRandom.count_rand, Distributions.Poisson(λ), n_tsamples)
end

println("testing ad random sampler")
for λ in [5.0, 10.0, 15.0, 20.0, 30.0]
  test_samples(PoissonRandom.ad_rand, Distributions.Poisson(λ), n_tsamples)
end

println("testing mixed random sampler")
for λ in [5.0, 10.0, 15.0, 20.0, 30.0]
  test_samples(pois_rand, Distributions.Poisson(λ), n_tsamples)
end


text:
function prec_mkl_ilu0( A::SparseMatrixCSC{Float64,Int64} )
  NNZ = size(A.nzval)[1]
  #println("NNZ = ", NNZ)
  bilu0 = zeros(NNZ)
  #for i = 1:size(A.colptr)[1]
  #  @printf("%8d %8d\n", i, A.colptr[i] )
  #end
  #for i = 1:size(A.rowval)[1]
  #  @printf("%8d %8d %18.10f\n", i, A.rowval[i], A.nzval[i] )
  #end
  ccall( (:mkl_ilu0, "../LF/extlibs/libmkl.so"), Void,
         (Int, Int, Ptr{Float64}, Ptr{Int}, Ptr{Int}, Ptr{Float64}),
         A.n, NNZ, A.nzval, A.colptr, A.rowval, bilu0 )
  return SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, bilu0)
end


text:
function center_test(x1, x2, x3, y1, y2, y3)
    arr1 = randn((x1, x2, x3))
    arr2 = zeros((y1, y2, y3))

    center_set!(arr2, arr1)
    arr3 = center_extract(arr2, (x1, x2, x3))
    @test arr1 ≈ arr3
end

 # test center set and center extract methods
@testset "center methods" begin
    center_test(4, 4, 4, 6,7,4)
    center_test(5, 4, 4, 7, 8, 4)
    center_test(5, 4, 4, 8, 8, 8)
    center_test(6, 4, 4, 7, 8, 8)

    @test 1 == center_pos(1)
    @test 2 == center_pos(2)
    @test 2 == center_pos(3)
    @test 3 == center_pos(4)
    @test 3 == center_pos(5)
    @test 513 == center_pos(1024)

    @test PhysicalOptics.get_indices_around_center((5), (2)) == (2, 3)
    @test PhysicalOptics.get_indices_around_center((5), (3)) == (2, 4)
    @test PhysicalOptics.get_indices_around_center((4), (3)) == (2, 4)
    @test PhysicalOptics.get_indices_around_center((4), (2)) == (2, 3)
end

@testset "rr methods" begin
    out = [1.4142135623730951 1.0 1.4142135623730951; 1.0 0.0 1.0; 1.4142135623730951 1.0 1.4142135623730951]
    @test out ≈ rr((3,3))
    @test [0] ≈ rr((1, 1))
    @test [2,1,0,1,2] ≈ rr((5, 1))
    @test [3,2,1,0,1,2] ≈ rr((6, 1))

    out = [ sqrt(2)  1.0; 1.0      0.0]
    @test out == rr((2,2))

    @test rr((4, 4)) == [2.8284271247461903 2.23606797749979 2.0 2.23606797749979; 2.23606797749979 1.4142135623730951 1.0 1.4142135623730951; 2.0 1.0 0.0 1.0; 2.23606797749979 1.4142135623730951 1.0 1.4142135623730951]

    @test rr((4, 4), L = 2) == [1.4142135623730951 1.118033988749895 1.0 1.118033988749895; 1.118033988749895 0.7071067811865476 0.5 0.7071067811865476; 1.0 0.5 0.0 0.5; 1.118033988749895 0.7071067811865476 0.5 0.7071067811865476]

    @test rr((5, 5), L = 4) == [2.8284271247461903 2.23606797749979 2.0 2.23606797749979 2.8284271247461903; 2.23606797749979 1.4142135623730951 1.0 1.4142135623730951 2.23606797749979; 2.0 1.0 0.0 1.0 2.0; 2.23606797749979 1.4142135623730951 1.0 1.4142135623730951 2.23606797749979; 2.8284271247461903 2.23606797749979 2.0 2.23606797749979 2.8284271247461903]

    @test rr((5, 5), L = (4, 10)) == [5.385164807134504 3.2015621187164243 2.0 3.2015621187164243 5.385164807134504; 5.0990195135927845 2.692582403567252 1.0 2.692582403567252 5.0990195135927845; 5.0 2.5 0.0 2.5 5.0; 5.0990195135927845 2.692582403567252 1.0 2.692582403567252 5.0990195135927845; 5.385164807134504 3.2015621187164243 2.0 3.2015621187164243 5.385164807134504]

    @test rr((5, 6), L = (20, 10)) == [11.180339887498949 10.540925533894598 10.137937550497032 10.0 10.137937550497032 10.540925533894598; 7.0710678118654755 6.009252125773315 5.270462766947299 5.0 5.270462766947299 6.009252125773315; 5.0 3.3333333333333335 1.6666666666666667 1.776356839400252e-16 1.6666666666666665 3.333333333333333; 7.0710678118654755 6.009252125773315 5.270462766947299 5.0 5.270462766947299 6.009252125773315; 11.180339887498949 10.540925533894598 10.137937550497032 10.0 10.137937550497032 10.540925533894598]

    out = [1.7320508075688772 1.4142135623730951; 1.4142135623730951 1.0; 1.4142135623730951 1.0; 1.0 0.0]
    out = reshape(out, (2,2,2))
    @test out ≈ rr_3D((2,2,2))

end

@testset "Test fftpos" begin
    @test fftpos(1, 10) == -0.5:0.1:0.4

    @test collect(fftpos(0.01, 3)) ≈ [-0.005, 0.0, 0.005]

    @test fftpos(0.1, 20) ≈ -0.05:0.005:0.045

    @test fftpos(0.1, 21) == -0.05:0.005:0.05

    for N = 2:30
        @test length(fftpos(123, N)) == N
    end
end

@testset "Test normabs2" begin
    x = [1.0, 1.0]
    @test normabs2(x) == x
    x = [1.0 3.0; 4.0 123123123.2]
    @test normabs2(x) ≈ abs2.(x) / maximum(abs2.(x))
end

@testset "Test apply_rot_symmetry" begin

    function test_asr(f, xpos, xpos_out, ypos_out, rtol)
        out_ref = [f(sqrt(x^2 + y^2)) for x in xpos_out, y in ypos_out]

        ev_1D = [f(x) for x in xpos]

        out = PhysicalOptics.apply_rot_symmetry(ev_1D, xpos, xpos_out, ypos_out)

        @test ≈(out_ref, out, rtol=rtol)
    end

    x = range(0, 2 * √2, length=250)
    x_out = range(0, 2, length=20)
    test_asr(sin, x, x_out, x_out, 1e-9)

    x = range(0, 10 * √2, length=1000)
    x_out = range(0, 10, length=20)
    test_asr(SpecialFunctions.jinc, x, x_out, x_out, 1e-7)

end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
@testset "Datasets availability" begin
    package_directory = joinpath(dirname(@__FILE__), "..", "data")

    dfs_found = vcat(map(readdir(package_directory)) do directory
        fns = readdir(joinpath(package_directory, directory))
        pkgs = fill(directory, length(fns))
        sets = replace.(fns, Ref(r"(\.(rda|csv|gz))+$" => ""))
        DataFrame(Package = pkgs, Dataset = sets)
    end...)

    dfs_documented = RDatasets.datasets()

    @test size(dfs_found, 1) == size(dfs_documented, 1)
    @test sort(dfs_documented[[:Package, :Dataset]]) == sort(dfs_found)

    # test caching works
    @test dfs_documented === RDatasets.datasets()
end

max_stars_repo_name: blockchainhelppro/Wadren-Analysis

max_issues_repo_name: blockchainhelppro/Wadren-Analysis

max_forks_repo_name: blockchainhelppro/Wadren-Analysis


text:
module WebDisplay
    using Logging

    include("web.jl")
    include("display.jl")

    function __init__()
        host = get(ENV, "WEB_DISPLAY_HOST", "0.0.0.0")
        port = try parse(Int, ENV["WEB_DISPLAY_PORT"]) catch; rand(6000:9000) end
        @async with_logger(SimpleLogger(stderr, Logging.Error)) do
            _web.listen(host, port)
        end
        pushdisplay(_display)
        @info "WebDisplay listening on $host:$port"
    end
end


text:
module SDMHelper

using DelimitedFiles
using PWAHelper

using LinearAlgebra

export SDM_to_pars, pars_to_SDM
export write_cmatrix, read_cmatrix
export constract_values
export enlarge_with_zeros!

function pars_to_SDM(pars, Bmat, block_masks)
    Tmap = get_parameter_map(block_masks, size(Bmat,1))
    pblocks = make_pblock_inds(block_masks)
    Nw = size(Bmat,1)
    nonormSDM = fill(0.0im, Nw, Nw)
    for bl in pblocks
        pars_bl = fill(0.0, length(pars))
        pars_bl[bl] .= pars[bl]
        sp = shrnk(pars_bl,Tmap)
        nonormSDM .+= [conj(sp[i])*sp[j] for i in 1:Nw, j in 1:Nw]
    end
    diagSHBM = real.(diag(Bmat))
    lSDM = [nonormSDM[i,j]*sqrt(diagSHBM[i]*diagSHBM[j]) for i in 1:Nw, j in 1:Nw]
    return lSDM
end

function SDM_to_pars(SDM, Bmat, block_masks)
    pars = sqrt.(diag(SDM))
    for i in 2:size(pars,1)
        el = SDM[i-1,i]
        (el != 0.0) && (pars[i] *= el/abs(el) * pars[i-1]/abs(pars[i-1]))
    end
    pars = pars ./ sqrt.(diag(Bmat))
    # convert to re im
    Tmap = get_parameter_map(block_masks, size(Bmat,1))
    fpars = fill(0.0, size(Tmap,2))
    for i in 1:size(Tmap,2)
        if Tmap[1,i] != 0
            (fpars[i] == 0.0) && (fpars[i] = real(pars[Tmap[1,i]]))
        else
            (fpars[i] == 0.0) && (fpars[i] = imag(pars[Tmap[2,i]]))
        end
    end
    return fpars
end

function write_cmatrix(sdm, fout)
    #io = open(fout,"w")
    #write(io,size(sdm,1))
    #write(io,size(sdm,1)+size(sdm,1))
    #write(io,[real(sdm) imag(sdm)])
    #close(io)

    writedlm(fout, [real(sdm) imag(sdm)])
end

function read_cmatrix(fin)
    #io = open(fin,"r")
    #r=read(io,Int64)
    #c=read(io,Int64)
    #C = Array{Float64}(undef,r,c)

    #ld = read!(io,C)

    ld = readdlm(fin)
    Nh = div(size(ld,2),2)
    ld[:,1:Nh] + 1im .* ld[:,(Nh+1):end]
end

function enlarge_with_zeros!(SDM_new,SDM_old,thresholds_mask)
    size(SDM_new,1) != length(thresholds_mask) && error("Check size of SDM_new")
    size(SDM_old,1) != sum(thresholds_mask) && error("Check size of SDM_old")
    SDM_new[thresholds_mask,thresholds_mask] .= SDM_old
    return
end

function enlarge_with_zeros!(SDM_old,thresholds_mask)
    size(SDM_old,1) != sum(thresholds_mask) && error("Check size of SDM_old")
    SDM_new = fill(0.0im, length(thresholds_mask), length(thresholds_mask))
    SDM_new[thresholds_mask,thresholds_mask] .= SDM_old
    return SDM_new
end

end


text:
using HTTP, JSON, Sockets, Dates

include("routes.jl")

struct TimeoutError <: Exception end

function create_channel(path, payload)
  # spawn=false for PyCall to avoid segfault
  spawn = path != "/gmeans" && path != "/gaussian_mixtures"
  Channel(spawn=spawn) do ch
    try
      answer = routes(path, payload)
      isopen(ch) && put!(ch, HTTP.Response(answer))
    catch e
      println("500 $e")
      Base.showerror(Base.stdout, e, catch_backtrace())
      isopen(ch) && put!(ch, HTTP.Response(500, "Error: $e"))
    end
  end
end

function run(host, port)
  println("Start to server: listen to $host:$port")
  HTTP.serve(host, port, reuseaddr=true, rate_limit=100000//1) do request::HTTP.Request
    path = request.target
    current_time = Dates.now()
    println("[$current_time] Request $path")
    payload = Nothing
    try
      payload = JSON.parse(IOBuffer(HTTP.payload(request)))
    catch e
      println("500 Internal Server Error: $e")
      Base.showerror(Base.stdout, e, catch_backtrace())
      return HTTP.Response(500, "Error: $e")
    end

    channel = create_channel(path, payload)
    println("running...")
    timedoutSeconds = if haskey(ENV, "MIANALYZER_TIMEOUT"); parse(Float64, ENV["MIANALYZER_TIMEOUT"]); else 5.0; end
    timedwait(() -> isready(channel), timedoutSeconds)
    println("finished")

    if !isready(channel)
      isopen(channel) && close(channel)
      println("408 Timeout")
      return HTTP.Response(408, "Error: Request Timeout")
    end
    result = take!(channel)
    isopen(channel) && close(channel)
    return result
  end
end

function main()
  host = "0.0.0.0"
  port = 8000
  run(host, port)
end; main()


text:
"GARCH LLH benchmark, Julia version."

import JSON
using Test
using Statistics

data = JSON.parsefile("data.json")
out = Dict()

function garchSim(ɛ²::Vector, ω, α, β)
    h = similar(ɛ²)
    h[1] = mean(ɛ²)
    for i = 2:length(ɛ²)
        @inbounds h[i] = ω + α*ɛ²[i-1] + β*h[i-1]
    end
    h
end

function garchLLH(y::Vector, params::Vector)
    ɛ² = y.^2
    T = length(y)
    h = garchSim(ɛ², params...)
    -(T-1)*log(2π)/2 - sum(log.(h) + (y./sqrt.(h)).^2)/2
end

function benchmark(n)
    #price = Array{Float32}(data["price"])
    ret = Array{Float32}(data["ret"])
    #ret = diff(log.(price))
    #ret = ret .- mean(ret)
    x0 = Array{Float32}(data["x0"])
    tmp = garchLLH(ret, x0)  # warm-up
    @test tmp ≈ data["llh"] atol=1e-3
    t = @elapsed for i = 1:n; tmp = garchLLH(ret, x0); end
end

n_ix = findfirst(ARGS.=="-n")
if n_ix < length(ARGS)
    n = tryparse(Int64, ARGS[n_ix+1])
    n==nothing && error("Wrong number of iterations!")
else
    error("Number of iterations is not set!")
end

out["julia"] = benchmark(n)
println(JSON.json(out))

max_stars_repo_name: AndreyKolev/MLLibMark

max_issues_repo_name: AndreyKolev/MLLibMark

max_forks_repo_name: AndreyKolev/MLLibMark


text:
using OpenQASM
using OpenQASM.Types
using OpenQASM.Tools
using MLStyle
using RBNF: Token
using Test

@testset "cmp_exp" begin
    @test cmp_exp(Neg(qasm_f64(0.2)), qasm_f64(-0.2))
    @test cmp_exp(qasm_f64(-0.2), Neg(qasm_f64(0.2)))
end

@testset "conversion" begin
    @test convert(String, Token{:str}("\"abc\"")) == "abc"
    @test convert(Symbol, Token{:id}("abc")) === :abc
    @test convert(Float64, Token{:float64}("0.23")) == 0.23
    @test convert(Int, Token{:int}("1")) == 1
end

@testset "qasm token helper" begin
    @test qasm_id("abc") ≈ Token{:id}("abc")
    @test qasm_id(:abc) ≈ Token{:id}("abc")
    @test qasm_int(2) ≈ Token{:int}("2")
    @test qasm_f64(2.3) ≈ Token{:float64}("2.3")
    @test qasm_str("abc") ≈ Token{:str}("\"abc\"")
end

@testset "Bit(::String[, addrs])" begin
    bit = Bit("qreg")
    @test bit.name == qasm_id("qreg")
    @test bit.address === nothing

    bit = Bit("qreg", 2)
    @test bit.name == qasm_id("qreg")
    @test bit.address == qasm_int(2)
end

@testset "qasm parser" begin

    qasm = """OPENQASM 2.0;
    include "qelib1.inc";
    gate custom(lambda) a {
        u1(sin(lambda) + 1) a;
    }
    // comment
    gate g a
    {
        U(0,0,0) a;
    }

    qreg q[4];
    creg c1[1];
    creg c2[1];
    U(-1.0, pi/2+3, 3.0) q[2];
    CX q[1], q[2];
    custom(0.3) q[3];
    barrier q;
    h q[0];
    measure q[0] -> c0[0];
    if(c0==1) z q[2];
    u3(0.1 + 0.2, 0.2, 0.3) q[0];
    """

    ast = OpenQASM.parse(qasm)
    println(ast)

    @testset "mainprogram" begin
        @test ast isa MainProgram
        @test ast.version == v"2.0.0"
    end

    @test @match ast begin
        MainProgram(v"2.0.0", progs) => true
        _ => false
    end

    @testset "include" begin
        @test ast.prog[1] isa Include
        @test ast.prog[1].file isa Token{:str}
        @test ast.prog[1].file.str == "\"qelib1.inc\""

        token = @match ast.prog[1] begin
            Include(file) => file
            _ => nothing
        end
        @test token.str == "\"qelib1.inc\""
    end

    @testset "gate custom(lambda) a" begin
        @test ast.prog[2] isa Gate
        custom = ast.prog[2]

        @test @match custom begin
            Gate(GateDecl(name, cargs, qargs), body) => true
            _ => false
        end

        @test custom.decl.name isa Token{:id}
        @test custom.decl.name.str == "custom"
        @test length(custom.decl.cargs) == 1
        @test custom.decl.cargs[1] isa Token{:id}
        @test custom.decl.cargs[1].str == "lambda"
        @test length(custom.decl.qargs) == 1
        @test custom.decl.qargs[1] isa Token{:id}
        @test custom.decl.qargs[1].str == "a"

        @test length(custom.body) == 1
        @test custom.body[1] isa Instruction
        inst = custom.body[1]
        @test inst.name == "u1"
        @test length(inst.cargs) == 1
        @test length(inst.qargs) == 1
        @test inst.cargs[1] isa Tuple
        @test inst.qargs[1] isa Bit
        carg = inst.cargs[1]
        qarg = inst.qargs[1]
        @test length(carg) == 3
        @test carg[1] isa Call
        @test carg[2] isa Token{:reserved}
        @test carg[3] isa Token{:int}
        @test carg[1].name === :sin
        @test carg[1].args.str == "lambda"
        @test carg[2].str == "+"
        @test carg[3].str == "1"
        @test qarg.name.str == "a"
        @test qarg.address === nothing
    end

    @testset "gate g" begin
        @test ast.prog[3] isa Gate
        g = ast.prog[3]
        @test isempty(g.decl.cargs)
        @test length(g.decl.qargs) == 1
        @test g.decl.qargs[1] isa Token{:id}
        @test g.decl.qargs[1].str == "a"
        @test g.decl.name isa Token{:id}
        @test g.decl.name.str == "g"

        @test length(g.body) == 1
        @test g.body[1] isa UGate
        ugate = g.body[1]
        @test ugate.qarg isa Bit
        @test ugate.qarg.name.str == "a"
        @test ugate.qarg.address === nothing
        @test ugate.y isa Token{:int}
        @test ugate.y.str == "0"
        @test ugate.z1.str == "0"
        @test ugate.z2.str == "0"
    end

    @testset "qreg" begin
        @test ast.prog[4] isa RegDecl
        reg = ast.prog[4]
        @test reg.name isa Token{:id}
        @test reg.name.str == "q"
        @test reg.size isa Token{:int}
        @test reg.size.str == "4"
        @test reg.type isa Token{:reserved}
        @test reg.type.str == "qreg"
    end

    @testset "creg" begin
        @test ast.prog[5] isa RegDecl
        reg = ast.prog[5]
        @test reg.name isa Token{:id}
        @test reg.name.str == "c1"
        @test reg.size isa Token{:int}
        @test reg.size.str == "1"
        @test reg.type isa Token{:reserved}
        @test reg.type.str == "creg"
    end

    @testset "U" begin
        @test ast.prog[7] isa UGate
        U = ast.prog[7]
        @test U.qarg isa Bit
        @test U.y isa Tuple
        @test U.z1 isa Neg
        @test U.z2 isa Token{:float64}
        @test U.qarg.name.str == "q"
        @test U.qarg.address.str == "2"
        @test U.y[2] isa Token{:reserved}
        @test U.y[2].str == "+"
        @test U.y[1] isa Tuple
        @test U.y[1][1] isa Token{:reserved}
        @test U.y[1][2] isa Token{:reserved}
        @test U.y[1][1].str == "pi"
        @test U.y[1][2].str == "/"
    end

    @testset "CX" begin
        @test ast.prog[8] isa CXGate
        cx = ast.prog[8]
        @test cx.ctrl isa Bit
        @test cx.qarg isa Bit
        @test cx.ctrl.name.str == "q"
        @test cx.qarg.name.str == "q"
        @test cx.ctrl.address.str == "1"
        @test cx.qarg.address.str == "2"
    end

    @testset "inst" begin
        @test ast.prog[9] isa Instruction
        inst = ast.prog[9]
        @test inst.name == "custom"
        @test inst.cargs[1].str == "0.3"
        @test inst.qargs[1].name.str == "q"
        @test inst.qargs[1].address.str == "3"
    end

    @testset "barrier" begin
        @test ast.prog[10] isa Barrier
        barrier = ast.prog[10]
        @test length(barrier.qargs) == 1
        @test barrier.qargs[1] isa Bit
        @test barrier.qargs[1].name.str == "q"
        @test barrier.qargs[1].address === nothing
    end

    @testset "h q[0]" begin
        @test ast.prog[11] isa Instruction
        @test ast.prog[11].name == "h"
        @test length(ast.prog[11].qargs) == 1
        @test length(ast.prog[11].cargs) == 0
    end

    @testset "measure" begin
        m = ast.prog[12]
        @test m isa Measure
        @test m.qarg isa Bit
        @test m.carg isa Bit
        @test m.qarg.name isa Token{:id}
        @test m.qarg.name.str == "q"
        @test m.carg.name isa Token{:id}
        @test m.carg.name.str == "c0"
        @test m.qarg.address isa Token{:int}
        @test m.carg.address isa Token{:int}
        @test m.qarg.address.str == "0"
        @test m.carg.address.str == "0"
    end

    @testset "ifstmt" begin
        ifstmt = ast.prog[13]
        @test ifstmt isa IfStmt
        @test ifstmt.left isa Token{:id}
        @test ifstmt.left.str == "c0"
        @test ifstmt.right isa Token{:int}
        @test ifstmt.right.str == "1"
        @test ifstmt.body isa Instruction
        @test ifstmt.body.name == "z"
        @test length(ifstmt.body.qargs) == 1
        @test length(ifstmt.body.cargs) == 0
        @test ifstmt.body.qargs[1].name.str == "q"
        @test ifstmt.body.qargs[1].address.str == "2"
    end

    @testset "u3" begin
        inst = ast.prog[14]
        @test inst isa Instruction
        @test inst.name == "u3"
        @test length(inst.cargs) == 3
        @test inst.cargs[1][1].str == "0.1"
        @test inst.cargs[1][2].str == "+"
        @test inst.cargs[1][3].str == "0.2"
        @test inst.cargs[2].str == "0.2"
        @test inst.cargs[3].str == "0.3"
    end

end

@testset "cmp_ast" begin
    qasm1 = """OPENQASM 2.0;
    include "qelib1.inc";
    gate custom(lambda) a {
        u1(sin(lambda) + 1) a;
    }
    // comment
    gate g a
    {
        U(0,0,0) a;
    }

    qreg q[4];
    creg c1[1];
    creg c2[1];
    U(-1.0, pi/2+3, 3.0) q[2];
    CX q[1], q[2];
    custom(0.3) q[3];
    barrier q;
    h q[0];
    measure q[0] -> c0[0];
    if(c0==1) z q[2];
    u3(0.1 + 0.2, 0.2, 0.3) q[0];
    """

    qasm2 = """OPENQASM 2.0;
    include "qelib1.inc";
    gate custom(lambda) a {
        u1(sin(lambda) + 1) a;
    }
    // comment
    gate g a
    {
        U(0,0,0) a;
    }

    qreg q[4];
    creg c1[1];
    creg c2[1];
    U(-1.0, pi/2+3, 3.0) q[2];
    CX q[1], q[2];
    custom(0.3) q[3];
    barrier q;
    h q[0];
    measure q[0] -> c0[0];
    if(c0==1) z q[2];
    """

    ast1 = OpenQASM.parse(qasm1)
    ast2 = OpenQASM.parse(qasm2)
    println(ast1)
    println(ast2)

    @test !(ast1 ≈ ast2)
    @test ast1 ≈ ast1
    @test ast2 ≈ ast2

    s = """
    gate test_gate(theta, phi) qreg_2, qreg_3, qreg_1 {
    x qreg_1;
    z qreg_2;
    rx(-(sin(theta))+2.0) qreg_3;
    rx(sin(theta)+tan(phi)) qreg_3;
    ry(cos(phi)-ln(phi)) qreg_3;
    ry(cos(phi)*sqrt(phi)) qreg_3;
    CX qreg_1, qreg_3;
    }
    """
    ast3 = OpenQASM.parse_gate(s)
    println(ast3)
    @test ast3 ≈ ast3
end


text:
import Dates
using Dates

"""
    mutable struct DoublesMatchPlayerSummary
        opponent1_player_id::Int64
        opponent2_player_id::Int64
        result::String
        score::String
    end

An object containing information regarding a player's match results within a tournament.

Fields
=========
- `opponent1_player_id`:
    the opposing team's player1 id
- `opponent2_player_id`:
    the opposing team's player2 id
- `result`:
    whether the player won or lost the match.
    options:
        1) "Won"
        2) "Lost"
- `score`
    the score of the match
"""

mutable struct DoublesMatchPlayerSummary
    opponent1_player_id::Int64
    opponent2_player_id::Int64
    result::String
    score::String
    function DoublesMatchPlayerSummary(;
        opponent1_player_id,
        opponent2_player_id,
        result,
        score,
    )
        return new(opponent1_player_id, opponent2_player_id, result, score)
    end
end

"""
    mutable struct DoublesTournamentPlayerSummary
        primary_player_id::Int64
        partner_player_id::Int64
        tournament_name::String
        tournament_type::String
        result::Int128
        matches::Vector{DoublesMatchPlayerSummary}
    end

An object containing information regarding a player's tournament results.

Fields
=========
- `primary_player_id`:
    the player's id.
- `partner_player_id`:
    the player's partner's id.
- `tournament_name`:
    the name of the tournament being played.
- `tournament_type`:
    the type of the tournament being played.
    options:
        1) Grandslam
        2) ATPNitto
        3) ATP1000
        4) ATP500
        5) ATP250
        6) ATPChallenger125
        7) ATPChallenger110
        8) ATPChallenger100
        9) ATPChallenger90
        10) ATPChallenger80
        11) ITFWorldTennisTour25k
        12) ITFWorldTennisTour15k
- `result`
    the finishing spot of the player in the tournament.
    options:
        1) "Winner"
        2) "Finals"
        3) "Semi-Finals"
        4) "Quarter-Finals"
        5) "Round-16"
        6) "Round-32"
        7) "Round-64"
        8) "Round-128"
        9) "Qualified"
- `matches`
    an array of information for matches played in the tournament.
"""

mutable struct DoublesTournamentPlayerSummary
    primary_player_id::Int64
    partner_player_id::Int64
    tournament_name::String
    tournament_type::String
    result::String
    matches::Vector{DoublesMatchPlayerSummary}
    function DoublesTournamentPlayerSummary(;
        primary_player_id,
        partner_player_id,
        tournament_name,
        tournament_type,
        result,
        matches,
    )
        return new(
            primary_player_id,
            partner_player_id,
            tournament_name,
            tournament_type,
            result,
            matches,
        )
    end
end

"""
    mutable struct Player
        id::Int64
        rankings::Dict{String, Int}
        matches_played::Int128
        tournaments_played::Int128
        doubles_results::Dict{String, DoublesTournamentPlayerSummary}
    end

A tennis player.

Fields
=========
- `id`:
    a unique identifier for the player.
- `rankings`:
    contains all latest rankings for the player
        ex:
            Dict(
                "Ground Truth" => 1,
                "Burn In" => 38,
                "ATP" => 2,
                "Elo" => 17
            )
- `matches_played`:
    a count of the number of matches a player has played.
- `tournaments_played`:
    a count of the number of tournaments a player has played.
- `doubles_results`:
    the results of tournaments a player has played in.
        ex:
            Dict(
                "Wimbledon_2018-07-02" => DoublesTournamentPlayerSummary(
                    primary_player_id = 1,
                    partner_player_id = 2,
                    tournament_name = "Wimbledon",
                    tournament_type = "Grandslam",
                    result = "Winner",
                    matches = Vector{DoublesMatchPlayerSummary}([])
                )
            )
"""

mutable struct Player
    id::Int64
    age::Int64
    rankings::Dict{String,Float64}
    matches_played::Int128
    tournaments_played::Int128
    doubles_results::Dict{String,DoublesTournamentPlayerSummary}
    function Player(;
        id::Int64,
        ground_truth_ranking::Int64,
        burn_in_ranking::Int64,
        ranking_schemes::Vector{String},
    )
        id = id
        age = -1
        matches_played = 0
        tournaments_played = 0
        doubles_results = Dict{String,DoublesTournamentPlayerSummary}()
        rankings::Dict{String,Int} = Dict(
            "Ground Truth" => ground_truth_ranking,
            "Burn In" => burn_in_ranking,
        )
        for rankingScheme in ranking_schemes
            rankings[rankingScheme] = 0
        end
        return new(
            id,
            age,
            rankings,
            matches_played,
            tournaments_played,
            doubles_results,
        )
    end
end

"""
    mutable struct Team
        player1::Player
        player2::Player
        ranking::Int64
    end

An object containing two players to form a doubles pair.

Fields
=========
- `player1`:
    the first player in the team.
- `opponent2_player_id`:
    the second player in the team.
- `ranking`:
    the overall ranking of the team.
"""

mutable struct Team
    player1::Player
    player2::Player
    ranking::Float64
    function Team(player1::Player, player2::Player, ranking_scheme::String)
        player1 = player1
        player2 = player2
        p1_rank = get(player1.rankings, ranking_scheme, -1)
        p2_rank = get(player2.rankings, ranking_scheme, -1)
        (p1_rank != -1 && p2_rank != -1) ||
            throw(DomainError(ranking_scheme, "ranking scheme unavailable"))
        ranking = p1_rank + p2_rank
        return new(player1, player2, ranking)
    end
end

"""
    mutable struct Tournament
        id::String
        name::String
        date::Date
        type::String
        drawsize::Int64
    end

A tournament played during the tennis season.

Fields
=========
- `id`:
    a unique identifier for the tournament combining name and date
- `name`:
    the name of the tournament being played
- `date`:
    the date when the tournament was played
- `type`
    the type of the tournament
    options:
        1) Grandslam
        2) ATPNitto
        3) ATP1000
        4) ATP500
        5) ATP250
        6) ATPChallenger125
        7) ATPChallenger110
        8) ATPChallenger100
        9) ATPChallenger90
        10) ATPChallenger80
        11) ITFWorldTennisTour25k
        12) ITFWorldTennisTour15k
- `drawsize`
    the draw size for the tournament
"""

mutable struct Tournament
    id::String
    name::String
    date::Date
    type::String
    drawsize::Int64
    function Tournament(;
        name::String,
        date::Date,
        type::String,
        drawsize::Int64,
    )
        id = string(name, "_", Dates.format(date, "yyyy-mm-dd"))
        return new(id, name, date, type, drawsize)
    end
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
SST_PATH = "/Users/ccolley/Documents/Research/SparseSymmetricTensors.jl/src/SparseSymmetricTensors.jl" # local path
#SST_PATH = "/homes/ccolley/Documents/Software/SparseSymmetricTensors.jl/src/SparseSymmetricTensors.jl" #Nilpotent path
include(SST_PATH)
using Main.SparseSymmetricTensors

#===============================================================================
                            Contraction Codes
===============================================================================#

#TODO: adapt to work for any order tensor
function implicit_contraction(A::COOTen,B::COOTen,x::Array{Float64,1})

    @assert length(x) == A.cubical_dimension*B.cubical_dimension
    m = A.cubical_dimension
    n = B.cubical_dimension
    y = similar(x)

	y .= 0

    ileave = (i,j) -> i + m*(j-1)

    for i in 1:length(A)

        for (i_1,i_2,i_3) in permutations(A.indices[i,:])

            for j in 1:length(B)
                j_1,j_2,j_3 = B.indices[j,:]

                y[ileave(i_1,j_1)] += 2*A.vals[i]*B.vals[j]*x[ileave(i_2,j_2)]*x[ileave(i_3,j_3)]
                y[ileave(i_2,j_2)] += 2*A.vals[i]*B.vals[j]*x[ileave(i_1,j_1)]*x[ileave(i_3,j_3)]
                y[ileave(i_3,j_3)] += 2*A.vals[i]*B.vals[j]*x[ileave(i_1,j_1)]*x[ileave(i_2,j_2)]

            end
        end
    end

    return y
end

#===============================================================================
                            Matching Code
===============================================================================#

#Computes the TAME score for this iterate by
function TAME_score(A::COOTen,B::COOTen,u::Array{Float64,1},v::Array{Float64,1})

    Match_mapping, weight = rank_one_matching(u,v)
    TAME_score(A,B,Match_mapping)

end

function TAME_score(A::COOTen,B::COOTen,U::Array{Float64,2},V::Array{Float64,2})

    Match_mapping = low_rank_matching(U,V)
    TAME_score(A,B,Match_mapping)

end

function TAME_score(A::COOTen,B::COOTen,X::Array{Float64,2};return_timings=false)

    if return_timings
         (_,_,matching,_) ,scoring_time = @timed bipartite_matching_primal_dual(X)
         (triangle_count, gaped_triangles,inverted_matching), matching_time = @timed TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
         return triangle_count, gaped_triangles,inverted_matching, matching_time, matching_time
     else
         (_,_,matching,_) ,scoring_time = @timed bipartite_matching_primal_dual(X)
         return TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
     end

 end

 function TAME_score(A::COOTen,B::COOTen,x::Array{Float64,1};return_timings=false)

    if return_timings
         (_,_,matching) ,matching_time = @timed bipartite_matching_primal_dual(reshape(x,A.cubical_dimension,B.cubical_dimension))
         (triangle_count, gaped_triangles,inverted_matching), scoring_time = @timed TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
         return triangle_count, gaped_triangles,inverted_matching, matching_time,scoring_time
     else
         _,_,matching = bipartite_matching_primal_dual(reshape(x,A.cubical_dimension,B.cubical_dimension))
         return TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
     end

 end

function TAME_score(A::COOTen,B::COOTen,Match_mapping::Dict{Int,Int})

    match_len = length(Match_mapping)

    Triangle_check = Dict{Array{Int,1},Int}()
    gaped_triangles = 0
    triangle_count = 0

    if A.unique_nnz > B.unique_nnz

        for i in 1:A.unique_nnz
            Triangle_check[A.indices[i,:]] = 1
        end

        #invert to map v indices to u
        Match_mapping = Dict(value => key for (key, value) in Match_mapping)

        for i in 1:B.unique_nnz
            v_i,v_j,v_k = B.indices[i,:]

            matched_triangle =
              sort([get(Match_mapping,v_i,-1),get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])
    #        println(B.indices[i,:]," -> ",matched_triangle)
            match = get(Triangle_check,matched_triangle,0)
            if match == 1
                triangle_count += 1
            else
                gaped_triangles += 1
            end
        end

    else
        for i in 1:B.unique_nnz
            Triangle_check[B.indices[i,:]] = 1
        end

        for i in 1:A.unique_nnz
            v_i,v_j,v_k = A.indices[i,:]
            matched_triangle =
               sort([get(Match_mapping,v_i,-1), get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])

            match = get(Triangle_check,matched_triangle,0)
            if match == 1
                triangle_count += 1
            else
                gaped_triangles += 1
            end
        end
    end

   # triangles, triangle_weight = count_triangles(sub_A,sub_B)

    return triangle_count, gaped_triangles, Match_mapping

end

function TAME_score(Triangle_Dict::Dict{Array{Int,1},Int},Input_tensor::COOTen,
    Match_mapping::Dict{Int,Int})

    triangle_count = 0
    gaped_triangles = 0

    for i in 1:Input_tensor.unique_nnz
        v_i,v_j,v_k = Input_tensor.indices[i,:]

        matched_triangle =
            sort([get(Match_mapping,v_i,-1),get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])

        match = get(Triangle_Dict,matched_triangle,0)
        if match == 1
            triangle_count += 1
        else
            gaped_triangles += 1
        end
    end
    return triangle_count, gaped_triangles, Match_mapping
end

function search_Krylov_space(A::COOTen,B::COOTen,U::Array{Float64,2},V::Array{Float64,2})

    best_score = -1
    best_i = -1
    best_j = -1
    best_matching = Dict{Int,Int}()

    Triangle_check = Dict{Array{Int,1},Int}()

    if A.unique_nnz > B.unique_nnz
        for i in 1:A.unique_nnz
            Triangle_check[A.indices[i,:]] = 1
        end
        Input_tensor = B
    else
        for i in 1:B.unique_nnz
            Triangle_check[B.indices[i,:]] = 1
        end
        Input_tensor = A
    end

    for i in 1:size(U,2)
       for j in 1:size(V,2)

            if A.unique_nnz > B.unique_nnz
                matched_tris, gaped_tris, matching = TAME_score(Triangle_check,Input_tensor,V[:,j],U[:,i])
            else
                matched_tris, gaped_tris, matching = TAME_score(Triangle_check,Input_tensor,U[:,i],V[:,j])
            end

            if matched_tris > best_score
                best_matching = matching
                best_score = matched_tris
                best_i = i
                best_j = j
            end
        end
    end
    return best_score, best_i, best_j, best_matching
end

#===============================================================================
                          TAME Implementations Code
===============================================================================#

#
#  Parameter Searching Code
#

function align_tensors(A::COOTen,B::COOTen;method::String="LambdaTAME",
                       no_matching=false,kwargs...)

    #put larger tensor on the left
    if B.n > A.n
        results = align_tensors(B,A;method = method, no_matching=no_matching,kwargs...)
        #flip the matchings if A and B were swapped
        if method == "LambdaTAME" ||  method == "LowRankTAME"
            best_TAME_PP_tris, max_triangle_match, U_best, V_best, best_matching = results
            return best_TAME_PP_tris, max_triangle_match, U_best, V_best, Dict((j,i) for (i,j) in best_matching)
        elseif method == "TAME"
            best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching = results
            return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, Dict((j,i) for (i,j) in best_matching)
        end
    end

    if method == "LambdaTAME"
        return ΛTAME_param_search(A,B;kwargs...)
    elseif method == "LowRankTAME"
        return LowRankTAME_param_search(A,B;no_matching = no_matching,kwargs...)
    elseif method == "TAME"
        return TAME_param_search(A,B;no_matching = no_matching,kwargs...)
    else
        throw(ArgumentError("method must be one of 'LambdaTAME', 'LowRankTAME',or 'TAME'."))
    end
end

function align_tensors_profiled(A::COOTen,B::COOTen; method::String="LambdaTAME",
                                no_matching=false,kwargs...)

    #put larger tensor on the left
    if B.n > A.n
        results =  align_tensors_profiled(B,A;method = method, no_matching=no_matching,kwargs...)
        #flip the matchings if A and B were swapped
        if method == "LambdaTAME" ||  method == "LowRankTAME"
            best_TAME_PP_tris, max_triangle_match, U_best, V_best, best_matching,profile = results
            return best_TAME_PP_tris, max_triangle_match, U_best, V_best, Dict((j,i) for (i,j) in best_matching), profile
        elseif method == "TAME"
            best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching,profile = results
            return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, Dict((j,i) for (i,j) in best_matching), profile
        end
    end

    if method == "LambdaTAME"
        return ΛTAME_param_search_profiled(A,B;kwargs...)
    elseif method == "LowRankTAME"
        return LowRankTAME_param_search_profiled(A,B;no_matching = no_matching,kwargs...)
    elseif method == "TAME"
        return TAME_param_search_profiled(A,B;no_matching = no_matching,kwargs...)
    else
        throw(ArgumentError("method must be one of 'LambdaTAME', 'LowRankTAME', or 'TAME'."))
    end
end

function ΛTAME_param_search_profiled(A::COOTen,B::COOTen; iter::Int = 15,tol::Float64=1e-6,
                                     alphas::Array{F,1}=[.5,1.0],
                                     betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001]) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    best_TAME_PP_tris = -1
    best_i  = -1
    best_j = -1
    best_matching = Dict{Int,Int}()

    if Ten == COOTen
        m = A.cubical_dimension
        n = B.cubical_dimension
    else
        m = A.n
        n = B.n
    end

    results = Dict(
        "TAME_timings" => Array{Float64,1}(undef,length(alphas)*length(betas)),
        "Krylov Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas))
    )
    exp_index = 1

    U = Array{Float64,2}(undef,m,iter)
    V = Array{Float64,2}(undef,n,iter)

    for α in alphas
        for beta in betas

            ((U,V),runtime) = @timed ΛTAME(A,B,beta,iter,tol,α)
            results["TAME_timings"][exp_index] = runtime

            #search the Krylov Subspace
            ((search_tris, i, j, matching),runtime) = @timed search_Krylov_space(A,B,U,V)
            results["Krylov Timings"][exp_index] = runtime
            exp_index += 1

            if search_tris > best_TAME_PP_tris
                best_matching = matching
                best_TAME_PP_tris = search_tris
                best_i = i
                best_j = j
            end

            println("α:$(α) -- β:$(beta) finished -- tri_match:$search_tris -- max_tris $(max_triangle_match) -- best tri_match: $best_TAME_PP_tris")
        end
    end

    println("best i:$best_i -- best j:$best_j")
    return best_TAME_PP_tris, max_triangle_match, U[:,best_i], V[:,best_j], best_matching, results

end

#add in SparseSymmetricTensors.jl function definitions
function TAME_param_search(A::COOTen,B::COOTen;iter::Int = 15,tol::Float64=1e-6,
                           alphas::Array{F,1}=[.5,1.0],
                           betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
                           kwargs...) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
    best_TAME_PP_tris::Int = -1
    best_matching = Dict{Int,Int}()

    m = A.cubical_dimension
    n = B.cubical_dimension

    best_TAME_PP_x = Array{Float64,2}(undef,m,n)

    for α in alphas
        for β in betas

            x, triangle_count, matching = TAME(A,B,β,iter,tol,α;W = ones(m,n),kwargs...)

            if triangle_count > best_TAME_PP_tris
                best_matching = matching
                best_TAME_PP_tris = triangle_count
                best_TAME_PP_x = copy(x)
            end
            println("α:$(α) -- β:$β finished -- tri_match:$(triangle_count) -- max_tris $(max_triangle_match) -- best tri_match:$(best_TAME_PP_tris)")
        end

    end

    return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching
end

function TAME_param_search_profiled(A::COOTen,B::COOTen;iter::Int = 15,tol::Float64=1e-6,
                                    alphas::Array{F,1}=[.5,1.0],
                                    betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
                                    profile::Bool=false,profile_aggregation="all",
                                    kwargs...) where {F <: AbstractFloat}
    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
    best_TAME_PP_tris = -1
    best_matching = Dict{Int,Int}()

    m = A.cubical_dimension
    n = B.cubical_dimension

    best_TAME_PP_x = Array{Float64,2}(undef,m,n)
    experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

            x, triangle_count, matching, experiment_profile = TAME_profiled(A,B,β,iter,tol,α;W = ones(m,n),kwargs...)
            push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

            if triangle_count > best_TAME_PP_tris
                best_matching = matching
                best_TAME_PP_tris = triangle_count
                best_TAME_PP_x = copy(x)
            end
            println("α:$(α) -- β:$β finished -- tri_match:$(best_TAME_PP_tris) -- max_tris $(max_triangle_match)")
        end
    end

    return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching, experiment_profiles

end

function LowRankTAME_param_search(A::COOTen,B::COOTen;iter::Int = 15,tol::Float64=1e-6,
                                  U_0::Array{Float64,2} = ones(A.n,1), V_0::Array{Float64,2} = ones(B.n,1),
                                  alphas::Array{F,1}=[.5,1.0],
                                  betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
                                  kwargs...) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
    best_TAME_PP_tris = -1
    best_matching = Dict{Int,Int}()

    m = A.cubical_dimension
    n = B.cubical_dimension

    best_TAME_PP_U = ones(m,1)
    best_TAME_PP_V = ones(n,1)

    for α in alphas
        for β in betas
            U, V, triangle_count,matching =
            LowRankTAME(A,B,U_0,V_0,β,iter,tol,α;kwargs...)

            if triangle_count > best_TAME_PP_tris
                best_TAME_PP_tris = triangle_count
                best_matching = matching
                best_TAME_PP_U = copy(U)
                best_TAME_PP_V = copy(V)
            end
            println("α:$(α) -- β:$(β) -- tri_match:$(triangle_count) -- max_tris:$(max_triangle_match) -- best tri match:$best_TAME_PP_tris")
        end
    end

    return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_U, best_TAME_PP_V, best_matching

end

#
#  Spectral Relaxation Code
#

function ΛTAME(A::COOTen, B::COOTen, β::Float64, max_iter::Int,
    tol::Float64,α::Float64;update_user::Int=-1)

    U = zeros(A.cubical_dimension,max_iter+1)
    V = zeros(B.cubical_dimension,max_iter+1) #store initial in first column

    U[:,1] = ones(A.cubical_dimension)
    U[:,1] /=norm(U[:,1])

    V[:,1] = ones(B.cubical_dimension)
    V[:,1] /=norm(U[:,1])

    sqrt_β = β^(.5)

    lambda = Inf
    i = 1

    while true

        U[:,i+1] = contract_k_1(A,U[:,i])
        V[:,i+1] = contract_k_1(B,V[:,i])

        lambda_A = (U[:,i+1]'*U[:,i])
        lambda_B = (V[:,i+1]'*V[:,i])
        new_lambda = lambda_A*lambda_B

        if β != 0.0
            U[:,i+1] .+= sqrt_β*U[:,i+1]
            V[:,i+1] .+= sqrt_β*V[:,i+1]
        end

        if α != 1.0
            U[:,i+1] = α*U[:,i+1] + (1 -α)*U[:,1]
            V[:,i+1] = α*V[:,i+1] + (1 -α)*V[:,1]
        end

        U[:,i+1] ./= norm(U[:,i+1])
        V[:,i+1] ./= norm(V[:,i+1])

        if update_user != -1 && i % update_user == 0
            println("iteration $(i)    λ_A: $(lambda_A) -- λ_B: $(lambda_B) -- newλ: $(new_lambda)")
        end

        if abs(new_lambda - lambda) < tol || i >= max_iter
            return U[:,1:i], V[:,1:i]
        else
            lambda = new_lambda
            i += 1
        end
    end
end

#=
  TODO: currently missing COOTen Implementations for
   - LowRankTAME
   - LowRankTAME_profiled
   - TAME
   - TAME_profiled
=#

#===============================================================================
                            PostProcessing Code
===============================================================================#
"""-----------------------------------------------------------------------------
   Computes a post processing routine on the optimal aligment produced by
   Λ-TAME. The algorithm runs a b-matching to find suitable replacements
   and computes whether or not to make the swap if the number of aligned
   triangles increases.

   Inputs:
   -------
   {A,B} - (COOTen)
     Coordinate representations of adjacency tensors for the motif induced
     hypergraphs.
   {u,v} - (Array{Float64,1})
     The components to the rank-1 alignment matrix which offered the best
     alignment scores.
   iterations - (Int)
     The number of iterations to run the post-processing algorithm for.
   b - (Int)
     The size of the b-matching.

-----------------------------------------------------------------------------"""
function post_process(A::COOTen,B::COOTen,u::Array{Float64,1},
                      v::Array{Float64,1},iterations::Int,b::Int)

    @assert length(u) <= length(v)

    Matching = rank_one_matching(u,v)

    A_incidence = produce_triangle_incidence(A.indices,A.cubical_dimension)
    B_incidence = produce_triangle_incidence(B.indices,B.cubical_dimension)

    potential_matchings = b_matching(u,v,b)

    u_sum = sum(u)
    v_sum = sum(v)

    match_q = Queue{Tuple{Int,Int,Float64}}()

    for _ in 1:iterations
        #sort by ∑_i x(ii') + ∑_i' x(ii')
        [enqueue!(q,z) for z in sort(potential_matchings,by=x->u[x[1]]*u_sum+ v[x[2]]*v_sum)]
        while length(q) != 0

            i,ip,_ = dequeue!(match_q)

            matched_tris = compute_triangles_aligned(A_incidence[i],B_incidence[ip],Matching)

            #TODO: optimize this later
            Pref_i = [k for (j,k,_) in potential_matchings if j == i]
            Pref_ip = [k for (j,k,_) in potential_matchings if j == ip]

        end
    end
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
"""
Optimisation algorithm following

C. Chen, D. Carlson, Z. Gan, C. Li, and L. Carin. Bridging the gap between
stochastic gradient MCMC and stochastic optimization. AISTATS, 2016.

Fields:
  -x: postion
  -v: velocity
  -stepsize
  -alpha
  -u
  -sigma: geometric averaging factor for grad.^2
  -lambda
"""

## santa SGD
type SantaState <: SamplerState
    x::Array{Float64}
    v::Array{Float64}
    stepsize::Float64
    alpha::Array{Float64}
    u::Array{Float64}
    sigma::Float64
    lambda::Float64

    function SantaState(x::Array{Float64};
                        v=zeros(length(x)),stepsize = 4e-11,
                        alpha=1000*sqrt(stepsize)*ones(length(x)),
                        u=sqrt(stepsize)*randn(length(x)),
                        sigma=0.1,lambda=1e-5)
        new(x,v,stepsize,alpha,u,sigma,lambda)
    end
end

function sample!(s::SantaState, grad)
  sigma=s.sigma
  stepsize=s.stepsize

  f=sgrad(s.x)
  s.v=sigma*s.v+(1-sigma)*f.^2
  g=1./(sqrt(s.lambda+sqrt(s.v)))
  s.u[:] = (1-s.alpha).*s.u[:] + stepsize*g.*f*60000/2
  s.x[:] += g.*s.u
end


text:
ia> [RpcApi.SimpleOrganism(0x000000000000464e,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         Creature.stepUp(o)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,25952,0x0104b17d,Dict{Int16,Int16}(),[156,113],12),RpcApi.SimpleOrganism(0x0000000000004656,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,25373,0x0104e81a,Dict{Int16,Int16}(),[137,143],11),RpcApi.SimpleOrganism(0x000000000000455a,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         Creature.stepUp(o)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,16782,0x0104b17d,Dict{Int16,Int16}(),[125,88],12),RpcApi.SimpleOrganism(0x000000000000468c,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         Creature.stepUp(o)
         Creature.stepLeft(o)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,15459,0x0104e81a,Dict{Int16,Int16}(),[53,89],13),RpcApi.SimpleOrganism(0x0000000000004671,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         Creature.stepUp(o)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,14872,0x0104b17d,Dict{Int16,Int16}(),[40,150],12),RpcApi.SimpleOrganism(0x0000000000004685,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         Creature.stepUp(o)
         return var_2
     end
     local var_1::Int8 = 48
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,13442,0x0104b17d,Dict{Int16,Int16}(),[100,47],12),RpcApi.SimpleOrganism(0x00000000000045ed,:(function (o::Creature.Organism,)
     function func_4(var_2::Int64=8688889982731189829,var_3::Bool=true)
         return var_2
     end
     local var_1::Int8 = 48
     Creature.eatDown(o,Int(var_1))
     if var_1 == var_1
     end
     Creature.stepUp(o)
     Creature.stepRight(o)
     Creature.eatRight(o,Int(var_1))
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     Creature.stepLeft(o)
     var_1 = var_1 - var_1
     return true
 end),[100,300,95,0,1,1,1],1,5000,1,12429,0x0104e81a,Dict{Int16,Int16}(),[88,143],11)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# name = "Couenne_jll"
# Couenne_jll = "f09e9e23-9010-5c9e-b679-9f1d8f79b85c"
# name = "AmplNLWriter"
# AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"

using AmplNLWriter, Couenne_jll

COUENNE_OPT = AmplNLWriter.Optimizer(Couenne_jll.amplexe)
# MOI.set(COUENNE_OPT, MOI.Silent(), true)
# MOI.set(COUENNE_OPT, MOI.TimeLimitSec(), 5)
# COUENNE_OPT = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(COUENNE_OPT)

push!(solvers_nlp, (opt = COUENNE_OPT, mode = BilevelJuMP.ProductMode(1e-5)))


text:
## Performance tests
import IntSets

const REPS = 100000

function test_density{T}(::Type{T}, N)
    sz = 10000

    times = Float64[]
    param = Float64[]
    test = Symbol[]
    sizehint!(times, N*10*11)
    sizehint!(param, N*10*11)
    sizehint!(test, N*10*11)
    for density=logspace(-3,0,7)
        println("Density: $density")
        elts = unique(round.(Int, 1:1/density:sz))
        set = sizehint!(T(),sz)
        scratch1 = sizehint!(T(), sz)
        scratch2 = sizehint!(T(), sz)
        empty = sizehint!(T(), sz)
        full = T(1:sz)
        for _=1:N
            # push!
            push!(times, @elapsed for e in elts push!(set, e); end)
            push!(param, density)
            push!(test, :push!)

            # Isempty
            push!(times, @elapsed isempty(set))
            push!(param, density)
            push!(test, :isempty)

            # Iteration (sum)
            push!(times, @elapsed sum(set))
            push!(param, density)
            push!(test, :sum)

            # Containment
            push!(times, @elapsed for e=1:sz e in set end)
            push!(param, density)
            push!(test, :in)

            # Intersection
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (intersect!(scratch1, full); intersect!(scratch2, set)))
            push!(param, density)
            push!(test, :intersect)

            # Union
            copy!(scratch1, set)
            copy!(scratch2, empty)
            push!(times, @elapsed (union!(scratch1, empty); union!(scratch2, set)))
            push!(param, density)
            push!(test, :union)

            # Setdiff
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (setdiff!(scratch1, full); setdiff!(scratch2, set)))
            push!(param, density)
            push!(test, :setdiff)

            # Symdiff
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (symdiff!(scratch1, full); symdiff!(scratch2, set)))
            push!(param, density)
            push!(test, :symdiff)

            # Comparison
            set2 = copy(set)
            push!(times, @elapsed set == set2)
            push!(param, density)
            push!(test, :eq)

            push!(times, @elapsed set < set2)
            push!(param, density)
            push!(test, :lt)

            push!(times, @elapsed set <= set2)
            push!(param, density)
            push!(test, :lte)

            # Popping
            push!(times, @elapsed for e in elts pop!(set, e); end)
            push!(param, density)
            push!(test, :pop!)
        end
    end
    (times,param,test)
end

function test_size{T}(::Type{T}, N)
    times = Float64[]
    param = Float64[]
    test = Symbol[]
    sizehint!(times, N*10*11)
    sizehint!(param, N*10*11)
    sizehint!(test, N*10*11)
    for sz=round.(Int,logspace(1,5,9))
        println("Size: $sz")
        elts = 1:sz
        set = sizehint!(T(),sz)
        scratch1 = sizehint!(T(), sz)
        scratch2 = sizehint!(T(), sz)
        empty = sizehint!(T(), sz)
        full = T(elts)
        for _=1:N
            # Isempty
            push!(times, @elapsed isempty(set))
            push!(param, sz)
            push!(test, Symbol("isempty (empty)"))

            # push!
            push!(times, @elapsed for e in elts push!(set, e); end)
            push!(param, sz)
            push!(test, :push!)

            # Isempty
            push!(times, @elapsed isempty(set))
            push!(param, sz)
            push!(test, Symbol("isempty (full)"))

            # Iteration (sum)
            push!(times, @elapsed sum(set))
            push!(param, sz)
            push!(test, :sum)

            # Containment
            push!(times, @elapsed for e=1:sz e in set end)
            push!(param, sz)
            push!(test, :in)

            # Intersection
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (intersect!(scratch1, full); intersect!(scratch2, set)))
            push!(param, sz)
            push!(test, :intersect)

            # Union
            copy!(scratch1, set)
            copy!(scratch2, empty)
            push!(times, @elapsed (union!(scratch1, empty); union!(scratch2, set)))
            push!(param, sz)
            push!(test, :union)

            # Setdiff
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (setdiff!(scratch1, full); setdiff!(scratch2, set)))
            push!(param, sz)
            push!(test, :setdiff)

            # Symdiff
            copy!(scratch1, set)
            copy!(scratch2, full)
            push!(times, @elapsed (symdiff!(scratch1, full); symdiff!(scratch2, set)))
            push!(param, sz)
            push!(test, :symdiff)

            # Comparison
            set2 = copy!(scratch1, set)
            push!(times, @elapsed set == set2)
            push!(param, sz)
            push!(test, :eq)

            push!(times, @elapsed set < set2)
            push!(param, sz)
            push!(test, :lt)

            push!(times, @elapsed set <= set2)
            push!(param, sz)
            push!(test, :lte)

            # Popping
            push!(times, @elapsed for e in elts pop!(set, e); end)
            push!(param, sz)
            push!(test, :pop!)
        end
    end
    (times,param,test)
end

##
test_density(Base.IntSet, 1)
test_density(IntSets.IntSet, 1)
test_size(Base.IntSet, 1)
test_size(IntSets.IntSet, 1)

##
(time, param, fn) = test_density(Base.IntSet, REPS)
base = Any[time param fn fill!(Array{Symbol}(size(time)), :Base)]
(time, param, fn) = test_density(IntSets.IntSet, REPS)
test = Any[time param fn fill!(Array{Symbol}(size(time)), :IntSets)]
df = [base; test]
writecsv("perf_density.csv", df)

(time, param, fn) = test_size(Base.IntSet, REPS)
base = Any[time param fn fill!(Array{Symbol}(size(time)), :Base)]
(time, param, fn) = test_size(IntSets.IntSet, REPS)
test = Any[time param fn fill!(Array{Symbol}(size(time)), :IntSets)]
df = [base; test]
writecsv("perf_size.csv", df)


text:
module LMDB_Common
    using LMDB
    using Test

    import LMDB.MDBValue

    @test LMDB.version()[1] >= v"0.9.15"

    # LMDBError
    ex = LMDBError(Cint(0))
    @test_throws LMDBError throw(ex)
    @test ex.code == 0

    # MDBValue
    val = "abcd" # string
    mdb_val_ref = Ref(MDBValue(val))
    mdb_val = mdb_val_ref[]
    # @test val == unsafe_string(convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)
    @test val == convert(String, mdb_val_ref)

    val = [1233] # dense array
    T = eltype(val)
    val_size = sizeof(val)
    mdb_val_ref = Ref(MDBValue(val))
    mdb_val = mdb_val_ref[]
    @test val_size == mdb_val.size
    nvals = floor(Int, mdb_val.size/sizeof(T))
    value = unsafe_wrap(Array, convert(Ptr{T}, mdb_val.data), nvals)
    @test val == value
    @test val == convert(Vector{Int}, mdb_val_ref)

    val = [0x0003, 0xff45]
    val_size = sizeof(val)
    T = eltype(val)
    mdb_val_ref = Ref(MDBValue(val))
    mdb_val = mdb_val_ref[]
    @test val_size == mdb_val.size
    nvals = floor(Int, mdb_val.size/sizeof(T))
    value = unsafe_wrap(Array, convert(Ptr{T}, mdb_val.data), nvals)
    @test val == value
    @test val == convert(Vector{UInt16}, mdb_val_ref)

    struct TestType
        i::Int
        j::Char
    end
    val = TestType(1,'a') # struct
    val_size = sizeof(val)
    T = typeof(val)
    @test_throws ErrorException MDBValue(val)
    val = [val]
    mdb_val_ref = Ref(MDBValue(val))
    mdb_val = mdb_val_ref[]
    @test val_size == mdb_val.size
    nvals = floor(Int, mdb_val.size/sizeof(T))
    value = unsafe_wrap(Array, convert(Ptr{T}, mdb_val.data), nvals)
    @test val == value
    @test val == convert(Vector{T}, mdb_val_ref)
    @test val[1] == convert(T, mdb_val_ref)
end


text:
using GraphMatFun, LsqFit, LinearAlgebra, Random
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=7
target_n=200 # Measure error wrt this discr
n=50;        # Optimize wrt this discr
rho=6.4;  # Increased SID rho m=7

its=5;
droptol0=1e-10;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;
base.params[:kickit_mode]=1;

oldrho="3_59";
oldrho2="2_7";
# First the optimization free
datadir=joinpath("..","..","data","exp");
filename=joinpath(datadir,"exp_sid_m$(m).cgr");
sid_org=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)

#filename=joinpath(datadir,"exp_bbc_m$(m)_rho$oldrho.cgr");
#bbc_org=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)

filename=joinpath(datadir,"exp_sastre_m6.cgr");
sastre_org=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true,
                            scale_and_square=true)

#sastre_org.cref=get_degopt_crefs(sastre_org.graph);

filename=joinpath(datadir,"exp_ps_m$(m)_rho$oldrho.cgr");
ps_org=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)

filename=joinpath(datadir,"exp_mono_m$(m)_rho$oldrho.cgr");
mono_org=init_state_file!(deepcopy(base),:mono,filename,showmeta=true)

filename=joinpath(datadir,"exp_mono_m$(m-1)_opt_rho$(oldrho2).cgr");
mono0=init_state_file!(deepcopy(base),:mono,filename,showmeta=true,scale_and_square=true)
(mono,simlist,commandlist)=
     run_sequence(mono0,"ddddddddddddddddddddddddddssSSGGGSSddSDDDSddSSddSddgggggggSSSGGGGGGGSSggggggSSSSS2SGGGGGG22S2ddddggggggggSSGGGGGGGGS2zq");

filename=joinpath(datadir,"exp_ps_m$(m)_opt_rho$oldrho.cgr");
ps0=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)
(ps,simlist,commandlist)=
     run_sequence(ps0,
     "sdddddddddddddddddSddSdddddSddSddSSGGGSdddddSSSddddSSSddSdSdgggggggggsSSssssGGGGGGGGGSSSkSSdggggggggSSddSgggggSSSSSSSGGGGSSSSSSSGGGGGSGGGGSSSssszsddSSSsSdgggggggggggSDSSGGGGGSGGGGGGSSddggggggggggSDsssGGGGSGdgggggggSSS2GGGSGGGGGGSGGGS22Szq");

filename=joinpath(datadir,"exp_sastre_m$(m)_opt_rho$oldrho.cgr");
sastre0=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true)
(sastre,simlist,commandlist)=
     run_sequence(sastre0,
                  "ddddddddddsdddddddddddssGGGssddddddSSSSddSS2ddSdggggggggdSSSSGGGGGGGGS22zSsSSddSSSSSz2ddSddSdSgggggggggdSSSGGGGGGGGGS222ddgggggggggSgggggSSDSSGGGGGGGSGGGSdGGGGDS22zq");

filename=joinpath(datadir,"exp_sid_m$(m)_opt_rho$oldrho.cgr");
#filename=joinpath(datadir,"exp_sid_m$(m).cgr");
sid0=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)
(sid,simlist,commandlist)=
     run_sequence(sid_org,
     "q")
#ddddsddddskssssGsssdddddddddssdddsssddsddssssddsddsddsssssGGssdddsdssssdgggggggggsdsssGGGsGGGGGGGgssdddggggggggssssssGGGGGGGGsssdddsddsddgggggggggsssssssGGGGGGGGGsssdddsddggggggssssGGGGGGssssssssztzt2tSddgggggggggsssssGGGGGGGSGGS222tztSStddggggggggggsssS2sssGGGGGSGGGGGS22StztSStzSSddSSSSSSStSddsdgggggggggSSSSSGGGGS2SS22GGGGG22SStT
#filename=joinpath(datadir,"exp_bbc_m$(m)_opt_rho$oldrho.cgr");
#bbc0=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)
#(bbc,simlist,commandlist)=
#     run_sequence(bbc0,
#     "sdddddsdddddssssssddddsssssssddddsssssssssssssdddssssssssssddddssssssssdsggggggsdssssssssssssssssGGGGGsssGGGsssq")
#

include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# A simple compiler implementation
using ClangCompiler

const CC = ClangCompiler

struct SimpleCompiler
    ctx::CC.Context
    instance::CC.CompilerInstance
end

function create_simple_compiler(src::String, args::Vector{String}; diag_show_colors=true)
    ctx = CC.Context()
    instance = CC.CompilerInstance()
    # diagnostics
    CC.setShowPresumedLoc(instance, true)
    CC.setShowColors(instance, diag_show_colors)
    CC.createDiagnostics(instance)
    diag = CC.getDiagnostics(instance)
    # invocation
    # do not emit `__dso_handle` etc.
    insert!(args, length(args), "-fno-use-cxa-atexit")
    invok = CC.create_compiler_invocation_from_cmd(src, args, diag)
    CC.setInvocation(instance, invok)
    CC.setTargetAndLangOpts(instance)
    # source
    CC.createFileManager(instance)
    CC.createSourceManager(instance)
    CC.setMainFileID(instance, src)
    # preprocessor & AST & sema
    CC.createPreprocessor(instance)
    CC.createASTContext(instance)
    CC.setASTConsumer(instance, CC.CreateLLVMCodeGen(instance, ctx))
    CC.createSema(instance)
    return SimpleCompiler(ctx, instance)
end

function dispose(x::SimpleCompiler)
    CC.dispose(x.instance)
    CC.dispose(x.ctx)
end

function compile(x::SimpleCompiler)
    CC.parse(x.instance) || error("failed to parse the source code.")
    m = CC.get_llvm_module(x.instance)
    m.ref == C_NULL && error("failed to generate IR.")
    return m
end

# run
src = joinpath(@__DIR__, "..", "sample.cpp")
args = get_compiler_args()
cc = create_simple_compiler(src, args)
m = compile(cc)
dispose(cc)


text:
function gradcheckGPU(netGPU;showgrad=false)
# It's a bit painful to write a gradchecker for the GPU.
# We'll therefore check whether the GPU gradient matches the CPU gradient (which we can independently check is OK).

println("GRADCHECK GPU")

#netCPU=deepcopy(netGPU)
#for i in netCPU.node
#    netCPU.value[i]=to_host(netGPU.value[i])
#    netCPU.gradient[i]=to_host(netGPU.gradient[i])
#end
netCPU=convert(netGPU,"CPU",Float64)

netCPU=compile(netCPU;gpu=false,eltype=Float64)

ADforward!(netCPU); ADbackward!(netCPU); #println("evaluated CPU")
ADforward!(netGPU); ADbackward!(netGPU); #println("evaluated GPU")

    tol=1e-5
    for n in netGPU.validnodes # loop over the nodes in the network
        println("------------------------------------------------------------------------------")
        println("node $n:")

        D=length(netCPU.value[n])
        diff=sum(abs(netCPU.value[n]-to_host(netGPU.value[n])))/D
        reldiff=sum(abs((netCPU.value[n]-to_host(netGPU.value[n]))./(realmin()+netCPU.value[n])))/D
        println("absolute difference between CPU and GPU value=$diff")
        println("relative difference between CUP and GPU value=$reldiff")

        if diff>tol
            if showgrad
                println("CPU value:"); println(netCPU.value[n])
                println("GPU value:"); println(netGPU.value[n])
            end
            print_with_color(:red,"failed: mismatch more than $tol\n")
        else
            print_with_color(:green,"passed\n")
        end

        if any(node[n].takederivative)
            diff=sum(abs(netCPU.gradient[n]-to_host(netGPU.gradient[n])))/D
            reldiff=sum(abs((netCPU.gradient[n]-to_host(netGPU.gradient[n]))./(realmin()+netCPU.gradient[n])))/D
            println("absolute difference between CPU and GPU gradient=$diff")
            println("relative difference between CUP and GPU gradient=$reldiff")

            if diff>tol
                if showgrad
                    println("CPU gradient:"); println(netCPU.gradient[n])
                    println("GPU gradient:"); println(netGPU.gradient[n])
                end
                print_with_color(:red,"failed: mismatch more than $tol\n")
            else
                print_with_color(:green,"passed\n")
            end

        end
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function create_graph()
    v = Video(10, 10)
    Background(1:10, (args...)->O)
    return JGraph(true, 200, 200)
end

function create_weighted_graph()
    g = SimpleGraph()
    wg = WeightedGraph(g)
    return nv(wg) == 0
end

function create_reference_graph()
    g = SimpleGraph()
    rg = ReferenceGraph(g)
    return rg.graph !== rg.adjacency_graph
end

function check_graph_build()
    g = create_graph()
    return g.width == 200
end

function check__global_property_limits()
end

function check__global_layout()
    g = create_graph()
    v1 = GraphVertex(1, 1:2)
    v2 = GraphVertex(2, 1:2)
    e1 = GraphEdge(1, 2, 1:2)
    g.object.frames = 1:70

    layout_x, layout_y = spring_layout(g.graph.adjacency_graph)

    JavisGraphs._global_layout(:video, g.object, 1)
    # for (idx, v) in enumerate(collect(zip(layout_x, layout_y)))
    for v in vertices(g)
        # cannot check exact position due to random seed
        # @test abs(((v.object.opts[:position][1].-(0.5, 0.5)).*(200, 200))[1]) == 1
        @test -1 <= ((v.object.opts[:position][1]/200)+0.5)[1] <= 1
    end
    JavisGraphs._global_layout(:video, g.object, 12)
    for v in vertices(g)
        @test -1 <= ((v.object.opts[:position][1]/200)+0.5)[1] <= 1
    end
end


text:
module HMC
using Flux

struct Linear
  W
  b
end

function Linear(; dims, bias::Bool=true)
  W = randn(dims)
  if bias
    b = randn(1, dims[2])
  else
    b = 0
  end
  Linear(param(W), param(b))
end

(m::Linear)(x) = x * m.W .+ m.b

function runtest(niters; num_obs=10000, lr=.01)
  # Simulate random predictors
  X = randn(num_obs, 1)

  # Set parameter truth
  W_true = 2
  b_true = 1

  # function to simulate response
  f(X, w, b) = w * X .+ b + randn(size(X, 1), 1) * .1

  # simulate response
  y = f(X, W_true, b_true)

  # numer of features
  num_features = size(X, 2)

  # build model
  model = HMC.Linear(dims=(num_features, 1))

  # loss function
  loss(x, y) = Flux.mean((model(x) - y) .^ 2)

  # parameters to track
  ps = Flux.Params([getfield(model, fname)
                    for fname in fieldnames(typeof(model))])

  # gradient descent optimizer
  opt = Descent(lr)

  # Training loop
  for i in 1:niters
    # NOTE: These versions are equivalent

    # Version I
    # Flux.train!(loss, ps, [(X, y)], opt)

    # Version II
    # gs = Tracker.gradient(() -> loss(X, y), ps)
    # Flux.Tracker.update!(opt, ps, gs)

    # Version III. Could be used in an HMC step.
    # See p.14 of https://arxiv.org/pdf/1206.1901.pdf
    gs = Tracker.gradient(() -> loss(X, y), ps)
    for p in ps
      # Update parameter `p` to have value `p += -gs[p]*lr`
      Flux.Tracker.update!(p, -gs[p] * lr)
    end
  end

  # Print result
  println("model: $(model) | truth: $((W_true, b_true)) | loss: $(loss(X, y))")
end

end  # end of module

#= Time the thing
@time HMC.runtest(100, num_obs=100, lr=.1)
=#