text:
"""Gradient noise generation (e.g. Perlin noises)
"""

function adjacent(p)
    fp1, fp2 = floor.(p)
    cp1, cp2 = ceil.(p)
    @SArray Int[fp1 cp1 fp1 cp1; fp2 fp2 cp2 cp2]
end

function lerp_weights(w)
    w1, w2 = w
    @SVector [1 - w1 - w2 + w1 * w2, w1 - w1 * w2, w2 - w1 * w2, w1 * w2]
end

function perlin!(A::AbstractArray, grid::AbstractArray)
    resolution, scale = size(A), size(grid)
    ranges = map.(Ref(x -> (1, x)), [resolution, scale])
    n = ndims(grid)
    neighbors = zeros(Int, n, 2^n)
    weights = zeros(2^n)
    # dot_prod_tmp = @SArray zeros(n)
    dot_prod_tmp = zeros(n)
    scaled_indices = zeros(n)
    # w = @SArray zeros(n)
    w = zeros(n)
    it_indices = collect(Base.Iterators.product(range.(Ref(1), resolution, step=1)...))

    # for indices in Base.Iterators.product(range.(Ref(1), resolution, step=1)...)
    #     scaled_indices .= remap.(indices, ranges...)
    #     @. w = fade(scaled_indices - floor(scaled_indices))
    #     lerp_weights!(weights, w)
    #     adjacent!(neighbors, scaled_indices)
    #     value = 0
    #     for i in 1:2^n
    #         neigh = neighbors[:, i]
    #         grad = getindex(grid, neigh...)
    #         @. dot_prod_tmp = (scaled_indices - neigh) * grad
    #         value += sum(dot_prod_tmp) * weights[i]
    #     end
    #     # distances = map(x -> scaled_indices .- x, neighbors)
    #     # dot_products = dot.(map(x -> getindex(grid, x...), neighbors), distances)
    #     # print(dot_products)
    #     # value = dot(dot_products, weights)
    #     setindex!(A, value, indices...)
    # end

    perlin_!(A, grid, weights, neighbors, scaled_indices, dot_prod_tmp, w, n, ranges, it_indices)

    A
end

function perlin_!(A, grid, weights, neighbors, scaled_indices, dot_prod_tmp, w, n, ranges, it_indices)
    for indices in it_indices
        scaled_indices .= remap.(indices, ranges[1], ranges[2])
        @. w = fade(scaled_indices - floor(scaled_indices))
        weights = lerp_weights(w)
        neighbors .= adjacent(scaled_indices)
        value = 0
        for i in 1:2^n
            @views neigh = neighbors[:, i]
            grad = getindex(grid, neigh[1], neigh[2])
            @. dot_prod_tmp = (scaled_indices - neigh) * grad
            value += sum(dot_prod_tmp) * weights[i]
        end
        # distances = map(x -> scaled_indices .- x, neighbors)
        # dot_products = dot.(map(x -> getindex(grid, x...), neighbors), distances)
        # print(dot_products)
        # value = dot(dot_products, weights)
        setindex!(A, value, indices...)
    end
end

# function perlin!(A::AbstractArray, grid::AbstractArray)
#     resolution, scale = size(A), size(grid)
#     ranges = map.(Ref(x -> (1, x)), [resolution, scale])
#     for indices in collect(Base.Iterators.product(range.(Ref(1), resolution, step=1)...))
#         scaled_indices = remap.(indices, ranges...)
#         subgrid_position = (scaled_indices .- floor.(scaled_indices))
#         w1, w2 = fade.(subgrid_position)
#         weights = (1 - w1 - w2 + w1 * w2, w1 - w1 * w2, w2 - w1 * w2, w1 * w2)
#         neighbors = adjacent(scaled_indices)
#         value = 0
#         for (neigh, weight) in zip(neighbors, weights)
#             grad = getindex(grid, neigh...)
#             value += dot(scaled_indices .- neigh, grad) * weight
#         end
#         # distances = map(x -> scaled_indices .- x, neighbors)
#         # dot_products = dot.(map(x -> getindex(grid, x...), neighbors), distances)
#         # print(dot_products)
#         # value = dot(dot_products, weights)
#         setindex!(A, value, indices...)
#     end

#     A
# end

function perlin(resolution, grid::AbstractArray)
    noise = zeros(resolution...)
    perlin!(noise, grid)
end

function create_gradients(scale)
    grid = zeros(scale...)
    grid = map(x -> normalize(rand(length(scale))), grid)
end

function perlin(resolution, scale)
    grid = create_gradients(scale)
    perlin(resolution, grid)
end

function perlin2d!(A::AbstractArray, grid::AbstractArray)
    resolution, scale = size(A), size(grid)
    for i in 1:resolution[1]
        for j in 1:resolution[2]
            x, y = remap.((i, j), map.(Ref(x -> (1, x)), [resolution, scale])...)
            x0 = Int(floor(x))
            x1 = Int(ceil(x))
            y0 = Int(floor(y))
            y1 = Int(ceil(y))
            sx = fade(x - x0)
            sy = fade(y - y0)

            n0 = dot_grid_gradient(grid, x0, y0, x, y)
            n1 = dot_grid_gradient(grid, x1, y0, x, y)
            ix0 = lerp([n0, n1], sx)

            n2 = dot_grid_gradient(grid, x0, y1, x, y)
            n3 = dot_grid_gradient(grid, x1, y1, x, y)
            ix1 = lerp([n2, n3], sx)
            # println((n0, n1, n2, n3))

            A[i, j] = lerp([ix0, ix1], sy)
        end
    end
    A
end

function dot_grid_gradient(grid, i, j, x, y)
    dx = x - i
    dy = y - j
    dot(grid[i, j], [dx, dy])
end
lerp(vals, w) = vals[1] * (1 - w) + vals[2] * w
fade(x) = 6x^5 - 15x^4 + 10x^3


text:
# ---
# title: 1020. Number of Enclaves
# id: problem1020
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/number-of-enclaves/description/>
# hidden: true
# ---
#
# Given a 2D array `A`, each cell is 0 (representing sea) or 1 (representing
# land)
#
# A move consists of walking from one land square 4-directionally to another
# land square, or off the boundary of the grid.
#
# Return the number of land squares in the grid for which we **cannot** walk off
# the boundary of the grid in any number of moves.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[0,0,0,0],[1,0,1,0],[0,1,1,0],[0,0,0,0]]
#     Output: 3
#     Explanation:
#     There are three 1s that are enclosed by 0s, and one 1 that isn't enclosed because its on the boundary.
#
# **Example 2:**
#
#
#
#     Input: [[0,1,1,0],[0,0,1,0],[0,0,1,0],[0,0,0,0]]
#     Output: 0
#     Explanation:
#     All 1s are either on the boundary or can reach the boundary.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 500`
#   2. `1 <= A[i].length <= 500`
#   3. `0 <= A[i][j] <= 1`
#   4. All rows have the same size.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
s=time()
spl=split(ARGS[1])
expr=spl[1]
a=float(spl[2])
b=float(spl[3])
n=Integer(float(spl[4]))
include_string(string("f(x)=",expr))
l=abs(b-a)
function serie(i)
	return f(a+(i*l)/(n))*(l/n)
end
r=0
for i in 1:n
	r+=serie(i-1)
end
println(r)
println(time()-s)

max_stars_repo_name: santiagohenao/misc

max_issues_repo_name: santiagohenao/misc

max_forks_repo_name: santiagohenao/misc


text:
# ---
# title: 1622. Fancy Sequence
# id: problem1622
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Design
# link: <https://leetcode.com/problems/fancy-sequence/description/>
# hidden: true
# ---
#
# Write an API that generates fancy sequences using the `append`, `addAll`, and
# `multAll` operations.
#
# Implement the `Fancy` class:
#
#   * `Fancy()` Initializes the object with an empty sequence.
#   * `void append(val)` Appends an integer `val` to the end of the sequence.
#   * `void addAll(inc)` Increments all existing values in the sequence by an integer `inc`.
#   * `void multAll(m)` Multiplies all existing values in the sequence by an integer `m`.
#   * `int getIndex(idx)` Gets the current value at index `idx` (0-indexed) of the sequence **modulo** `109 + 7`. If the index is greater or equal than the length of the sequence, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["Fancy", "append", "addAll", "append", "multAll", "getIndex", "addAll", "append", "multAll", "getIndex", "getIndex", "getIndex"]
#     [[], [2], [3], [7], [2], [0], [3], [10], [2], [0], [1], [2]]
#     **Output**
#     [null, null, null, null, null, 10, null, null, null, 26, 34, 20]
#
#     **Explanation**
#     Fancy fancy = new Fancy();
#     fancy.append(2);   // fancy sequence: [2]
#     fancy.addAll(3);   // fancy sequence: [2+3] -> [5]
#     fancy.append(7);   // fancy sequence: [5, 7]
#     fancy.multAll(2);  // fancy sequence: [5*2, 7*2] -> [10, 14]
#     fancy.getIndex(0); // return 10
#     fancy.addAll(3);   // fancy sequence: [10+3, 14+3] -> [13, 17]
#     fancy.append(10);  // fancy sequence: [13, 17, 10]
#     fancy.multAll(2);  // fancy sequence: [13*2, 17*2, 10*2] -> [26, 34, 20]
#     fancy.getIndex(0); // return 26
#     fancy.getIndex(1); // return 34
#     fancy.getIndex(2); // return 20
#
#
#
#
# **Constraints:**
#
#   * `1 <= val, inc, m <= 100`
#   * `0 <= idx <= 105`
#   * At most `105` calls total will be made to `append`, `addAll`, `multAll`, and `getIndex`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
##############################################################################
##
## Utils
##
##############################################################################

macro isok(A)
    :($A ==  Int32(1) || throw(CHOLMODException("")))
end

# Update B as A'
function transpose_unsym_!{Tv<: VTypes}(A::Sparse{Tv}, values::Integer, B::Sparse{Tv})
    @isok ccall((@cholmod_name("transpose_unsym", SuiteSparse_long),:libcholmod),
        Cint,
            (Ptr{C_Sparse{Tv}}, Cint, Ptr{SuiteSparse_long}, Ptr{SuiteSparse_long}, Csize_t, Ptr{C_Sparse{Tv}}, Ptr{UInt8}),
            A.p, values, C_NULL, C_NULL, 0, B.p, common())
    return B
end

# Update B as Sparse(A)
function Sparse!{Tv<:VTypes,Ti<:ITypes}(A::SparseMatrixCSC{Tv,Ti}, B::Sparse{Tv})
    s = unsafe_load(B.p)
    unsafe_copy!(s.x, pointer(A.nzval), length(A.nzval))
    @isok check_sparse(B)
    return B
end

##############################################################################
##
## Constructor
##
##############################################################################

type SparseCholeskyAllocatedSolver{Tv, Ti <: Integer, Tx <: AbstractVector} <: AbstractAllocatedSolver
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    v::Tx
    sparseJ::Sparse{Tv}
    sparseJt::Sparse{Tv}
    F::Factor{Tv}
    cm::Array{UInt8, 1}
    function SparseCholeskyAllocatedSolver(colptr, rowval, v, sparseJ, sparseJt, F, cm)
        size(sparseJ) == (size(sparseJt, 2),  size(sparseJt, 1)) || throw(DimensionMismatch("J and Jt should have transposed dimension"))
        new(colptr, rowval, v, sparseJ, sparseJt, F, cm)
    end
end

function SparseCholeskyAllocatedSolver{Tv, Ti <: Integer, Tx <: AbstractVector}(colptr::Vector{Ti}, rowval::Vector{Ti}, v::Tx, sparseJ::Sparse{Tv}, sparseJt::Sparse{Tv}, F::Factor{Tv}, cm::Array{UInt8, 1})
    SparseCholeskyAllocatedSolver{Tv, Ti, Tx}(colptr, rowval, v, sparseJ, sparseJt, F, cm)
end

function AbstractAllocatedSolver{Tx, Ty, Tf, TJ <: SparseMatrixCSC , Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg, solver::Cholesky)
    colptr = deepcopy(nls.J.colptr)
    rowval = deepcopy(nls.J.rowval)
    sparseJ = Sparse(nls.J)
    sparseJt = transpose_(sparseJ, 2)
    cm = defaults(common())
    set_print_level(cm, 0)
    unsafe_store!(common_final_ll, 1)
    F = analyze(sparseJt, cm)
    return SparseCholeskyAllocatedSolver(colptr, rowval, _zeros(nls.x), sparseJ, sparseJt, F, cm)
end

##############################################################################
##
## solve J'J \ J'y by Cholesky
##
##############################################################################

function A_ldiv_B!(x::AbstractVector, J::SparseMatrixCSC, y::AbstractVector, A::SparseCholeskyAllocatedSolver)
    colptr, rowval, v, sparseJ, sparseJt, F, cm =
    A.colptr, A.rowval, A.v, A.sparseJ, A.sparseJt, A.F, A.cm

    # check symbolic structure is the same
    if colptr != J.colptr || rowval != J.rowval
        error("The symbolic structure of the Jacobian has been changed. Either (i) rewrite g! so that it does not modify the structure of J (see Julia issue #9906) (ii) use solver = :iterative rather than solver = :cholesky")
    end

    Sparse!(J, sparseJ)
    transpose_unsym_!(sparseJ, 2, sparseJt)
    factorize_p!(sparseJt, 0, F, cm)
    Ac_mul_B!(v, J, y)
    # !! there is a memory allocation here
    copy!(x, F \ v)
    return x, 1
end


text:
#!/usr/bin/env -S julia -t 10
using FASTX,BioAlignments,BioSequences,ArgParse,CodecZlib

"""
    getSNPListFromAlignment(aln)

Takes an alignment and produces a list of the SNPs and indels in the alignment.
"""
function getSNPAndListFromAlignment(aln,chrom,chrStart)
	index = 1

	# This is a bit much, but everything else in the alignment
	# interface is super clean...
	#refstart=aln.a.aln.anchors[1].refpos+1
	refstart=aln.a.aln.anchors[1].refpos+chrStart

	snplist=Any[]
	for (q,r) in collect(aln)
		# matches and insertions in query advance reference index
		if r == q || q == DNA_Gap
			index+=1
		# Mismatches are saved unless reference is a gap.
		elseif r != q && r!= DNA_Gap
			push!(snplist,(chrom,index+refstart-1,r,q))
			index+=1
		else
			println("DEBUG:",r,q)
		end
	end
	return(snplist)
end

function parseargs()
	s = ArgParseSettings()
	s.description=
	"""
	Align two sequences from multiple-fasta file against each other.
	"""
	@add_arg_table s begin
		"--aseq","-a"
			help = "ID of first sequence"
			required=true
		"--bseq","-b"
			help = "ID of second sequence"
			required=true
		"--querydb","-q"
			help = "Query database of sequences"
			nargs = '+'
			required=true
		"--chromstart","-c"
			help = "Chromosome start"
			default = "0"
			required = false
	end
	return parse_args(s)
end

function readRecordDict(files,name2seq)
	SEQ_TYPE = FASTA
	for fileName in files
		if occursin("fq",fileName) || occursin("fastq",fileName)
			SEQ_TYPE = FASTQ
		end

		reader = SEQ_TYPE.Reader(GzipDecompressorStream(open(fileName)))
		print(stderr,"Reading ",fileName,"...")
		for record in reader
			name2seq[SEQ_TYPE.identifier(record)]=record
		end
		println(stderr,length(name2seq)," sequences read.")
	end

	return(SEQ_TYPE)
end

function main()
	args = parseargs()

	allelename2seq = Dict()
	SEQ_TYPE= readRecordDict(args["querydb"],allelename2seq)

	scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=-1);
	# Align hla type sequence to reference region
	arecord=allelename2seq[args["aseq"]]
	aseq = SEQ_TYPE.sequence(arecord)
	brecord=allelename2seq[args["bseq"]]
	bseq = SEQ_TYPE.sequence(brecord)

	# Don't know how they are oriented, so align both ways and take better score...
	respos = pairalign(LocalAlignment(),aseq,bseq,scoremodel)
	posscore = score(respos)

	revaseq = reverse_complement(aseq)
	resneg = pairalign(LocalAlignment(),revaseq,bseq,scoremodel)
	negscore = score(resneg)

	println(stderr,"Alignment score:        \t",posscore)
	println(stderr,"Reverse alignment score:\t",negscore)

	println(args["aseq"])
	println(args["bseq"])
	if negscore > posscore
		println(resneg)
	else
		println(respos)
#		chrom="chr6"
#		chrStart = parse(Int,args["chromstart"])
#		if chrStart != 0
#			println("SNPS:")
#			snps = getSNPListFromAlignment(alignment(respos),chrom,chrStart)
#			for (chrom,loc,r,q) in snps
#				println(chrom,":",string(loc),":",string(r),":",string(q))
#			end
#		end
	end
end

main()

# bedfile B Start: 31351875
# bedfile B End: 31359179
# Size:
# minimap B start: 31353599
# minimap B end:  31357442

# 50 sequences:
# real	0m4.318s
# 100 alignments:
# real	0m6.165s
# 500 alignments:
# real	0m10.336s.  Not bad

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
import Vulkan

import Vulkan: handle

const Vk = Vulkan
using .Vk.ResultTypes: unwrap

instance = Instance(Vk.InstanceCreateInfo([], []))
physical_device = first(unwrap(Vk.enumerate_physical_devices(instance)))
device = Device(physical_device, Vk.DeviceCreateInfo([Vk.DeviceQueueCreateInfo(0, [1.0])], [], []))

sem = Semaphore(device, Vk.SemaphoreCreateInfo())
event = Event(device)


text:
using NeuralVerification, LazySets, Test, LinearAlgebra
import NeuralVerification: ReLU, Id
include("runtime_aux.jl")

# Group 1

# Small Problem: problem_small_HH
# Solvers: ExactReach(), Ai2(), MaxSens(resolution = 0.6)

# ACAS Problem: problem_acas_HH
# Solvers: ExactReac() - timesout, Ai2(), MaxSens(1.0, false)

# Group 2

# Small Problem: problem_small_HrH
# Solvers: NSVerify(GLPKSolverMIP(), 1000.0), MIPVerify(GLPKSolverMIP()),
# 	ILP(GLPKSolverMIP(), 1),

# ACAS Problem: problem_acas_HrH
# Solvers: NSVerify(), MIPVerify(), ILP()

# Group 3

# Small Problem: problem_small_HH
# Solvers: Duality(), ConvDual(),

# ACAS Problem: problem_acas_HrH
# Solvers: Duality(), ConvDual()

# Group 4

# Small Problem: problem_small_HrH
# Solvers: FastLin(), FastLip(), MIPVerify(), ILP()

# ACAS Problem: problem_acas_HrH
# Solvers: FastLin(), FastLip(), MIPVerify(), ILP()

# Group 5

# Small Problem: problem_small_RR
# Solvers: ReluVal(), DLV(), Sherlock(), BaB()

# ACAS Problem: problem_acas1_RR
# Solvers: ReluVal(), DLV(), Sherlock(), BaB()

# Group 6

# Small Problem: problem_small_HH
# Solvers: Planet(), Reluplex(), ReluVal()

# ACAS Problem: problem_acas_HrH
# Solvers: Planet(), Reluplex(), ReluVal()

test_runtime(problem_small_HH, ExactReach())


text:
@json struct Plan
    name::String
    space::Int
    private_repos::Int
    collaborators::Int
end

@json struct User
    login::String
    id::Int
    node_id::String
    avatar_url::String
    gravatar_id::String
    url::String
    html_url::String
    followers_url::String
    following_url::String
    gists_url::String
    starred_url::String
    subscriptions_url::String
    organizations_url::String
    repos_url::String
    events_url::String
    received_events_url::String
    type::String
    site_admin::Bool
    name::String
    company::String
    blog::String
    location::String
    email::String
    hireable::Bool
    bio::String
    public_repos::Int
    public_gists::Int
    followers::Int
    following::Int
    created_at::DateTime
    updated_at::DateTime
    private_gists::Int
    total_private_repos::Int
    owned_private_repos::Int
    disk_usage::Int
    collaborators::Int
    two_factor_authentication::Bool
    plan::Plan
end

endpoint(::GitHubAPI, ::typeof(get_user)) = Endpoint(:GET, "/user")
@not_implemented(::GitHubAPI, ::typeof(get_user), ::Int64)
endpoint(::GitHubAPI, ::typeof(get_user), name::AStr) = Endpoint(:GET, "/users/$name")
into(::GitHubAPI, ::typeof(get_user)) = User

endpoint(::GitHubAPI, ::typeof(get_users)) = Endpoint(:GET, "/users")
into(::GitHubAPI, ::typeof(get_users)) = Vector{User}

@not_implemented(api::GitHubAPI, ::typeof(update_user), id::Integer)

endpoint(::GitHubAPI, ::typeof(update_user)) = Endpoint(:PATCH, "/user")
into(::GitHubAPI, ::typeof(update_user)) = User

@not_implemented(::GitHubAPI, ::typeof(create_user))

@not_implemented(::GitHubAPI, ::typeof(delete_user), id::Integer)


text:
export JonesMorari

"""
    JonesMorari

A simple 4-state 2-input discrete time system.

This system has simple upper and lower bounds on the states and inputs.
It was originally proposed in Section 6 of [1].

[1] C. N. Jones and M. Morari, ‘The Double Description Method for the Approximation of Explicit MPC Control Laws’, in 47th IEEE Conference on Decision and Control (CDC), Cancun, Mexico, 2008, pp. 4724–4730, doi: 10.1109/CDC.2008.4738794.
"""
struct JonesMorari <: ControlBenchmarkProblem end

@addbenchmark JonesMorari "lin02" "Schur-stable system (Jones-Morari)" "Simple 4-state 2-input schur-stable system" Set([:linear, :stable, :discrete])

"""
    controlbenchmark( ::JonesMorari ) -> ( sys, stateBounds, inputBounds )

Create the components of the [`JonesMorari`](@ref) benchmark problem.

* `sys` is the StateSpace system.
* `stateBounds` are the limits on the state variables.
* `inputBounds` are the limits on the input variables.
"""
function controlbenchmark( ::JonesMorari )

    # System matrices
    A = [ 0.7 -0.1  0.0  0.0;
          0.2 -0.5  0.1  0.0;
          0.0  0.1  0.1  0.0;
          0.5  0.0  0.5  0.5 ]

    B = [ 0.0  0.1;
          0.1  1.0;
          0.1  0.0;
          0.0  0.0 ]

    C = Matrix{Float64}( I, 4, 4 )

    D = fill( 0.0, (4, 2) )

    # Bounds
    stateBounds = Bounds( [ Bound( -5.0, 5.0 );
                            Bound( -5.0, 5.0 );
                            Bound( -5.0, 5.0 );
                            Bound( -5.0, 5.0 ) ] )

    inputBounds = Bounds( [ Bound( -5.0, 5.0 );
                            Bound( -5.0, 5.0 ) ] )

    return (; :sys => StateSpace( A, B, C, D, 0.01 ),
              :stateBounds => stateBounds,
              :inputBounds => inputBounds )
end


text:
@testset "Utils" begin
    s = @inferred(LightGraphs.sample!([1:10;], 3))
    @test length(s) == 3
    for e in s
        @test 1 <= e <= 10
    end

    s = @inferred(LightGraphs.sample!([1:10;], 6, exclude=[1, 2]))
    @test length(s) == 6
    for e in s
        @test 3 <= e <= 10
    end

    s = @inferred(LightGraphs.sample(1:10, 6, exclude=[1, 2]))
    @test length(s) == 6
    for e in s
        @test 3 <= e <= 10
    end

    # tests if isbounded has the correct behaviour
    bounded_int_types = [Int8, Int16, Int32, Int64, Int128,
                         UInt8, UInt16, UInt32, UInt64, UInt128,
                         Int, Bool]
    unbounded_int_types = [BigInt, Signed, Unsigned, Integer, Union{Int8, UInt8}]
    for T in bounded_int_types
        @test LightGraphs.isbounded(T) == true
        @test LightGraphs.isbounded(T(0)) == true
    end
    for T in unbounded_int_types
        @test LightGraphs.isbounded(T) == false
        if isconcretetype(T)
            @test LightGraphs.isbounded(T(0)) == false
        end
    end
end

@testset "Generate Reduce" begin

    function make_vec(g::AbstractGraph{T}) where T<:Integer
        return Vector{T}(undef, nv(g))
    end

    g1 = StarGraph(5)

    for g in testgraphs(g1)
        s = @inferred(LightGraphs.generate_max_set(g, make_vec, 5))
        @test length(s) == 5
    end

    for g in testgraphs(g1)
        s = @inferred(LightGraphs.generate_max_set(g, make_vec, 5))
        @test length(s) == 5
    end
end

@testset "Greedy Contiguous Partition" begin

    p = @inferred(LightGraphs.greedy_contiguous_partition([1, 1, 1, 3], 2))
    @test p == [1:3, 4:4]

    p = @inferred(LightGraphs.greedy_contiguous_partition([1, 2, 3, 4, 5, 100, 1, 3, 1, 1], 3))
    @test p == [1:5, 6:6, 7:10]

    p = @inferred(LightGraphs.greedy_contiguous_partition([1, 1, 1, 1], 4))
    @test p == [1:1, 2:2, 3:3, 4:4]
end

@testset "Optimal Contiguous Partition" begin

    p = @inferred(LightGraphs.optimal_contiguous_partition([1, 1, 1, 3], 2))
    @test p == [1:3, 4:4]

    p = @inferred(LightGraphs.optimal_contiguous_partition([1, 2, 3, 4, 5, 100, 1, 3, 1, 1], 3))
    @test p == [1:5, 6:6, 7:10]

    p = @inferred(LightGraphs.optimal_contiguous_partition([1, 1, 1, 1], 4))
    @test p == [1:1, 2:2, 3:3, 4:4]
end


text:
# DTU02435 Decision Making under Uncertainty
# Assignment 3, Task 2
# Edward J. Xu
# May 5th, 2019
push!(LOAD_PATH, "$(homedir())/Desktop/StochasticAirplaneRent")
cd("$(homedir())/Desktop/StochasticAirplaneRent")
using JuMP
# using GLPKMathProgInterface
# using Gurobi
using CPLEX
using CSV

function processData()
    datf_numMax = CSV.read("Data/MaxPlanes.csv", header = ["p", "num"], delim = ',')
    datf_cost = CSV.read("Data/OperationCosts.csv", header = true, delim = ',')
    datf_cap = CSV.read("Data/PlaneCapacity.csv", header = ["p", "cap"], delim = ',')
    datf_route = CSV.read("Data/DemandsRoute.csv", header = true, delim = ',')
    datf_costAdmin = CSV.read("Data/CostAdquisition.csv", header = ["p", "cost"], delim = ',')
    datf_prob = CSV.read("Data/Probability.csv", header = ["prob"], delim = ',')
    datf_demand = CSV.read("Data/Scenario.csv", header = true, delim = ',')
    datf_expect = CSV.read("Data/ExpectedValue.csv", header = ["value"], delim = ',')
    # Calculate the demand matrix
    mat4_demand = zeros(4, 10, 12, 10)
    for i = 1: 4
        for j = 1: 10
            for w = 1: 10
                for t = 1: 12
                    mat4_demand[i, j, t, w] = 1000 * datf_demand[t, w] * datf_route[i, j+1]
                end
            end
        end
    end
    mat3_demandExpect = zeros(4, 10, 12)
    for i = 1: 4
        for j = 1: 10
            for t = 1: 12
                mat3_demandExpect[i, j, t] = 1000 * datf_expect.value[t] * datf_route[i, j+1]
            end
        end
    end
    return (datf_numMax.num, datf_cost, datf_cap.cap, datf_costAdmin.cost, mat4_demand,
            datf_prob.prob, mat3_demandExpect)
end

function solvePrintResult(model, variable_y, variable_m)
    println("---------------------------- 1/2. Begin Optimization ---------------------------")
    timeStart = time()
    solve(model)
    println("Elapsed time: $(time()-timeStart) seconds.")
    ## 2,  Print the Result
    objResult = getobjectivevalue(model)
    vec_yResult = getvalue(variable_y)
    vec_mResult = getvalue(variable_m)
    println("Result: obj = $(objResult)\n",
            "Result: vec_y = $(vec_yResult)\n",
            "Result: total unfulfilled passenger demand = $(sum(vec_mResult))\n",
            "---------------------------------- 2/2. Ending ---------------------------------")
    return vec_yResult
end

function optim_determin(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat3_demandExpect)
    ## 1,  Begin Optimization
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_y[1: 10] >= 0, Int)
    @variable(model, mat4_z[1: 4, 1: 10, 1: 10, 1: 12] >= 0, Int)
    @variable(model, mat3_m[1: 4, 1: 10, 1: 12] >= 0)
    @objective(model, Min, sum(vec_y[p] * vec_costAdmin[p] for p = 1: 10) +
        sum(datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j), (datf_cost[:p] .== p))][1] *
        mat4_z[i, j, p, t] for p = 1: 10, i = 1: 4, j = 1: 10, t = 1: 12) +
        sum(mat3_m[i, j, t] * 1000 for i = 1: 4, j = 1: 10, t = 1: 12))
    @constraint(model, [p = 1: 10], vec_y[p] <= vec_numMax[p])
    @constraint(model, [i = 1: 4, j = 1: 10, t = 1: 12], mat3_demandExpect[i, j, t] -
        sum([mat4_z[i, j, p, t] * vec_cap[p]] for p = 1: 10)[1] <= mat3_m[i, j, t])
    @constraint(model, [p = 1: 10, t = 1: 12], sum(mat4_z[i, j, p, t] for i = 1: 4, j = 1: 10) <= vec_y[p])
    vec_yResult = solvePrintResult(model, vec_y, mat3_m)
    return vec_yResult
end

function optim(vec_y, datf_cost, vec_cap, vec_costAdmin, mat_demand)
    ## Most basic function to optimize the arrangement of plane for any month in any scenario, with fixed fleet.
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, mat3_z[1:4, 1:10, 1:10] >= 0, Int)
    @variable(model, mat_m[1:4, 1:10] >= 0)
    @objective(model, Min, sum(datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j),
        (datf_cost[:p] .== p))][1] * mat3_z[i, j, p] for p = 1:10, i = 1:4, j = 1:10) + sum(mat_m) * 1000)
    @constraint(model, [i = 1:4, j = 1:10], mat_demand[i, j] - sum(mat3_z[i, j, p] * vec_cap[p] for p = 1:10)
        <= mat_m[i, j])
    @constraint(model, [p = 1:10], sum(mat3_z[i, j, p] for i = 1:4, j = 1:10) <= vec_y[p])
    solve(model)
    objResult = getobjectivevalue(model)
    vec_mResult = getvalue(mat_m)
    return (objResult, sum(vec_mResult))
end

function optim_determinStochas(vec_yResult, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    # Function to calculate the obj when first stage variable is fixed by solutions from deterministic program, using
    # stochastic model. The program is decomposed into programs for different months.
    mat_objResult = zeros(10, 12)
    vec_objResult = zeros(10)
    mat_mResult = zeros(10, 12)
    vec_mResult = zeros(10)
    ## 1,  Begin Optimization
    println("---------------------------- 1/2. Begin Calculation ----------------------------\n")
    timeStart = time()
    for w = 1:10
        for t = 1:12
            (mat_objResult[w, t], mat_mResult[w, t]) = optim(vec_yResult, datf_cost, vec_cap, vec_costAdmin,
                mat4_demand[:, :, t, w])
        end
        vec_objResult[w] = sum(mat_objResult[w, tt] for tt = 1:12) +
            sum(vec_yResult[p] * vec_costAdmin[p] for p = 1: 10)
        vec_mResult[w] = sum(mat_mResult[w, tt] for tt = 1:12)
    end
    objResult = sum(vec_objResult[w] * vec_prob[w] for w = 1:10)
    mResult = sum(vec_mResult)
    println("Elapsed time: $(time()-timeStart) seconds.\n",
            "Result: obj = $(objResult)\n",
            "Result: vec_y = $(vec_yResult)\n",
            "Result: total unfulfilled passenger demand = $(mResult)\n",
            "---------------------------------- 2/2. Ending ---------------------------------")
end

function optim_stochas(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    ## 1,  Begin Optimization
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_y[1: 10] >= 0, Int)
    @variable(model, mat5_z[1: 4, 1: 10, 1: 10, 1: 12, 1: 10] >= 0, Int)
    @variable(model, mat4_m[1: 4, 1: 10, 1: 12, 1: 10] >= 0)
    @objective(model, Min, sum(vec_y[p] * vec_costAdmin[p] for p = 1: 10) +
        sum((datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j), (datf_cost[:p] .== p))][1] *
        mat5_z[i, j, p, t, w] + mat4_m[i, j, t, w] * 1000) * vec_prob[w] for w = 1: 10, p = 1: 10, i = 1: 4,
        j = 1: 10, t = 1: 12))
    @constraint(model, [p = 1: 10], vec_y[p] <= vec_numMax[p])
    @constraint(model, [i = 1: 4, j = 1: 10, t = 1: 12, w = 1: 10], mat4_demand[i, j, t, w] -
        sum([mat5_z[i, j, p, t, w] * vec_cap[p]] for p = 1: 10)[1] <= mat4_m[i, j, t, w])
    @constraint(model, [p = 1: 10, t = 1: 12, w = 1: 10], sum(mat5_z[i, j, p, t, w] for i = 1: 4, j = 1: 10) <=
        vec_y[p])
    solvePrintResult(model, vec_y, mat4_m)
end

# function optim_stochas_2(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
#     ## 1,  Begin Optimization
#     model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
#     # model = Model(solver = GurobiSolver(Presolve=0))
#     # model = Model(solver = GLPKSolverMIP())
#     @variable(model, mat_y[1: 10, 1: 10] >= 0, Int)
#     @variable(model, mat5_z[1: 4, 1: 10, 1: 10, 1: 12, 1: 10] >= 0, Int)
#     @variable(model, mat4_m[1: 4, 1: 10, 1: 12, 1: 10] >= 0)
#     @objective(model, Min, sum(mat_y[p, w] * vec_costAdmin[p] +
#         (datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j), (datf_cost[:p] .== p))][1] *
#         mat5_z[i, j, p, t, w] + mat4_m[i, j, t, w] * 1000) * vec_prob[w] for w = 1: 10, p = 1: 10, i = 1: 4,
#         j = 1: 10, t = 1: 12))
#     @constraint(model, [p = 1:10, w = 1:10], mat_y[p, w] <= vec_numMax[p])
#     @constraint(model, [i = 1:4, j = 1:10, t = 1:12, w = 1:10], mat4_demand[i, j, t, w] -
#         sum([mat5_z[i, j, p, t, w] * vec_cap[p]] for p = 1:10)[1] <= mat4_m[i, j, t, w])
#     @constraint(model, [p = 1:10, t = 1: 12, w = 1: 10], sum(mat5_z[i, j, p, t, w] for i = 1: 4, j = 1: 10) <=
#         mat_y[p, w])
#     @constraint(model, [p = 1:10, w = 2:10], mat_y[p, 1] == mat_y[p, w])
#     solvePrintResult(model, mat_y, mat4_m)
# end

function optim_waitSeeSub(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat3_demand)
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_y[1:10] >= 0, Int)
    @variable(model, mat4_z[1:4, 1:10, 1:10, 1:12] >= 0, Int)
    @variable(model, mat3_m[1:4, 1:10, 1:12] >= 0)
    @objective(model, Min, sum(vec_y[p] * vec_costAdmin[p] for p = 1:10) +
        sum(datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j), (datf_cost[:p] .== p))][1] *
        mat4_z[i, j, p, t] + mat3_m[i, j, t] * 1000 for p = 1: 10, i = 1: 4, j = 1: 10, t = 1: 12))
    @constraint(model, [p = 1:10], vec_y[p] <= vec_numMax[p])
    @constraint(model, [i = 1:4, j = 1:10, t = 1:12], mat3_demand[i, j, t] -
        sum(mat4_z[i, j, p, t] * vec_cap[p] for p = 1:10) <= mat3_m[i, j, t])
    @constraint(model, [p = 1: 10, t = 1: 12], sum(mat4_z[i, j, p, t] for i = 1: 4, j = 1: 10) <= vec_y[p])
    solve(model)
    objResult = getobjectivevalue(model)
    vec_yResult = getvalue(vec_y)
    vec_mResult = getvalue(mat3_m)
    println("Result: obj = $(objResult)\n",
            "Result: vec_y = $(vec_yResult)\n",
            "Result: unfulfilled passenger demand = $(sum(vec_mResult))")
    return objResult
end

function optim_waitSee(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    ## 1,  Begin Optimization
    vec_objResult = zeros(10)
    for w = 1:10
        println("----- $(w)-th Sub-Problem -----")
        vec_objResult[w] = optim_waitSeeSub(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand[:,:,:,w])
    end
    objWeightedSum = sum(vec_prob[i] .* vec_objResult[i] for i = 1:10)
    println("----- Final Result -----\n",
            "Result: weighted sum of obj = $(objWeightedSum)")
end

function optim_robust(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    ## 1,  Begin Optimization
    model = Model(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    # model = Model(solver = GurobiSolver(Presolve=0))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_y[1: 10] >= 0, Int)
    @variable(model, mat5_z[1: 4, 1: 10, 1: 10, 1: 12, 1: 10] >= 0, Int)
    @variable(model, mat4_m[1: 4, 1: 10, 1: 12, 1: 10] >= 0)
    @variable(model, alpha >= 0)
    @objective(model, Min, sum(vec_y[p] * vec_costAdmin[p] for p = 1: 10) + alpha)
    @constraint(model, [w = 1: 10],
        sum((datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j), (datf_cost[:p] .== p))][1] *
        mat5_z[i, j, p, t, w] + mat4_m[i, j, t, w] * 1000) for p = 1: 10, i = 1: 4, j = 1: 10, t = 1: 12) <= alpha)
    @constraint(model, [p = 1: 10], vec_y[p] <= vec_numMax[p])
    @constraint(model, [i = 1: 4, j = 1: 10, t = 1: 12, w = 1: 10], mat4_demand[i, j, t, w] -
        sum([mat5_z[i, j, p, t, w] * vec_cap[p]] for p = 1: 10)[1] <= mat4_m[i, j, t, w])
    @constraint(model, [p = 1: 10, t = 1: 12, w = 1: 10], sum(mat5_z[i, j, p, t, w] for i = 1: 4, j = 1: 10) <=
        vec_y[p])
    solvePrintResult(model, vec_y, mat4_m)
end

function main()
    (vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob, mat3_demandExpect) = processData()
    println("################################################################################\n",
            "######################## 1/5, Deterministic Programming ########################\n",
            "################################################################################")
    vec_yResult = optim_determin(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat3_demandExpect)
    # vec_yResult = [0.0, 1.0, 2.0, 2.0, 0.0, 12.0, 5.0, 4.0, 8.0, 8.0]
    # optim_determinStochas(vec_yResult, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    # println("################################################################################\n",
    #         "########################## 2/5, Stochastic Programming #########################\n",
    #         "################################################################################")
    # optim_stochas(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    # println("################################################################################\n",
    #         "########################### 3/5, Wait-See Programming ##########################\n",
    #         "################################################################################")
    # optim_waitSee(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob)
    # println("################################################################################\n",
    #         "########################### 4/5, Robust Programming ############################\n",
    #         "################################################################################")
    # optim_robust(vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand)
    println("################################################################################\n",
            "################################### 5/5, End ###################################\n",
            "################################################################################")
end

main()


text:
# --------------
# Initialization
# --------------
using LLLplus
using Test
using LinearAlgebra
using Documenter

# --------------
# tests with small matrices
# --------------
@testset "tests with small matrices..." begin

    # Matrix from http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=442
    H =[1   9   1   2;
        1   8   8   3;
        7   4   5   1;
        2   6   7   1];
    B1=[2   3  -2  -4
        3  -1   2   1
        1   1   6  -4
        1   3  -1   3];
    T1=[0   0   1  -1
        0   1  -1   0
        0   0   0   1
        1  -3   3  -2];
    (B,T) = lll(H);
    @test sum(abs.(B-B1))==0
    @test sum(abs.(T-T1))==0
end

# --------------
# Run tests with random matrices
# --------------

@testset "random matrices" begin
    println("\nIn all the following tests, the first time includes the "*
            "JIT compilation; \nfor the second execution the "*
            "compilation is already done and the time\n"*
            "should be faster.\n")

    # Time LLL decomposition of a 1000x1000 real matrix with randn entries
    N = 1000;
    H = randn(N,N);
    println("Testing LLL on $(N)x$(N) real matrix...")
    @time (B,T,Q,R) = lll(H);
    @time (B,T,Q,R) = lll(H);

    # Time LLL, Seysen, VBLAST decompositions of a 10x10 complex matrix with
    # randn entries
    N = 10;
    println("\nTesting LLL on $(N)x$(N) complex matrix...")
    H = randn(N,N) + im*randn(N,N);
    @time (B,T,Q,R) = lll(H);
    @time (B,T,Q,R) = lll(H);
    println("Testing Seysen on same $(N)x$(N) complex matrix...")
    @time (B,T) = seysen(H);
    @time (B,T) = seysen(H);
    println("Testing VBLAST on same $(N)x$(N) complex matrix...")
    @time (W,P,B) = vblast(H);
    @time (W,P,B) = vblast(H);
    println("Testing Brun on real part of same $(N)x$(N) matrix...")
    @time (B,T) = brun(real.(H));
    @time (B,T) = brun(real.(H));
end

# Test CVP
@testset "CVP" begin
    N=500;
    println("\nTesting cvp on $(N)x$(N) BPSK system...")
    H=randn(N,N);
    Q,Rtmp=qr(H); R=UpperTriangular(Rtmp);
    u=Int.(rand([-1,1],N));
    y=H*u+rand(N)/100;
    uhat=cvp(Q'*y,R,Val(false),-1,1);
    errRate = sum(abs.(u-uhat))
    println("Error Rate is $(errRate). It should be zero or very small.\n")
end

# --------------
# doctests
# --------------
if VERSION >= v"1.1.0"
    println("\nRunning doctests...")
    DocMeta.setdocmeta!(LLLplus, :DocTestSetup,
                        :(using LLLplus, LinearAlgebra,Random); recursive=true)
    doctest(LLLplus)
end

max_stars_repo_name: christianpeel/LLLplus

max_issues_repo_name: christianpeel/LLLplus

max_forks_repo_name: christianpeel/LLLplus


text:
mutable struct Path
    group::Int
    θ::Vector{Float64}
    x::Vector{Int}
    t::Float64
    a::Vector{Float64}
    fly::Vector{Float64}
    ra::Vector{Int}
    dm::Float64
    w::Float64
    Path(group::Int,θ::Vector{Float64}, x::Vector{Int}, t::Float64, a::Vector{Float64}, fly::Vector{Float64}, ra::Vector{Int}, dm::Float64, w::Float64) = new(group, θ, x, t, a, fly, ra, dm, w)
end

Path(group::Int,θ::Vector{Float64}, x::Vector{Int}, t::Float64, ra::Vector{Int}) = Path(group, θ, copy(x), t, zeros(Float64, length(ra)), Float64.(ra), ra, 0.0, 1.0)
Path(group::Int,θ::Vector{Float64}, x::Vector{Int}, t::Float64, ra::Vector{Int}, dm::Float64) = Path(group, θ, copy(x), t, zeros(Float64, length(ra)), Float64.(ra), ra, dm, 1.0)

function Base.show(io::IO, p::Path)
    print("Path { \n group: ", p.group, " \n θ: ", p.θ, " \n x: ", p.x, ", \n t: ", p.t,  ", \n reaction firings: ", p.ra, ", \n objective value: ", p.dm, "}")
end

function getCostVal(p::Path)
    minimum(p.dm)
end
function getCostVal(p::Path, i)
    p.dm[i]
end

function getRates(p::Path)
    p.θ
end

function getGroup(p::Path)
    p.group
end

function getGroupCounts(ens::Vector{Path})
    counts(getGroup.(ens))
end

function overwritePath!(p1::Path, p2::Path)
    p1.group = p2.group
    p1.θ = copy(p2.θ)
    p1.x = copy(p2.x)
    p1.t = p2.t
    p1.a = copy(p2.a)
    p1.fly = copy(p2.fly)
    p1.ra = copy(p2.ra)
    p1.dm = copy(p2.dm)
end

function copyPath(p::Path)
    group = p.group
    θ = copy(p.θ)
    x = copy(p.x)
    t = p.t
    a = copy(p.a)
    fly = copy(p.fly)
    ra = copy(p.ra)
    dm = copy(p.dm)
    w = copy(p.w)
    Path(group, θ, x, t, a, fly, ra, dm, w)
end
function copyPath(p::Path, rw::Int)
    group = p.group
    θ = copy(p.θ)
    x = copy(p.x)
    t = p.t
    a = copy(p.a)
    fly = copy(p.fly)
    ra = copy(p.ra)
    dm = copy(p.dm)
    w = copy(p.w) / rw
    Path(group, θ, x, t, a, fly, ra, dm, w)
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
ENV["GKSwstype"] = 100 # workaround for gr segfault on GH actions
using LowLevelParticleFilters
import LowLevelParticleFilters.resample
using Test, Random, LinearAlgebra, Statistics, StaticArrays, Distributions, Plots
using MonteCarloMeasurements
gr(show=false)
Random.seed!(0)

@testset "LowLevelParticleFilters" begin
    @testset "logsumexp" begin
        w = randn(10); we = similar(w)
        ll = logsumexp!(w,we)
        @test sum(we) ≈ 1
        @test sum(exp.(w)) ≈ 1
        wc = copy(w)
        I = ones(10)
        logsumexp!(I,we)
        @test w ≈ wc.-log(sum(exp, wc))
        @test I ≈ fill(log(1/10),10)
        w = randn(10); we = similar(w)
        wc = copy(w)
        LowLevelParticleFilters.expnormalize!(we,w)
        @test sum(we) ≈ 1
        @test w ≈ wc atol=1e-15
        LowLevelParticleFilters.expnormalize!(w)
        @test sum(w) ≈ 1
    end

    @testset "weigthed_mean" begin
        x = [randn(3) for i = 1:10000]
        w = ones(10000); we = similar(w)
        logsumexp!(w,we)
        @test sum(abs, weigthed_mean(x,we)) < 0.06
    end

    @testset "resample" begin
        s = PFstate(10)
        @test effective_particles(s.we) ≈ 10
        logsumexp!(s.w, s.we)
        @test resample(s.we) ≈ 1:10

        w = [1.,1,1,2,2,2,3,3,3]; we=similar(w)
        logsumexp!(w,we)
        @test we |> resample |> sum >= 56
        @test length(resample(we)) == length(we)
        for i = 1:10
            w = randn(100); we = randn(100)
            logsumexp!(w,we)
            j = resample(we)
            @test maximum(j) <= 100
            @test minimum(j) >= 1
        end
    end

    @testset "static distributions" begin
        x = @SVector ones(2)
        d = MvNormal(2,2)
        @test logpdf(d,x) == logpdf(d,Vector(x))
        d = Product([Normal(0,2), Normal(0,2)])
        @test logpdf(d,x) == logpdf(d,Vector(x))

        dt = LowLevelParticleFilters.TupleProduct((Normal(0,2), Normal(0,2)))
        @test logpdf(dt,x) == logpdf(dt,Vector(x)) == logpdf(d,x)
        @test_nowarn rand(dt)
        @test var(dt) == var(d)
        @test cov(dt) == cov(d)
        @test entropy(dt) == entropy(d)
        @test rand(dt) isa Vector{Float64}

    end

    @testset "End to end" begin
        eye(n) = Matrix{Float64}(I,n,n)
        n = 2 # Dinemsion of state
        m = 2 # Dinemsion of input
        p = 2 # Dinemsion of measurements

        dg = MvNormal(p,1.0)          # Dynamics noise Distribution
        df = MvNormal(n,0.1)          # Measurement noise Distribution
        d0 = MvNormal(randn(n),2.0)   # Initial state Distribution

        # Define random linenar state-space system
        Tr = randn(n,n)
        A = SMatrix{n,n}([0.99 0.1; 0 0.2])
        B = @SMatrix randn(n,m)
        C = SMatrix{p,p}(eye(p))
        # C = SMatrix{p,n}([1 1])

        dynamics(x,u,t) = A*x .+ B*u
        measurement(x,u,t) = C*x

        N     = 1000 # Number of particles
        T     = 200 # Number of time steps
        M     = 100 # Number of smoothed backwards trajectories
        pf    = ParticleFilter(N, dynamics, measurement, df, dg, d0)
        pfa   = AuxiliaryParticleFilter(N, dynamics, measurement, df, dg, d0)
        @test !shouldresample(pf)
        @test !shouldresample(pfa)
        du    = MvNormal(2,1) # Control input distribution
        xp,up,yp = LowLevelParticleFilters.simulate(pf,T,du,100)
        @test xp[1] isa MonteCarloMeasurements.Particles{Float64,100}
        @test size(xp) == (T,n)
        x,u,y = LowLevelParticleFilters.simulate(pf,T,du)

        xf, wf, wef, ll = forward_trajectory(pf, u, y)
        parts = Particles(xf,wef)
        @test size(parts) == (T,n)
        @test length(parts[1].particles) == N

        xm = reduce(hcat,x)
        tosvec(y) = reinterpret(SVector{length(y[1]),Float64}, reduce(hcat,y))[:] |> copy
        x,u,y = tosvec.((x,u,y))

        xpf = mean_trajectory(pf,u,y)[1]

        xb,ll = smooth(pf, M, u, y) # Sample smooting trajectories
        @test size(xb) == (M,T)
        xbm = smoothed_mean(xb)     # Calculate the mean of smoothing trajectories
        @test mean(abs2, xm - xbm) < 5

        xb,ll = smooth(pfa, M, u, y)
        xbma = smoothed_mean(xb)
        @test mean(abs2, xm - xbma) < 5

        @show mean(abs2, xm - xbm)
        @show mean(abs2, xm - xbma)

        xbc = smoothed_cov(xb)      # And covariance
        @test all(tr(C) < 2 for C in xbc)
        maximum(tr(C) for C in xbc)
        xbt = smoothed_trajs(xb)    # Can't remember what this does

        kf     = KalmanFilter(A, B, C, 0, 0.01eye(n), eye(p), d0)
        # x,u,y = simulate(kf,T,du)
        xf,xt,R,Rt,ll = forward_trajectory(kf, u, y)
        xT,R,lls = smooth(kf, u, y)

        # @test_broken mean(abs2, xm) > mean(abs2, xm - reduce(hcat,xf)) > mean(abs2, xm - reduce(hcat,xt)) > mean(abs2, xm - reduce(hcat,xT))
        # @test_broken mean(abs2, xm) > mean(abs2, xm - reduce(hcat,xpf))  > mean(abs2, xm - reduce(hcat,xt)) > mean(abs2, xm - reduce(hcat,xT))
        # @test_broken mean(abs2, xm) > mean(abs2, xm - reduce(hcat,xpf))  > mean(abs2, xm - xbm) > mean(abs2, xm - reduce(hcat,xT))
        # plot(xm', layout=2)
        # plot!(reduce(hcat,xf)')
        # plot!(reduce(hcat,xt)')
        # plot!(reduce(hcat,xT)')

        svec = exp10.(LinRange(-1,2,22))
        llspf = map(svec) do s
            df = MvNormal(n,s)
            pfs = ParticleFilter(N, dynamics, measurement, df, dg, d0)
            loglik(pfs,u,y)
        end

        llspfa = map(svec) do s
            df = MvNormal(n,s)
            pfs = AuxiliaryParticleFilter(N, dynamics, measurement, df, dg, d0)
            loglik(pfs,u,y)
        end

        llskf = map(svec) do s
            kfs = KalmanFilter(A, B, C, 0, s^2*eye(n), eye(p), d0)
            loglik(kfs,u,y)
        end
        # plot(svec, [llspf llspfa llskf], xscale=:log10, lab=["PF" "APF" "KF"])

        @testset "Metropolis" begin
            N = 1000
            function filter_from_parameters(θ,pf=nothing)
                pf === nothing && (return ParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0))
                ParticleFilter(pf.state, dynamics, measurement, df,dg, d0)
            end
            # The call to `exp` on the parameters is so that we can define log-normal priors
            priors = [Normal(1,2),Normal(1,2)]
            ll     = log_likelihood_fun(filter_from_parameters,priors,u,y)
            θ₀ = log.([1.,1.]) # Starting point
            # We also need to define a function that suggests a new point from the "proposal distribution". This can be pretty much anything, but it has to be symmetric since I was lazy and simplified an equation.
            draw = θ -> θ .+ rand(MvNormal(0.1ones(2)))
            burnin = 200
            theta, lls = metropolis(ll, 20, θ₀, draw)

        end

        @testset "example_lineargaussian" begin
            include("../src/example_lineargaussian.jl")
        end

    end

    @testset "debugplot" begin
        eye(n) = Matrix{Float64}(I,n,n)
        n = 2 # Dinemsion of state
        m = 1 # Dinemsion of input
        p = 1 # Dinemsion of measurements

        dg = MvNormal(p,1.0)          # Dynamics noise Distribution
        df = MvNormal(n,0.1)          # Measurement noise Distribution
        d0 = MvNormal(randn(n),2.0)   # Initial state Distribution

        # Define random linenar state-space system
        Tr = randn(n,n)
        A = SMatrix{n,n}([0.99 0.1; 0 0.2])
        B = @SMatrix [0;1]
        C = @SMatrix [1 0]
        # C = SMatrix{p,n}([1 1])

        dynamics(x,u,t) = A*x .+ B*u
        measurement(x,u,t) = C*x

        N     = 100 # Number of particles
        T     = 3   # Number of time steps
        pf    = ParticleFilter(N, dynamics, measurement, df, dg, d0)
        pfa   = AuxiliaryParticleFilter(N, dynamics, measurement, df, dg, d0)
        du    = MvNormal(1,1) # Control input distribution
        x,u,y = LowLevelParticleFilters.simulate(pf,T,du)
        ##

        debugplot(pf,u,y, runall=true, xreal=x, density=true)
        debugplot(pf,u,y, runall=true, xreal=x, density=false)
        debugplot(pf,u,y, runall=true, xreal=x, leftonly=false)
        debugplot(pf,u,y, runall=true, xreal=x, leftonly=false, density=false)
        # debugplot(pfa,x,u,y, runall=true)

        # commandplot(pf,x,u,y)
        # commandplot(pfa,x,u,y)
    end

@testset "UKF" begin

    eye(n) = Matrix{Float64}(I,n,n)
    n = 2 # Dinemsion of state
    m = 2 # Dinemsion of input
    p = 2 # Dinemsion of measurements

    dg = MvNormal(p,1.0)          # Dynamics noise Distribution
    df = MvNormal(n,0.1)          # Measurement noise Distribution
    d0 = MvNormal(randn(n),2.0)   # Initial state Distribution
    du = MvNormal(2,1) # Control input distribution

    # Define random linenar state-space system
    Tr = randn(n,n)
    A = SA[0.99 0.1; 0 0.2]
    B = @SMatrix randn(n,m)
    C = SMatrix{p,p}(eye(p))
    # C = SMatrix{p,n}([1 1])

    dynamics(x,u,t) = A*x .+ B*u
    measurement(x,u,t) = C*x

    T     = 200 # Number of time steps
    kf   = KalmanFilter(A, B, C, 0, eye(n), eye(p), d0)
    ukf  = UnscentedKalmanFilter(dynamics, measurement, eye(n), eye(p), d0)
    x,u,y = LowLevelParticleFilters.simulate(kf,T,du) # Simuate trajectory using the model in the filter
    @test_nowarn LowLevelParticleFilters.simulate(ukf,T,du)
    tosvec(y) = reinterpret(SVector{length(y[1]),Float64}, reduce(hcat,y))[:] |> copy
    x,u,y = tosvec.((x,u,y))

    reskf = forward_trajectory(kf, u, y) # filtered, prediction, pred
    resukf = forward_trajectory(ukf, u, y)

    norm(mean(x .- reskf[1]))
    norm(mean(x .- resukf[1]))

    norm(mean(x .- reskf[2]))
    norm(mean(x .- resukf[2]))
    @test norm(mean(x .- reskf[2])) < norm(mean(x .- reskf[1])) # Filtered should be better than prediction
    @test norm(mean(x .- resukf[2])) < norm(mean(x .- resukf[1]))
    @test norm(mean(x .- reskf[2])) ≈ norm(mean(x .- resukf[2])) atol=5e-2
    # @test norm(mean(x .- reskf[2])) < norm(mean(x .- resukf[2]))  # KF should be better than UKF
    # @test norm(mean(x .- reskf[1])) < norm(mean(x .- resukf[1]))  # KF should be better than UKF
    @test norm(mean(x .- reskf[2])) < 0.2
end

@testset "Advanced filters" begin

    eye(n) = Matrix{Float64}(I,n,n)
    n = 2 # Dinemsion of state
    m = 2 # Dinemsion of input
    p = 2 # Dinemsion of measurements

    dg = MvNormal(p,1.0)          # Dynamics noise Distribution
    df = MvNormal(n,0.1)          # Measurement noise Distribution
    d0 = MvNormal(randn(n),2.0)   # Initial state Distribution
    du = MvNormal(2,1)            # Control input distribution

    # Define random linenar state-space system
    A = SMatrix{n,n}([0.99 0.1; 0 0.2])
    B = @SMatrix randn(n,m)
    C = SMatrix{p,p}(eye(p))
    # C = SMatrix{p,n}([1 1])

    dynamics(x,u,t,noise=false) = A*x .+ B*u + noise*rand(df)
    measurement(x,u,t,noise=false) = C*x .+ noise*rand(dg)
    measurement_likelihood(x,u,y,t) = logpdf(dg, measurement(x,u,t)-y)
    tosvec(y) = reinterpret(SVector{length(y[1]),Float64}, reduce(hcat,y))[:] |> copy

    T     = 200 # Number of time steps
    N     = 500
    Random.seed!(0)
    apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0)
    x,u,y = LowLevelParticleFilters.simulate(apf,T,du) # Simuate trajectory using the model in the filter
    x,u,y = tosvec.((x,u,y))
    @time resapf,ll = mean_trajectory(apf, u, y)
    sf = SigmaFilter(N,dynamics, measurement, measurement_likelihood, df, d0)
    @test_nowarn LowLevelParticleFilters.simulate(sf,T,du)
    @time ressf,ll = mean_trajectory(sf, u, y)

    norm(mean(x))
    norm(mean(x .- resapf))
    norm(mean(x .- ressf))

    @test norm(mean(x .- ressf)) ≈ norm(mean(x .- resapf)) atol=2e-1
    @test_skip norm(mean(x .- ressf)) < norm(mean(x .- resapf))  # SF should be better than PF since we now covs are Gaussian, this test is not robust enough to include
end

@testset "ekf" begin
    @info "Testing ekf"
    include("test_ekf.jl")
end

end


text:
"""
    is_same_except(m1::MLJType, m2::MLJType, exceptions::Symbol...)

Returns `true` only the following conditions all hold:

- `m1` and `m2` have the same type.

- `m1` and `m2` have the same undefined fields.

- Corresponding fields agree, or are listed as
  `exceptions`, or have `AbstractRNG` as values (one or both)

Here "agree" is in the sense of "==", unless the objects are
themselves of `MLJType`, in which case agreement is in the sense of
`is_same_except` with no exceptions allowed.

Note that Base.== is overloaded such that `m1 == m2` if and only if
`is_same_except(m1, m2)`.

"""
is_same_except(x1, x2) = ==(x1, x2)
function is_same_except(m1::M1, m2::M2,
            exceptions::Symbol...) where {M1<:MLJType,M2<:MLJType}
    if typeof(m1) != typeof(m2)
        return false
    end
    defined1 = filter(fieldnames(M1)|>collect) do fld
        isdefined(m1, fld) && !(fld in exceptions)
    end
    defined2 = filter(fieldnames(M1)|>collect) do fld
        isdefined(m2, fld) && !(fld in exceptions)
    end
    if defined1 != defined2
        return false
    end
    same_values = true
    for fld in defined1
        same_values = same_values &&
            (is_same_except(getfield(m1, fld), getfield(m2, fld)) ||
             getfield(m1, fld) isa AbstractRNG ||
             getfield(m2, fld) isa AbstractRNG)
    end
    return same_values
end

import Base.==

==(m1::M1, m2::M2) where {M1<:MLJType,M2<:MLJType} = is_same_except(m1, m2)

# for using `replace` or `replace!` on collections of MLJType objects
# (eg, Model objects in a learning network) we need a stricter
# equality:
MLJBase.isequal(m1::MLJType, m2::MLJType) = (m1 === m2)

## TODO: Do we need to overload hash here?

# Note: To prevent julia crash, it seems we mysteriously need to
# annotate the type of itr:
function special_in(x, itr)::Union{Bool,Missing}
    for y in itr
        ismissing(y) && return missing
        y === x && return true
    end
    return false
end
Base.in(x::MLJType, itr::Set) = special_in(x, itr)
Base.in(x::MLJType, itr::AbstractVector) = special_in(x, itr)
Base.in(x::MLJType, itr::NTuple) = special_in(x, itr)


text:
module Runtests
using Unregistered
using Test

@test true

end


text:
using DynamicGrids, Dates, DimensionalData, Setfield, Unitful, Test
using Unitful: d
using DynamicGrids: SimData, Extent, _calc_auxframe, _cyclic_index
const DG = DynamicGrids

@testset "Aux" begin
    @testset "sequence cycling" begin
        @test _cyclic_index(-4, 2) == 2
        @test _cyclic_index(-3, 2) == 1
        @test _cyclic_index(-2, 2) == 2
        @test _cyclic_index(-1, 2) == 1
        @test _cyclic_index(0, 2) == 2
        @test _cyclic_index(1, 2) == 1
        @test _cyclic_index(2, 2) == 2
        @test _cyclic_index(3, 2) == 1
        @test _cyclic_index(4, 2) == 2
        @test _cyclic_index(20, 10) == 10
        @test _cyclic_index(21, 10) == 1
        @test _cyclic_index(27, 10) == 7
    end

    @testset "aux sequence" begin
        a = cat([0.1 0.2; 0.3 0.4], [1.1 1.2; 1.3 1.4], [2.1 2.2; 2.3 2.4]; dims=3)

        @testset "the correct frame is calculated for aux data" begin
            dimz = X(1:2), Y(1:2), Ti(15d:5d:25d)
            seq = DimensionalArray(a, dimz)
            init = zero(seq[Ti(1)])
            sd = SimData(Extent(init=init, aux=(seq=seq,), tspan=1d:1d:100d), Ruleset())
            @test DynamicGrids.boundscheck_aux(sd, Aux{:seq}()) == true
            tests = (1, 1), (4, 1), (5, 2), (6, 2), (9, 2), (10, 3), (11, 3), (14, 3), (15, 1),
                    (19, 1), (20, 2), (25, 3), (29, 3), (30, 1), (34, 1), (35, 2)
            for (f, ref_af) in tests
                @set! sd.currentframe = f
                af = _calc_auxframe(sd).seq
                @test af == ref_af
            end
        end

        @testset "boundscheck_aux" begin
            seq = DimensionalArray(a, (X(1:2), Y(1:2), Ti(15d:5d:25d)))
            bigseq = DimensionalArray(zeros(5, 2, 3), (X(1:5), Y(1:2), Ti(15d:5d:25d)))
            init = zero(seq[Ti(1)])
            sd = SimData(Extent(init=init, aux=(seq=seq, bigseq=bigseq), tspan=1d:1d:100d), Ruleset())
            @test DynamicGrids.boundscheck_aux(sd, Aux{:seq}()) == true
            @test_throws ErrorException DynamicGrids.boundscheck_aux(sd, Aux{:bigseq}())
            @test_throws ErrorException DynamicGrids.boundscheck_aux(sd, Aux{:missingseq}())
        end

        @testset "correct values are returned by get" begin
            dimz = X(1:2), Y(1:2), Ti(Date(2001, 1, 15):Day(5):Date(2001, 1, 25))
            seq = DimensionalArray(a, dimz)
            init = zero(seq[Ti(1)])
            data = SimData(Extent(init=init, aux=(seq=seq,), tspan=Date(2001):Day(1):Date(2001, 3)), Ruleset())
            data1 = DG._updatetime(data, 1)
            @test data1.auxframe == (seq = 1,)
            @test get(data1, Aux(:seq), (1, 1)) == 0.1
            data2 = DG._updatetime(data, 5)
            @test data2.auxframe == (seq = 2,)
            @test get(data2, Aux(:seq), 1, 1) == 1.1
            data3 = DG._updatetime(data, 10)
            @test data3.auxframe == (seq = 3,)
            @test get(data3, Aux(:seq), (1, 1)) == 2.1
            data4 = DG._updatetime(data, 15)
            @test data4.auxframe == (seq = 1,)
            @test get(data4, Aux(:seq), CartesianIndex(1, 1)) == 0.1
        end

        @testset "errors" begin
            output = ArrayOutput(zeros(3, 3); tspan=1:3)
            @test_throws ArgumentError  DynamicGrids.aux(output, Aux{:somekey}())
        end
    end
end

# Use copyto to test all parametersources, as well as testing CopyTo itself
@testset "CopyTo" begin
    @testset "Copy construction" begin
        rule = CopyTo(7)
        rule2 = @set rule.from = Aux{:a}()
        @test rule2.from == Aux{:a}()
    end

    init = [0 0]

    @testset "CopyTo from value" begin
        ruleset = Ruleset(CopyTo(7))
        output = ArrayOutput(init; tspan=1d:1d:3d)
        sim!(output, ruleset)
        @test output == [[0 0], [7 7], [7 7]]
    end

    @testset "CopyTo from Grid" begin
        @test CopyTo(Aux(:l)) === CopyTo(; from=Aux(:l))
        @test CopyTo{:a}(; from=Aux(:l)) === CopyTo{:a}(Aux(:l))
        ruleset = Ruleset(CopyTo(Aux(:l)))
        output = ArrayOutput(init; tspan=1:3, aux=(l=[3 4],))
        sim!(output, ruleset)
        @test output == [[0 0], [3 4], [3 4]]

        da = DimArray(cat([1 2], [3 4]; dims=3) , (X(), Y(), Ti(4d:1d:5d)))
        output = ArrayOutput(init; tspan=1d:1d:3d, aux=(l=da,))
        sim!(output, ruleset)
        @test output == [[0 0], [1 2], [3 4]]
    end

    @testset "CopyTo from Grid" begin
        ruleset = Ruleset(Cell{:s,:s}((d, x, I) -> x + 1), CopyTo{:d}(from=Grid(:s)))
        output = ArrayOutput((s=[1 3], d=[0 0],); tspan=1d:1d:3d)
        sim!(output, ruleset)
        @test output == [(s=[1 3], d=[0 0]), (s=[2 4], d=[2 4]), (s=[3 5], d=[3 5])]

        ruleset = Ruleset(Cell{:s,:s}((d, x, I) -> x + 1), CopyTo{Tuple{:d1,:d2}}(from=Grid{:s}()))
        output = ArrayOutput((s=[1 3], d1=[0 0], d2=[-1 -1],); tspan=1d:1d:3d)
        sim!(output, ruleset)
        @test output == [(s=[1 3], d1=[0 0], d2=[-1 -1]),
                         (s=[2 4], d1=[2 4], d2=[2 4]),
                         (s=[3 5], d1=[3 5], d2=[3 5])]
    end

    @testset "CopyTo from Delay" begin
        ruleset = Ruleset(Cell{:s,:s}((d, x, I) -> x + 1), CopyTo{:d}(from=Delay{:s}(1d)))
        @test DynamicGrids.hasdelay(rules(ruleset)) == true
        output = ArrayOutput((s=[1 3], d=[0 0],); tspan=1d:1d:4d)
        sim!(output, ruleset)
        @test output == [
            (s=[1 3], d=[0 0]),
            (s=[2 4], d=[1 3]),
            (s=[3 5], d=[2 4]),
            (s=[4 6], d=[3 5])
        ]

        ruleset = Ruleset(Cell{:s,:s}((d, x, I) -> x + 1), CopyTo{:d}(from=Delay{:s}(Month(2))))
        @test DynamicGrids.hasdelay(rules(ruleset)) == true
        output = ArrayOutput((s=[1 3], d=[0 0]); tspan=Date(2001):Month(1):Date(2001, 6))
        sim!(output, ruleset)
        @test output == [
            (s=[1 3],   d=[0 0]),
            (s=[2 4],   d=[1 3]),
            (s=[3 5],   d=[1 3]),
            (s=[4 6],   d=[2 4]),
            (s=[5 7],   d=[3 5]),
            (s=[6 8],   d=[4 6]),
        ]

    end

    @testset "CopyTo from Lag" begin
        ruleset = Ruleset(Cell{:s,:s}((d, x, I) -> x + 1), CopyTo{:d}(from=Lag{:s}(1)))
        @test DynamicGrids.hasdelay(rules(ruleset)) == true
        output = ArrayOutput((s=[1 3], d=[0 0],); tspan=1d:1d:4d)
        sim!(output, ruleset)
        @test output == [
            (s=[1 3], d=[0 0]),
            (s=[2 4], d=[1 3]),
            (s=[3 5], d=[2 4]),
            (s=[4 6], d=[3 5])
        ]

    end
end


text:
using Pkg
Pkg.activate("TestPackage1")

#%%
using TestPackage1

TestPackage1.greet()

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
using Images
using ImageView
using VideoIO
using PlutoUI
using ImageTransformations

north = load("/home/shefi/Pictures/flat_earth_north.png");
south = load("/home/shefi/Pictures/flat_earth_south.png");

# Scale the images
scale = 384  # 128*3
scale = 1024  # 128*6
north_s = imresize(north, scale, scale);
south_s = imresize(south, scale, scale);

println("North: $(size(north))  small version $(size(north_s))")
println("South: $(size(south))  small version $(size(south_s))")

function az_to_lat(az::Int)
    az = az % 360
    if az < 180
        "E $(Int(round(az, digits=0)))°"
    else
        "W $(Int(round(180-az, digits=0)))°"
    end
end

function flat_vertical_map(az::Int)
    rad = az * (2π/360);

    # keeps the axes - this is great for a cicular image
    north_rot = imrotate(north_s, rad + π/2, axes(north_s));
    south_rot = imrotate(south_s, -rad + π/2, axes(south_s));

    flat = vcat(north_rot,south_rot)
end

function flat_horizontal_map(az::Int, north_s::Array{T,2}, south_s::Array{T,2}) where T<:Colorant
    rad = az * (2π/360);

    # keeps the axes - this is great for a cicular image
    north_rot = imrotate(north_s, rad, axes(north_s));
    south_rot = imrotate(south_s, -rad, axes(south_s));

    flat = hcat(north_rot,south_rot)
end

function quad_flat(az::Int)
    rad = az * (2π/360);
    tl = imrotate(north_s, rad, axes(north_s));
    tr = imrotate(south_s, -rad, axes(south_s));
    bl = imrotate(south_s, -rad + π, axes(south_s));
    br = imrotate(north_s, rad + π, axes(north_s));

    flat = hcat(vcat(tl, bl), vcat(tr, br))
end

frames = []
for az in 100:1:520
    flat_earth = flat_vertical_map(az)
    frame = convert.(RGB, flat_earth)
    # annotate!(1.,1.,az_to_lat(az))
    push!(frames, frame)
end

vid_file = "/tmp/test.mp4"

# Using ffmpeg - did not work yet
# open(`ffmpeg -f rawvideo -pix_fmt gray -s:v 512x512 -r 70 -i pipe:0 $(vid_file)`, "w") do out
#     for frame in frames
#         # img = convert(Image{Gray{U8}}, frame)
#         # img = convert(Image{fmt}, frame)
#         img = frame
#         write(out, reinterpret(UInt8, data(img)))
#      end
# end

# using VideoIO
props = [:priv_data => ("crf"=>"22","preset"=>"medium")]
# encodevideo("video.mp4",frames,framerate=15,AVCodecContextProperties=props)
encodevideo("flat_map.mp4", frames, framerate=25, AVCodecContextProperties=props)


text:
module LibXLS

# Load libreadstat from our deps.jl
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibXLS not installed properly, run Pkg.build(\"LibXLS\"), restart Julia and try again")
end
include(depsjl_path)

include("c.jl")
include("types.jl")
include("workbook.jl")
include("worksheet.jl")

end


text:
"""
Type: StoppingMeta
Methods: no

Common stopping criteria for algorithms:
    - absolute and relative tolerance
    - threshold for unboundedness
    - time limit to let the algorithm run
    - maximum number of function (and derivatives) evaluations

It's a mutable struct therefore we can modified elements of a StoppingMeta.
	- The nb_of_stop is incremented everytime stop! or update_and_stop! is called
	- The optimality0 is modified once at the beginning of the algorithm (start!)
    - The start_time is modified once at the beginning of the algorithm (start!)
      if not precised before.
	- The different status: fail_sub_pb, unbounded, unbounded_pb, tired, stalled,
      iteration_limit, resources, optimal, main_pb, domainerror, suboptimal, infeasible

Note: fail_sub_pb, suboptimal, and infeasible are modified by the algorithm.
"""
mutable struct StoppingMeta <: AbstractStoppingMeta

 # problem tolerances
 atol                :: Number # absolute tolerance
 rtol                :: Number # relative tolerance
 optimality0         :: Number # value of the optimality residual at starting point
 tol_check           :: Function #function of atol, rtol and optimality0
                                 #by default: tol_check = max(atol, rtol * optimality0)
                                 #other example: atol + rtol * optimality0

 unbounded_threshold :: Number # beyond this value, the problem is declared unbounded
 unbounded_x         :: Number # beyond this value, x is unbounded

 # fine grain control on ressources
 max_f               :: Int    # max function evaluations allowed

 # global control on ressources
 max_eval            :: Int    # max evaluations (f+g+H+Hv) allowed
 max_iter            :: Int    # max iterations allowed
 max_time            :: Number # max elapsed time allowed

 #intern Counters
 nb_of_stop :: Int
 #intern start_time
 start_time :: Float64

 # stopping properties status of the problem)
 fail_sub_pb         :: Bool
 unbounded           :: Bool
 unbounded_pb        :: Bool
 tired               :: Bool
 stalled             :: Bool
 iteration_limit     :: Bool
 resources           :: Bool
 optimal             :: Bool
 infeasible          :: Bool
 main_pb             :: Bool
 domainerror         :: Bool
 suboptimal          :: Bool

 function StoppingMeta(;atol                :: Number   = 1.0e-6,
                        rtol                :: Number   = 1.0e-15,
                        optimality0         :: Number   = 1.0,
                        tol_check           :: Function = (atol,rtol,opt0) -> max(atol,rtol*opt0),
                        unbounded_threshold :: Number   = 1.0e50,
                        unbounded_x         :: Number   = 1.0e50,
                        max_f               :: Int      = typemax(Int),
                        max_eval            :: Int      = 20000,
                        max_iter            :: Int      = 5000,
                        max_time            :: Number   = 300.0,
                        start_time          :: Float64  = NaN,
                        kwargs...)

   fail_sub_pb = false

   try
       tol_check(1,1,1)
   catch
       throw("tol_check must have 3 arguments")
   end

   unbounded       = false
   unbounded_pb    = false
   tired           = false
   stalled         = false
   iteration_limit = false
   resources       = false
   optimal         = false
   infeasible      = false
   main_pb         = false
   domainerror     = false
   suboptimal      = false

   nb_of_stop = 0

   return new(atol, rtol, optimality0, tol_check, unbounded_threshold, unbounded_x,
              max_f, max_eval, max_iter, max_time, nb_of_stop, start_time,
              fail_sub_pb, unbounded, unbounded_pb, tired, stalled,
              iteration_limit, resources, optimal, infeasible, main_pb,
              domainerror, suboptimal)
 end
end


text:
module EbitODEServer

export start_ode_server
using Sockets
using Main.EbitSolver
using Main.EbitODEMessages
using ProtoBuf

function create_proto_msg(buffer)
    @debug "Creating EbitODEMessages.Message() from buffer"
    return readproto(PipeBuffer(buffer), EbitODEMessages.Message())
end

function create_proto_err_msg(msg::AbstractString)
    return EbitODEMessages.Message(msg_type=EbitODEMessages.MessageType.Error,
                                   err = EbitODEMessages.ErrorEncountered(msg=msg))
end

function read_msg_from_stream(socket::IO)
    read_no_bytes = read(socket, UInt32)
    @debug "Read message size: $read_no_bytes bytes"
    buffer = read(socket, read_no_bytes)
    msg = create_proto_msg(buffer)
    @debug "Read message buffer and created Message"
    return msg
end

function process(msg::EbitODEMessages.Message, report_progress)
    try
        if msg.msg_type == EbitODEMessages.MessageType.SolveODE
            ret = EbitSolver.solve_ode(msg.ode_problem, report_progress)
            return ret
        end
        return create_proto_err_msg("Unable to handle msg type: $msg.msg_type")
    catch e
        @warn "Catched exception" e
        return create_proto_err_msg("exception: $e")
    end
end

function send_proto_msg_to_stream(message, socket)
    iob = PipeBuffer()
    @debug "Created buffer"
    len = writeproto(iob, message)
    @debug "Writing return size: $len"
    write(socket, convert(UInt32, len))
    @debug "Writing message"
    write(socket, iob)
    @debug "Written message"
end

server = Sockets.TCPServer()

function ode_bind_server(port)
    global server
    if ( server.status == Base.StatusActive
         || server.status == Base.StatusOpen
         || server.status == Base.StatusConnecting )
        close(server)
    end
    new_server = Sockets.TCPServer()

    !bind(new_server, IPv4(UInt(0)), port) && error("cannot bind to port; may already be in use or access denied")
    global server = new_server
end

progress_message = EbitODEMessages.Message(msg_type=EbitODEMessages.MessageType.StatusUpdate,
                                           status = EbitODEMessages.Status.SolivingODEInProgress,
                                           progress = EbitODEMessages.Progress(time = 0.0))

function send_progress(socket, t, abort)
    progress_message.progress.time = t
    send_proto_msg_to_stream(progress_message, socket)
    return abort()
end

@noinline function start_ode_server(port=2000)
    ode_bind_server(port)
    listen(server)
    @info "Created server" port
    try
        while true
            socket = accept(server)
            @info "Accepting connection on port: $port"
            @async while isopen(socket)
                try
                    abort = false
                    msg = read_msg_from_stream(socket)

                    @async while isopen(socket)
                        msg = read_msg_from_stream(socket)
                        if (msg.msg_type == EbitODEMessages.MessageType.StopServer)
                            abort = true
                        end
                    end

                    ret_val = process(msg, t -> send_progress(socket, t, () -> abort))
                    send_proto_msg_to_stream(ret_val, socket)
                catch e
                    if isa(e, EOFError)
                        nothing
                    else
                        @warn "Caught exception" e
                    end
                end
            end
        end
    catch e
        @debug "Caught exception" e
    end
    @debug "Closing EbitODEServer on port $port"
end

function async_server_start(port=2000)
    @async start_ode_server(port)
    return server
end

function restart(port=2000)
    include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitSolver.jl")
    include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitODEServer.jl")
    close(server)
    EbitODEServer.async_server_start(port)
end

end


text:
function quat(q0,q1,q2,q3)

  q_length=sqrt(q0*q0+q1*q1+q2*q2+q3*q3);

  q0=q0/q_length;
  q1=q1/q_length;
  q2=q2/q_length;
  q3=q3/q_length;

  phi=atan2(2*(q0*q1+q2*q3),1-2*(q1*q1+q2*q2));
  theta=asin(2*(q0*q2-q3*q1));
  psii=atan2(2*(q0*q3+q1*q2),1-2*(q2*q2+q3*q3));

  return phi, theta, psii
end

phi, theta, psi=quat(0.99619,0,0,0.08716)
psi*180/pi

function Quaternion_to_EulerianAngle(x, y, z, w)
  ysqr = y*y

  t0 = +2.0 * (w * x + y*z)
  t1 = +1.0 - 2.0 * (x*x + ysqr)
  X = atan2(t0, t1)*180/pi

  t2 = +2.0 * (w*y - z*x)
  t2 = ((t2 > 1.0) ? 1.0 : t2);
  t2 = ((t2 < -1.0) ? -1.0 : t2);
  Y = asin(t2)*180/pi

  t3 = +2.0 * (w * z + x*y)
  t4 = +1.0 - 2.0 * (ysqr + z*z)
  Z = atan2(t3, t4)*180/pi

  return X, Y, Z
end

type	Quaternion
  x
  y
  z
  w
end

function

function toQuaternion(roll, pitch, yaw)

	t0 = cos(yaw/2)
	t1 = sin(yaw/2)
	t2 = cos(roll/2)
	t3 = sin(roll/2)
	t4 = cos(pitch/2)
	t5 = sin(pitch/2)

  w = t0*t2*t4 + t1*t3*t5
	x = t0*t3*t4 - t1*t2*t5
	y = t0*t2*t5 + t1*t3*t4
	z = t1*t2*t4 - t0*t3*t5

  q = Quaternion(x,y,z,w)

	return q
end

roll = 10
pitch = 0
yaw = 0
toQuaternion(roll, pitch, yaw)
phi, theta, psi=Quaternion_to_EulerianAngle(0.99619,0,0,0.08716)

using PyCall

@pyimport tf.transformations as tf

roll = 0
pitch = 0
yaw = pi/2

Q = tf.quaternion_from_euler(roll, pitch, yaw)

roll, pitch, yaw = tf.euler_from_quaternion(Q)


text:
using Test, DiffEqOperators, SpecialFunctions

@testset "None BC DerivativeOperator" begin
    N = 10
    d_order = 1
    approx_order = 2
    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:3] ≈ [-1.5, 2.0, -0.5]
    for row = 2:N-1
        @test Amat[row,row-1:row+1] ≈ [-0.5, 0.0, 0.5]
    end
    @test Amat[end,end-2:end] ≈ [0.5, -2.0, 1.5]

    d_order = 2
    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:4] ≈ [2.0, -5.0, 4.0, -1.0]
    for row = 2:N-1
        @test Amat[row,row-1:row+1] ≈ [1.0, -2.0, 1.0]
    end
    @test Amat[end,end-3:end] ≈ [-1.0, 4.0, -5.0, 2.0]

    N = 10
    d_order = 1
    approx_order = 3
    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:4] ≈ [-11/6, 3, -3/2, 1/3]
    for row = 3:N-2
        @test Amat[row,row-1:row+1] ≈ [-1/2, 0, 1/2]
    end
    @test Amat[end,end-3:end] ≈ [-1/3, 3/2, -3, 11/6]

end

@testset "None BC FiniteDifference Operator" begin
    N = 10
    d_order = 1
    approx_order = 2
    A = FiniteDifference{Float64}(d_order,approx_order,ones(N-1),N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:3] ≈ [-1.5, 2.0, -0.5]
    for row = 2:N-1
        @test Amat[row,row-1:row+1] ≈ [-0.5, 0.0, 0.5]
    end
    @test Amat[end,end-2:end] ≈ [0.5, -2.0, 1.5]

    d_order = 2
    A = FiniteDifference{Float64}(d_order,approx_order,ones(N-1),N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:4] ≈ [2.0, -5.0, 4.0, -1.0]
    for row = 2:N-1
        @test Amat[row,row-1:row+1] ≈ [1.0, -2.0, 1.0]
    end
    @test Amat[end,end-3:end] ≈ [-1.0, 4.0, -5.0, 2.0]

    N = 10
    d_order = 1
    approx_order = 3
    A = FiniteDifference{Float64}(d_order,approx_order,ones(N-1),N,:None,:None)
    Amat = convert(Array,A)

    @test Amat[1,1:4] ≈ [-11/6, 3, -3/2, 1/3]
    for row = 3:N-2
        @test Amat[row,row-1:row+1] ≈ [-1/2, 0, 1/2]
    end
    @test Amat[end,end-3:end] ≈ [-1/3, 3/2, -3, 11/6]

end

@testset "FiniteDifference Operator sin(x)==d/dx^4 sin(x)" begin
    N = 10
    n = 1000
    x = erf.(range(-2,stop=2,length=n))*2π
    x .+=abs(x[1])
    dx = diff(x)

    A = DiffEqOperators.FiniteDifference{Float64}(4,4,dx,length(x),:None,:None)
    boundary_points = A.boundary_point_count

    y = @. sin(x)+2 #avoid 0.0
    res = A*y
    @test res[boundary_points[1] .+ 1: N - boundary_points[2]] .+ 2.0 ≈ y[boundary_points[1] .+ 1: N - boundary_points[2]] atol=10.0^-1;

end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

__clz_array=[0x0004,0x0003,0x0002,0x0002,
             0x0001,0x0001,0x0001,0x0001,
             0x0000,0x0000,0x0000,0x0000,
             0x0000,0x0000,0x0000,0x0000]

function clz_nobound(n::UInt64)
  (n == 0) && return 64
  res::UInt16 = 0
  #use the binary search method
  (n & 0xFFFF_FFFF_0000_0000 == 0) && (n <<= 32; res += 0x0020)
  (n & 0xFFFF_0000_0000_0000 == 0) && (n <<= 16; res += 0x0010)
  (n & 0xFF00_0000_0000_0000 == 0) && (n <<= 8;  res += 0x0008)
  (n & 0xF000_0000_0000_0000 == 0) && (n <<= 4;  res += 0x0004)
  res + __clz_array[(n >> 60) + 1]
end

function clz_inbound(n::UInt64)
  (n == 0) && return 64
  res::UInt16 = 0
  #use the binary search method
  (n & 0xFFFF_FFFF_0000_0000 == 0) && (n <<= 32; res += 0x0020)
  (n & 0xFFFF_0000_0000_0000 == 0) && (n <<= 16; res += 0x0010)
  (n & 0xFF00_0000_0000_0000 == 0) && (n <<= 8;  res += 0x0008)
  (n & 0xF000_0000_0000_0000 == 0) && (n <<= 4;  res += 0x0004)
  @inbounds res += __clz_array[(n >> 60) + 1]
  res
end

count = 10000

alist = [rand(UInt64) for idx = 1:count]

function testwith(f, a)
  for idx=1:count
    f(a[idx])
  end
end

@time testwith(clz_nobound, alist)
@time testwith(clz_inbound, alist)
@time testwith(clz_nobound, alist)
@time testwith(clz_inbound, alist)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using NCDatasets, DataStructures
ds = Dataset("diffusivity_2_random_phase.nc", "c")
initial_data = θ_A
final_data = ld
ensemble_data = θ_F

defDim(ds, "x", size(initial_data)[1])
defDim(ds, "y", size(initial_data)[2])

v = defVar(ds, "t0", Float32, ("x", "y"), attrib=OrderedDict(
    "units" => " "))
v.attrib["comments"] = "The initial tracer concentration"
v[:, :] = initial_data

w = defVar(ds, "kappa_2e-13", Float32, ("x", "y"), attrib=OrderedDict(
    "units" => " "))
w.attrib["comments"] = "Tracer concentration at t=20"
w[:, :] = final_data

we = defVar(ds, "kappa_2e-8", Float32, ("x", "y"), attrib=OrderedDict(
    "units" => " "))
we.attrib["comments"] = "Tracer concentration at t=20"
we[:, :] = ensemble_data

close(ds)

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
#builds mcmc xsections, currently only used in testing
function xsectionmcmc(Θ::VolumePartsMCMC,::Type{TM}, ::Type{TV}, ::Type{T}) where {TM,TV,T}
  @unpack Xv,ts = Θ

  #=D = (ws=Xv.ws', RLws=Xv.RLws', Ã=Θ.expand[:Ã])
  xsection = dataxsection(D,Θ.ts,AbstractMatrix{T}, AbstractVector{T}, T)
  LD = G=Θ.expand[:G],σ²=Θ.expand[:σ²], ω=Θ.expand[:ω]=#

  dims = (T=size(Θ.Ã,2), K=size(Θ.Ã,1))
  xsections = Vector{NamedTuple}(undef, dims.T-1)
  for t ∈ 2:dims.T
    #count the number of stocks in each window (t and t-1)
    Nt = sum(ts.==t)

    #repeat t for each stock
    #Ainds = repeat([t],Nt)
    #LAinds = repeat([t-1], Nt)
    xinds = ts .== t
    xsections[t-1] = (
      ws = view(Θ.Xv.ws, xinds, :),
      RLws = view(Θ.Xv.RLws, xinds, :),
      v = view(Θ.Xv.v, xinds),
      #X̃ = view(Θ.X̃, xinds, :),
      #ṽ  = view(Θ.ṽ, xinds),
      Ã = view(Θ.expand[:Ã], :, xinds),
      LÃ = view(Θ.expand[:LÃ], :, xinds),
      G = view(Θ.expand[:G], :, xinds),
      σ²   = view(Θ.expand[:σ²], xinds),
      ω = view(Θ.expand[:ω], xinds),
      )
  end
  xsection::NamedTuple = (
    ws = (x->x.ws).(xsections),
    RLws = (x->x.RLws).(xsections),
    v = (x->x.v).(xsections),
    #X̃ = (x->x.X̃).(xsections),
    #ṽ = (x->x.ṽ).(xsections),

    Ã = (x->x.Ã).(xsections),
    LÃ = (x->x.LÃ).(xsections),
    G = (x->x.G).(xsections),
    σ²   = (x->x.σ²).(xsections),
    ω = (x->x.ω).(xsections),
  )

  return xsection
end

function testcapacitychains()
  hyper=(A = (a=1.1,), V=(v1=2.1, v2=2.2), M=(m1=rand(2,2), m2=3.2))
  testparams = (A = 2.3, V=rand(5), M=rand(2,5))
  cc = CapacityChains(Float64, Nburn=25, Ndraws=10, paramstemplate=testparams, hyper=hyper)
  Nparams = 16
  @assert size(cc) === (10,16)

  ###test the indices
  #prior to looping
  for i ∈ 1:10
    @assert i == cc.idx
    @assert i == cc.i
    increment!(cc)
  end

  #after the first loop of the buffer
  for i ∈ 1:15
    @assert (i-1)%10 + 1 == cc.idx
    @assert i + 10 == cc.i
    increment!(cc)
  end

  #after the first reset
  @assert cc.idx == 1
  @assert cc.i == 26
  @assert cc.previdx ==10

  ###test the state vector
  Ψ::CCState = CCState(cc)
  v = rand(Nparams)
  Ψ .= v
  @assert all(v .== Ψ)
  updateandincrement!(cc, Ψ)
  @assert all(cc.value[cc.previdx,:,1] .== Ψ)
  @assert cc.i == 27
  @assert cc.idx == 2
  @assert cc.previdx == 1
  for (n,p) ∈ pairs(Ψ.params) #test each params
    @assert all(reshape(cc.params[n][cc.previdx,:,1],size(p)) .== p)

    #=(all(reshape(cc.params[n][cc.previdx, :, 1],size(p)) .== p)) || error(
      "mismatch!!! cc.idx: $(cc.idx), previdx: $(cc.previdx)
      reshape(cc.params[$n][$(cc.previdx), :, 1],size(p)): $(reshape(cc.params[n][:,cc.previdx,1],size(p)))
      p: $(p)")=#
  end

  #test the params individually
  Ψ.A[] = testparams.A
  Ψ.V .= testparams.V
  Ψ.M .= testparams.M
  @assert all([(v->[v...;]).(values(testparams))...;] .== Ψ)
  updateandincrement!(cc,Ψ)
  @assert (cc.i, cc.idx, cc.previdx) === (28,3,2)
  for (n,p) ∈ pairs(Ψ.params) #test each params
    if n ≡ :A
      @assert reshape(cc.params[:A][cc.previdx,:,1],size(p))[] === testparams.A
    elseif n ≡ :V
      @assert all(reshape(cc.params[:V][cc.previdx,:,1],size(p)) .=== testparams.V)
    elseif n ≡ :M
      @assert all(reshape(cc.params[:M][cc.previdx,:,1],size(p)) .=== testparams.M)
    else
      @assert false
    end
  end

end

function testA₀(Θ::AbstractVolumePartsMCMC{TM,TV,T},
  Ntest::Int=PARAM[:mcmcntest]) where {TM,TV,T}

  @unpack Ψ, Xv, xsection, Ã = deepcopy(Θ)
  @unpack A₀, G, ω, σ² = Ψ
  #@info "hyper keys=$(keys(Θ.cc.hyper))"
  @unpack μA₀, VA₀inv = Θ.cc.hyper[:A₀]

  dims = (K=size(Ã, 1), T=size(Ã, 2), )

  A₀ver = Matrix(undef, Ntest, dims.K)
  A₀test = similar(A₀ver)

  function verA₀!()
    Ã .= hcat(ones(T, dims.K), cumprod(G, dims=2))

    #store the cross-sectional matrices
    #Λₜ = Vector{TM}(undef, dims.T-1)
    #Āₜ = Vector{TV}(undef, dims.T-1)
    VA₀invmat = Matrix(VA₀inv, dims.K, dims.K)
    μA₀vec = fill(μA₀, dims.K)
    Λₜ = [Matrix{T}(undef, dims.K, dims.K) for t ∈ 1:(dims.T-1)]
    ΛĀₜ = [Vector{T}(undef, dims.K) for t ∈ 1:(dims.T-1)]
    τₜ = 1 ./ σ²
    Threads.@threads for t ∈ 1:(dims.T-1)
      X = abs.(xsection[:ws,t] .* xsection[:Ã,t]' .- xsection[:RLws,t] .* xsection[:LÃ,t]')
      XtX = X' * X
      ṽ = xsection[:v,t] .- xsection[:ω,t]
      Â = cholesky(XtX)\(X'*ṽ)
      ΛĀₜ[t] .= τₜ[t] .* (XtX*Â .+ VA₀inv*μA₀vec)
      Λₜ[t] .= τₜ[t] .* (XtX .+ VA₀invmat)

      #=if !(Λₜ[t] ≈ Λₜ[t]') #fixes some roundoff error
        println("XtX: ")
        printmln(XtX)
        println("Λₜ[t]: ")
        printmln(Λₜ[t])
        error("$t: something is wrong with Λₜ[t]: !(Λₜ[t] ≈ Λₜ[t]')")
      end=#
    end

    ΛĀ = sum(ΛĀₜ)
    Λ = sum(Λₜ)

    Σ = Λ\I

    if Σ ≈ Σ' #fixes some roundoff error
      Σ .= (Σ .+ Σ') ./ 2
    else
      printmln(Σ)
      error("something is wrong with Σ: !(Σ ≈ Σ')")
    end

    Ā = Σ*ΛĀ
    #@info sum(abs.(Σ .- Σ'))
    d = MvNormal(Ā, Σ)
    return rand(d)

    #=for k ∈ shuffle(1:dims.K)
      notk = k .≠ 1:dims.K
      Σkkinv = cholesky(Σ[notk,notk])\I

      μk = Ā[k] + Σ[k:k, notk] * Σkkinv * (A₀[notk] .- Ā[notk])=#

  end

  #verA₀!()
  #@assert false

  if PARAM[:mcmcbenchtest]
    print("time updateA₀!: ")
    @btime updateA₀!($Θ)
    print("time verA₀!: ")
    @btime $verA₀!()
  end

  for i ∈ 1:Ntest
    A₀test[i,:] .= updateA₀!(Θ)
    A₀ver[i,:] .= verA₀!()
  end

  #halfNtest = Ntest ÷ 2
  #A₀vermean = mean(A₀ver[halfNtest:end, :], dims=1)
  #A₀testmean = mean(A₀test[halfNtest:end, :], dims=1)

  A₀vermean = mean(A₀ver, dims=1)
  A₀testmean = mean(A₀test, dims=1)

  Δ = abs.(A₀vermean .- A₀testmean)
  @info "
    A₀vermean: $A₀vermean
    A₀testmean: $A₀testmean
    Δ: $Δ"

end

function testposteriors(Θ::AbstractVolumePartsMCMC)
  testA₀(Θ)
end

#we use a lot of views in VolumeParts- this tests them
function testmcmcexpand(Θ::AbstractVolumePartsMCMC{TM,TV,T}, xsectionver) where {TM, TV, T}
  dims = (K=size(Θ.Ã,1), T=size(Θ.Ã,2))
  @unpack ts=Θ

  #these sort the views as though they were first divided into xsections
  function sorttvert(M::AbstractMatrix, ts)
    allts = sort!(unique(ts))
    reduce(vcat, (t->view(M, ts .== t, :)).(allts))
  end
  function sorttvert(M::AbstractVector, ts)
    allts = sort!(unique(ts))
    reduce(vcat, (t->view(M, ts .== t)).(allts))
  end

  function sortthoriz(M::AbstractMatrix, ts)
    tst = ts'
    allts = sort!(unique(ts))
    reduce(hcat, (t->view(M, :, vec(tst .== t))).(allts))
  end

  #concatenate the cross-sections in xsection to verify the Θ.expand views
  vdata = [:ws, :RLws, :v, :σ², :ω, ]
  hdata = [:Ã, :LÃ, :G, ]
  for s ∈ [vdata; hdata]
    local catfunc
    local sortfunc
    if s ∈ vdata #check if we arranged the data vertically or horizontally
      catfunc = vcat
      sortfunc = sorttvert
    elseif s ∈ hdata
      catfunc = hcat
      sortfunc = sortthoriz
    else
      @assert false
    end

    #check the complete expanded view
    vermatrix = reduce(catfunc, xsectionver[s])
    testmatrix = sortfunc(Θ.expand[s], ts)
    testxmatrix = reduce(catfunc, (t->Θ.xsection[s,t]).(1:(dims.T-1)))
    #@assert size(vermatrix) === size(testmatrix) === size(testxmatrix)
    all(vermatrix .== testmatrix .== testxmatrix)||error(
      "!all(vermatrix .== testmatrix .== testxmatrix)
      vermatrix[1:3,1:3]: $(vermatrix[1:3,1:3])
      testmatrix[1:3,1:3]: $(testmatrix[1:3,1:3])
      testxmatrix[1:3,1:3]: $(testxmatrix[1:3,1:3])
      "
    )

  #concatenate the cross-sections in xsection to verify the Θ.expand views (now for each t)
    for t ∈ 1:(dims.T-1)

      partts = @view ts[ts .≥ t + 1]
      vermatrix = reduce(catfunc, xsectionver[s][t:(dims.T-1)])
      testmatrix = sortfunc(Θ.expand[s,t], partts)
      testxmatrix = reduce(catfunc, (t->Θ.xsection[s,t]).(t:(dims.T-1)))
      @assert size(vermatrix) === size(testmatrix) === size(testxmatrix)
      @assert all(vermatrix .== testmatrix .== testxmatrix)

      @assert all(Θ.xsection[s,t] .== xsectionver[s][t])
    end
  end
return nothing
end

function testmcmc(panel::AbstractDataFrame,  zs::ZSpec,
    ::Type{T}=PARAM[:mcmctype],
    ::Type{TV}=PARAM[:mcmcgpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:mcmcgpu] ? CUDA.CuMatrix{T} : Matrix{T};
    prior::Symbol=PARAM[:mcmcprior],
    Nburn::Int=PARAM[:mcmcnburn],
    Ndraws::Int=PARAM[:mcmcndraws]
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  testcapacitychains()

  @unpack Zms = zs
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, zs, T, TV, TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)

  dims = p.dims
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  v::TV = p.dat[:Fvol]
  ts::Vector{Int} = p.ts
  Ascale::Vector{T} = p.Ascale
  Nexpanded::Int = length(ts)
  #this creates the factored objects- which are accurate given that initially G=1
  absΔ(wsᵢ::CuArray,RLwsᵢ)=CUDA.abs(wsᵢ-RLwsᵢ)
  absΔ(wsᵢ,RLwsᵢ)=abs(wsᵢ-RLwsᵢ)

  #now test the level versions
  Xv = XYIter(ws ,RLws, v, ts, Ascale)
  hyper=PARAM[prior]
  initial=(
    A₀ = TV(rand(T, dims.K)),
    G = TM(ones(T, dims.K, dims.T-1)),
    σ² = rand(T, dims.T-1),
    ω = rand(T, dims.T-1)
    )
  cc = CapacityChains(T, Nburn=Nburn, Ndraws=Ndraws, paramstemplate=initial, hyper=hyper)

  Θ = VolumePartsMCMC(TM, TV, T; initial, cc, Xv, ts, )

  xsection = xsectionmcmc(Θ, TM, TV, T)

  testmcmcexpand(Θ,xsection)

  testposteriors(Θ)

#=  testG(Θ, Xvlevel, _Xv, xsectionlevel)
  @assert Θ.G |> Matrix! ≈ Θ.A[:,2:end] ./ Θ.A[:,1:(end-1)] |> Matrix!

  #WARNING WARNING WARNING uncomment below!!!
  testprediction(Θ, Xvlevel, deepcopy(_Xv))
  testpredictionzygote(deepcopy(Θ),Xvlevel,deepcopy(_Xv))
  testiterloss(deepcopy(Θ),Xvlevel,deepcopy(_Xv))

  if !(TM<:CuArray)

    Xvdemean::XYIterDemean = XYIterDemean(ws ,RLws, v, ts, Ascale)
    Ddemean = (ws=Xvdemean.ws, RLws=Xvdemean.RLws, ṽ = Xvdemean.ṽ,)
    xsectiondemean = dataxsection(Ddemean,Θ.ts,TM, TV, T)
    validatexsection(Θ, Xvdemean, _Xv, p, xsectiondemean)

    testG(Θ, Xvdemean, _Xv, xsectiondemean)
    @assert Θ.G |> Matrix! ≈ Θ.A[:,2:end] ./ Θ.A[:,1:(end-1)] |> Matrix!

    testXv(Θ, Xvlevel, Xvdemean, _Xv)

    testprediction(Θ, Xvdemean, deepcopy(_Xv))
    testpredictionzygote(deepcopy(Θ),Xvdemean,deepcopy(_Xv))
    testiterloss(deepcopy(Θ),Xvdemean,deepcopy(_Xv))
  else
    #create cpu versions
    pcpu = modelparts(panel, zs,T, Vector{T}, Matrix{T},
      weightdata = (Fw=Fw, FRLw=FRLw),
      Fvol=Zms.Fvol)

    wscpu::Matrix{T} = pcpu.dat[:Fw]
    RLwscpu::Matrix{T} = pcpu.dat[:FRLw]
    vcpu::Vector{T} = pcpu.dat[:Fvol]

    Xvlevelcpu::XYIterLevel = XYIterLevel(wscpu ,RLwscpu, vcpu, ts, Ascale)
    Θcpu = VolumePartsIter(dims.T, dims.K, ts, T, Matrix{T}, Vector{T})
    Θcpu.G .= Θ.G |> Matrix
    Θcpu.A .= Θ.A |> Matrix
    _Xvcpu = XYIterAlloc(Θ, Xvlevelcpu)

    testcudazygote(deepcopy(Θ),Xvlevel,deepcopy(_Xv), Θcpu, Xvlevelcpu,  _Xvcpu)
  end=#

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


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)
end

@testset "interpolate methods" begin
    x = [12,2,2,1,2,4,3,1]
    y = [12, 7, 2, 2, 2, 1.5, 1, 1.5, 2, 3, 4, 3.5, 3, 2, 1]
    @test y ≈ my_interpolate(x, (15))

    x = collect(0:0.05:3)
    y = my_interpolate(x, (2 * size(x)[1]))
    @test mean(exp.(x)) ≈ mean(exp.(x))

    x = sin.(collect(0:0.001:3))
    y = my_interpolate(x, size(x)[1] * 3 + 2)
    y2 = my_interpolate(y, size(x)[1])
    @test isapprox(y2, x, rtol=1e-6)

    x = [12,2,2,1,2,4,3,1]
    y = [12, 7, 2, 2, 2, 1.5, 1, 1.5, 2, 3, 4, 3.5, 3, 2, 1]
    y_2d = reshape(y, 1, :)
    x_2d = my_interpolate(reshape(x, 1, :), size(y_2d))
    @test isapprox(x_2d, y_2d, rtol=1e-6)

end

@testset "Generate downsample" begin
    ds = generate_downsample(2, [1,2], 2)
    @test [2.5] ≈ ds([1 2; 3 4])

    ds = generate_downsample(2, [2], 2)
    @test [1.5; 3.5; 5.5; 7.5] ≈ ds([1 2; 3 4; 5 6; 7 8])

    ds = generate_downsample(2, [1], 2)
    @test [2.0 3.0; 6.0 7.0] ≈ ds([1 2; 3 4; 5 6; 7 8])
end

@testset "Generate PSF method" begin
    # large aperture is delta peak
    out = zeros((5, 5))
    out[1,1] = 1

    @test out ≈ generate_psf((5, 5), 100)

    # pinhole aperture
    out = ones((10, 10))
    out ./= sum(out)
    @test out ≈ generate_psf((10, 10), 0.01)

    # normalized
    @test 1 ≈ sum(generate_psf((100, 100), 10))
end

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

    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 ≈ DeconvOptim.rr_3D((2,2,2))

end

@testset "next fast fft size" begin
    @test DeconvOptim.next_fast_fft_size(23) == 24
    @test DeconvOptim.next_fast_fft_size((23, 46)) == (24, 48)
end


text:
# ---
# title: 61. Rotate List
# id: problem61
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Two Pointers
# link: <https://leetcode.com/problems/rotate-list/description/>
# hidden: true
# ---
#
# Given the `head` of a linked list, rotate the list to the right by `k` places.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/rotate1.jpg)
#
#
#
#     Input: head = [1,2,3,4,5], k = 2
#     Output: [4,5,1,2,3]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/roate2.jpg)
#
#
#
#     Input: head = [0,1,2], k = 4
#     Output: [2,0,1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is in the range `[0, 500]`.
#   * `-100 <= Node.val <= 100`
#   * `0 <= k <= 2 * 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@info "Installing packages. This may take up to a few minutes."

using Pkg
Pkg.activate(".")
Pkg.instantiate()

@info "Loading packages"

# Standard packages
using Random
using Printf
using Test
using DelimitedFiles

using JuMP
const MOI = JuMP.MathOptInterface

using GLPK
using Cbc

@info "All packages successfully imported."


text:
using Test

using MarkovChains
@testset "solve" begin
    @testset "irreducible1" begin
        chain = ContMarkovChain()
        n0 = add_state!(chain)
        n1 = add_state!(chain)
        n2 = add_state!(chain)
        n3 = add_state!(chain)
        add_transition!(chain, n0, n1, 1.0)
        add_transition!(chain, n1, n2, 1.0)
        add_transition!(chain, n2, n3, 1.0)
        add_transition!(chain, n3, n2, 3.0)
        add_transition!(chain, n2, n1, 2.0)
        add_transition!(chain, n1, n0, 1.0)
        init_prob = [0.1, 0.9, 0.0, 0.0]
        ts_sol = solve(chain, init_prob, 0.0)
        @test ts_sol.prob ≈ [0.1, 0.9, 0.0, 0.0]
        @test ts_sol.cumtime ≈ [0.0, 0.0, 0.0, 0.0]

        ss_sol = solve(chain, init_prob, Inf)
        @test ss_sol.prob ≈ [0.375, 0.375, 0.1875, 0.0625]
        @test ss_sol.cumtime == [Inf, Inf, Inf, Inf]
    end
    @testset "acyclic1" begin
    chain = ContMarkovChain()
    n1 = add_state!(chain)
    n2 = add_state!(chain)
    n3 = add_state!(chain)
    add_transition!(chain, n1, n2, 2.0)
    add_transition!(chain, n2, n3, 4.0)
    init_prob = [1.0, 0.0, 0.0]
    ts_sol = solve(chain, init_prob, 0.0)
    @test ts_sol.prob ≈ [1.0, 0.0, 0.0]

    ts_sol = solve(chain, init_prob, 1.0)
    @test state_prob(ts_sol, 1) < 1.0
    @test state_prob(ts_sol, 2) > 0.0
    @test state_prob(ts_sol, 3) > 0.0

    ss_sol = solve(chain, init_prob, Inf)
    @test state_cumtime(ss_sol, 1) == 0.5
    @test state_cumtime(ss_sol, 2) == 0.25
    @test state_cumtime(ss_sol, 3) == Inf
    @test ss_sol.prob ≈ [0.0, 0.0, 1.0]

    end
end


text:
module Messaging

using ..Pages
import Pages: connections

export broadcast, message, block, add_library

"""
Broadcast a message to all connected web pages to be interpreted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function broadcast(route,args::Dict)
    if haskey(connections,route)
        for (cid,client) in connections[route]
            if isopen(client)
                write(client, JSON.json(args))
            end
        end
    end
end
broadcast(r,t,d) = broadcast(r,Dict("type"=>t,"data"=>d))

function broadcast(args::Dict)
    for route in keys(connections)
        for (cid,client) in connections[route]
            if isopen(client)
                write(client, JSON.json(args))
            end
        end
    end
end
broadcast(t,d) = broadcast(Dict("type"=>t,"data"=>d))
broadcast(d) = broadcast(Dict("type"=>"say","data"=>d))

"""
Send a message to the specified connection to be interpreted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function message(route,args::Dict)
    id = pop!(args,"id");
    if haskey(connections,route) && haskey(connections[route],id)
        client = connections[route][id]
        if isopen(client)
            write(client,JSON.json(args))
        end
    end
end
message(route,mid,mtype,mdata) = message(route,Dict("id" => mid, "type" => mtype, "data" => mdata))
message(route,mid,mdata) = message(route,Dict("id" => mid, "type" => "say", "data" => mdata))

"""
Block Julia control flow until until callback["notify"](name) is called.
"""
function block(f::Function,name)
    conditions[name] = Condition()
    f()
    wait(conditions[name])
    delete!(conditions,name)
    return nothing
end

end


text:
######################### Oblivious Algorithm ##################################
# chooses the k features most correlated with the target
function oblivious(A::AbstractMatOrFac, b::AbstractVector, k::Int)
    nzind = partialsortperm(abs.(A'b), 1:k, rev = true)
    x = spzeros(size(b))
    x[nzind] = @view(A[:, nzind]) \ b
    return x
end


text:
module Tentacles

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Tentacles" => Ahorn.EntityPlacement(
        Maple.Tentacles,
        "point",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + 32, Int(entity.data["y"]))]
        end
    )
)

# Maple definition doesn't have "Fake keys"
const fearDistance = Dict{String, String}(
    "None" => "",
    "Close" => "close",
    "Medium" => "medium",
    "Far" => "far"
)

Ahorn.nodeLimits(entity::Maple.Tentacles) = 1, -1

Ahorn.editingOptions(entity::Maple.Tentacles) = Dict{String, Any}(
    "fear_distance" => fearDistance
)

function Ahorn.selection(entity::Maple.Tentacles)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.Rectangle(x - 12, y - 12, 24, 24)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.Rectangle(nx - 12, ny - 12, 24, 24))
    end

    return res
end

function drawTentacleIcon(ctx::Ahorn.Cairo.CairoContext, x, y)
    Ahorn.drawImage(ctx, Ahorn.Assets.tentacle, x - 12, y - 12)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Tentacles)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        drawTentacleIcon(ctx, nx, ny)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Tentacles, room::Maple.Room)
    x, y = Ahorn.position(entity)
    drawTentacleIcon(ctx, x, y)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn