text:
using LinearAlgebra;

function blocking(A, delta = 0.1)
    #BLOCKING  Produce blocking pattern for block Parlett recurrence in FUNM.
    #   M = BLOCKING(A, DELTA, SHOWPLOT) accepts an upper triangular matrix
    #   A and produces a blocking pattern, specified by the vector M,
    #   for the block Parlett recurrence.
    #   M(i) is the index of the block into which A(i,i) should be placed,
    #   for i=1:LENGTH(A).
    #   DELTA is a gap parameter (default 0.1) used to determine the blocking.
    #
    #   For A coming from a real matrix it should be posible to take
    #   advantage of the symmetry about the real axis.  This code does not.

    a = diag(A); n = length(a);
    m = zeros(LinearAlgebra.BlasInt, n); maxM = 0;

    for i = 1 : n

        if m[i] == 0
            m[i] = maxM + 1; # If a(i) hasn`t been assigned to a set
            maxM = maxM + 1; # then make a new set and assign a(i) to it.
        end

        for j = i+1 : n
            if m[i] != m[j]    # If a(i) and a(j) are not in same set.
                if abs(a[i]-a[j]) <= delta

                    if m[j] == 0
                        m[j] = m[i]; # If a(j) hasn`t been assigned to a
                        # set, assign it to the same set as a(i).
                    end

                    for j = i+1 : n
                        if m[i] != m[j]    # If a(i) and a(j) are not in same set.
                            if abs(a[i]-a[j]) <= delta

                                if m[j] == 0
                                    m[j] = m[i]; # If a(j) hasn`t been assigned to a
                                    # set, assign it to the same set as a(i).
                                else
                                    p = max(m[i],m[j]);
                                    q = min(m[i],m[j]);

                                    m[m .== p] .= q; # If a(j) has been assigned to a set
                                    # place all the elements in the set
                                    # containing a(j) into the set
                                    # containing a(i) (or vice versa).
                                    m[m .> p] .= m[m .> p] .- 1;
                                    maxM = maxM - 1;
                                    # Tidying up. As we have deleted set
                                    # p we reduce the index of the sets
                                    # > p by 1.
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    return m
end

function largest_block(ind)
    n = length(ind);
    v = zeros(LinearAlgebra.BlasInt, n);

    for i = 1 : n
        v[ind[i]] += 1;
    end

    return maximum(v)
end

function merge_blocks(ord, ind, block_size)
    #MERGE_BLOCKS returns [ord, ind]

    ind2 = [];
    cluster = [];

    for j = 1 : length(ind)
        if length(cluster) < block_size
            cluster = [ cluster ind{j} ];
        else
            ind2 = [ ind2 cluster ];
            ord[cluster] = j;
            cluster = [ ind[j] ];
        end
    end

    if cluster != []
        ind2 = [ ind2 cluster ];
        ord[cluster] = j;
    end

    ind = ind2;

end

function swapping(m)
    #SWAPPING  Choose confluent permutation ordered by average index.
    #   [MM,IND] = SWAPPING(M) takes a vector M containing the integers
    #   1:K (some repeated if K < LENGTH(M)), where M(J) is the index of
    #   the block into which the element T(J,J) of a Schur form T
    #   should be placed.
    #   It constructs a vector MM (a permutation of M) such that T(J,J)
    #   will be located in the MM(J)'th block counting from the (1,1) position.
    #   The algorithm used is to order the blocks by ascending
    #   average index in M, which is a heuristic for minimizing the number
    #   of swaps required to achieve this confluent permutation.
    #   The cell array vector IND defines the resulting block form:
    #   IND{i} contains the indices of the i'th block in the permuted form.

    mmax = maximum(m); mm = zeros(LinearAlgebra.BlasInt, size(m));
    g = zeros(mmax); h = zeros(mmax);

    for i = 1:mmax
        p = findall(m -> m==i, m);
        h[i] = length(p);
        g[i] = sum(p)/h[i];
    end

    y = sortperm(g);
    h = [0 ; cumsum(h[y])];

    ind = [];
    for i = 1 : mmax
        mm[m .== y[i]] .= i;
        rg = convert(Array{LinearAlgebra.BlasInt}, h[i]+1 : h[i+1])
        ind = [ ind ; [ rg ] ];
    end

    return mm, ind;
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon the Sudoku module being defined first.

using Test

function puzzle_2()
    puzzle_2_reference = Array{Int16}(undef,(4,4))
    puzzle_2_reference[:] = [1,2,3,4, 3,4,1,2, 2,3,4,1, 4,1,2,3][:]
    @test Sudoku.puzzle_n(2) == puzzle_2_reference
end

function puzzle_rank_array()
    @test Sudoku.get_rank(Sudoku.puzzle_n(2)) == 2
    @test Sudoku.get_rank(Sudoku.puzzle_n(3)) == 3
    @test Sudoku.get_rank(Sudoku.puzzle_n(4)) == 4
end

function puzzle_3_entry()
    rank = 3
    rank_squared = rank*rank
    test_entry = Sudoku.make_entry(rank_squared,0)
    @test Sudoku.get_value(test_entry) == 0
    for i = 1:rank_squared
        test_entry = Sudoku.make_entry(rank_squared,i)
        @test Sudoku.get_value(test_entry) == i
    end
end

function solvable_puzzle_2_construction()
    puzzle_2_reference = Array{Int16}(undef,(4,4))
    puzzle_2_reference[:] = [1,2,3,4, 3,4,1,2, 2,3,4,1, 4,1,2,3][:]
    p2 = Sudoku.SolvablePuzzle(2)
    @test length(p2.grid) == 16
    for i = 1:16
        @test Sudoku.get_value(p2.grid[i]) == puzzle_2_reference[i]
        @test p2.grid[i].possibilities == (0x1 << (puzzle_2_reference[i]-1))
    end

    # Verify text conversion
    puzzle_2_reference_text = Array{String}(undef,(4,4))
    puzzle_2_reference_text[:] = ["1","2","3","4", "3","4","1","2", "2","3","4","1", "4","1","2","3"][:]
    @test Sudoku.as_text(p2) == puzzle_2_reference_text

    # Verify zero uncertainty
    @test Sudoku.uncertainty(p2) == 0
end

function solvable_puzzle_2_modification()
    puzzle_2_reference = Array{Int8}(undef,(4,4))
    puzzle_2_reference[:] = [0,2,3,4, 3,0,1,2, 2,3,0,1, 4,1,2,0][:]
    p2 = Sudoku.SolvablePuzzle(Int8(2))
    # Make the diagonal unknowns
    for i = 1:4
        Sudoku.set_unknown(p2,i,i)
    end

    for i = 1:4
        for j = 1:4
            @test Sudoku.get_value(p2.grid[i,j]) == puzzle_2_reference[i,j]
            if i == j
                @test p2.grid[i,j].possibilities == 0xf
            else
                @test p2.grid[i,j].possibilities == 0x1 << (puzzle_2_reference[i,j]-1)
            end
        end
    end

    # Verify value extraction
    puzzle_2_reference_values = Array{UInt8}(undef,(4,4))
    puzzle_2_reference_values[:] = [0,2,3,4, 3,0,1,2, 2,3,0,1, 4,1,2,0][:]
    @test Sudoku.as_values(p2) == puzzle_2_reference_values

    # Verify possibility extraction
    puzzle_2_reference_possibilities = Array{UInt16}(undef,(4,4))
    puzzle_2_reference_possibilities[:] = [0xf,0x2,0x4,0x8, 0x4,0xf,0x1,0x2, 0x2,0x4,0xf,0x1, 0x8,0x1,0x2,0xf][:]
    @test Sudoku.as_possibilities(p2) == puzzle_2_reference_possibilities

    # Verify text conversion
    puzzle_2_reference_text = Array{String}(undef,(4,4))
    puzzle_2_reference_text[:] = [" ","2","3","4", "3"," ","1","2", "2","3"," ","1", "4","1","2"," "][:]
    @test Sudoku.as_text(p2) == puzzle_2_reference_text

    # Verify 4 totally unknown entries
    @test Sudoku.uncertainty(p2) == 16
end

function solvable_puzzle_2_assignment()
    puzzle_2_reference = Array{Int16}(undef,(4,4))
    puzzle_2_reference[:] = [1,2,0,3, 0,3,1,2, 2,0,3,1, 3,1,2,0][:]
    p2 = Sudoku.SolvablePuzzle(2)
    Sudoku.assign_values!(p2,puzzle_2_reference)
    @test length(p2.grid) == 16
    for i = 1:16
        @test Sudoku.get_value(p2.grid[i]) == puzzle_2_reference[i]
        if puzzle_2_reference[i] == 0
            @test p2.grid[i].possibilities == 0xf
        else
            @test p2.grid[i].possibilities == (0x1 << (puzzle_2_reference[i]-1))
        end
    end
end

ext: jl

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
export bachelierVolatility

const alpha4 = 0.15
const aLFK4n = 7
const aLFK4d = 5

const aLFK4 = [0.06155371425063157, 2.723711658728403, 10.83806891491789, 301.0827907126612, 1082.864564205999, 790.7079667603721, 109.330638190985, 0.1515726686825187, 1.436062756519326, 118.6674859663193, 441.1914221318738, 313.4771127147156, 40.90187645954703]

const cLFK4n = 9
const cLFK4d = 7

const bLFK4 = [0.6409168551974357, 788.5769356915809, 445231.8217873989, 149904950.4316367, 32696572166.83277, 4679633190389.852, 420159669603232.9, 2.053009222143781e+16, 3.434507977627372e+17, 2.012931197707014e+16, 644.3895239520736, 211503.4461395385, 42017301.42101825, 5311468782.258145, 411727826816.0715, 17013504968737.03, 247411313213747.3]
const cLFK4 = [0.6421106629595358, 654.5620600001645, 291531.4455893533, 69009535.38571493, 9248876215.120627, 479057753706.175, 9209341680288.471, 61502442378981.76, 107544991866857.5, 63146430757.94501, 437.9924136164148, 90735.89146171122, 9217405.224889684, 400973228.1961834, 7020390994.356452, 44654661587.93606, 76248508709.85633]
const dLFK4 = [0.936024443848096, 328.5399326371301, 177612.3643595535, 8192571.038267588, 110475347.0617102, 545792367.0681282, 1033254933.287134, 695066365.5403566, 123629089.1036043, 756.3653755877336, 173.9755977685531, 6591.71234898389, 82796.56941455391, 396398.9698566103, 739196.7396982114, 493626.035952601, 87510.31231623856]

const betaStart = -log(alpha4)
const betaEnd = -log(floatmin(Float64))
function eta(u::T) where {T}
	if u < 1e-2
		return 1 - u*(0.5+u*(1.0/12+u*(1.0/24+u*(19.0/720+u*(3.0/160+u*(863.0/60480+u*275.0/24192))))))
	end
	return -u / log1p(-u)
end

function ratval(cofs::Vector{Float64}, nn::Int, dd::Int, x::T)  where {T}
	sumn = cofs[nn]
	sumd = 0.0
	for j = nn - 1:-1:1
		sumn = sumn*x + cofs[j]
	end
	for j = nn + dd - 1:-1:nn
		sumd = (sumd + cofs[j]) * x
	end
	return sumn / (1.0 + sumd)
end

"Implied volatility of a vanilla option under the Bachelier/Normal model"
function bachelierVolatility(price::T, isCall::Bool, strike::T, forward::T, tte::T, df::T) where {T}
	price/=df
	sign = 1.0
	if !isCall
		sign = -1.0
	end
	x = (forward - strike) * sign
	if abs(x) < eps()
		return price * sqrt(2*pi/tte)
	end
	z = 0.0
	if x > 0
		#ITM in compute otm
		z = (price - x) / x
	else
		z = -price / x
	end
	vol = 0.0
	if z < alpha4  #highly OTM options up to 1e-100
		u = -(log(z)+betaStart)/(betaEnd-betaStart)
		hz = 0.0
		if u < 0.0091
			num = bLFK4[1+0] + u*(bLFK4[1+1]+u*(bLFK4[1+2]+u*(bLFK4[1+3]+u*(bLFK4[1+4]+u*(bLFK4[1+5]+u*(bLFK4[1+6]+u*(bLFK4[1+7]+u*(bLFK4[1+8]+u*bLFK4[1+9]))))))))
			den = 1.0 + u*(bLFK4[1+10]+u*(bLFK4[1+11]+u*(bLFK4[1+12]+u*(bLFK4[1+13]+u*(bLFK4[1+14]+u*(bLFK4[1+15]+u*(bLFK4[1+16])))))))
			hz = num / den
		elseif u < 0.088
			num = cLFK4[1+0] + u*(cLFK4[1+1]+u*(cLFK4[1+2]+u*(cLFK4[1+3]+u*(cLFK4[1+4]+u*(cLFK4[1+5]+u*(cLFK4[1+6]+u*(cLFK4[1+7]+u*(cLFK4[1+8]+u*cLFK4[1+9]))))))))
			den = 1.0 + u*(cLFK4[1+10]+u*(cLFK4[1+11]+u*(cLFK4[1+12]+u*(cLFK4[1+13]+u*(cLFK4[1+14]+u*(cLFK4[1+15]+u*(cLFK4[1+16])))))))
			hz = num / den
		 else
			num = dLFK4[1+0] + u*(dLFK4[1+1]+u*(dLFK4[1+2]+u*(dLFK4[1+3]+u*(dLFK4[1+4]+u*(dLFK4[1+5]+u*(dLFK4[1+6]+u*(dLFK4[1+7]+u*(dLFK4[1+8]+u*dLFK4[1+9]))))))))
			den = 1.0 + u*(dLFK4[1+10]+u*(dLFK4[1+11]+u*(dLFK4[1+12]+u*(dLFK4[1+13]+u*(dLFK4[1+14]+u*(dLFK4[1+15]+u*(dLFK4[1+16])))))))
			hz = num / den
		end
		vol = abs(x) / (sqrt(hz * tte))
	else
		putPrice = price
		if x <= 0
			putPrice = price - x
		end
		z = abs(x) / putPrice
		u = eta(z)
		num = aLFK4[1+0] + u*(aLFK4[1+1]+u*(aLFK4[1+2]+u*(aLFK4[1+3]+u*(aLFK4[1+4]+u*(aLFK4[1+5]+u*(aLFK4[1+6]+u*(aLFK4[1+7])))))))
		den = 1.0 + u*(aLFK4[1+8]+u*(aLFK4[1+9]+u*(aLFK4[1+10]+u*(aLFK4[1+11]+u*(aLFK4[1+12])))))
		hz = num / den

		vol = sqrt(1/tte) * putPrice * hz
	end
	return vol
end


text:
@testset "Restriction" begin
    A = reshape([UInt16(i) for i = 1:60], 4, 5, 3)
    B = restrict(A, (1,2))
    Btarget = cat(   [  0.96875   4.625   5.96875;
                        2.875    10.5    12.875;
                        1.90625   5.875   6.90625],
                     [  8.46875  14.625  13.46875;
                       17.875    30.5    27.875;
                        9.40625  15.875  14.40625],
                     [ 15.96875  24.625  20.96875;
                       32.875    50.5    42.875;
                       16.90625  25.875  21.90625], dims=3)
    @test B ≈ Btarget
    Argb = reinterpretc(RGB, reinterpret(N0f16, permutedims(A, (3,1,2))))
    B = restrict(Argb)
    Bf = permutedims(reinterpretc(eltype(eltype(B)), B), (2,3,1))
    # isapprox no longer lies, so atol is now serious
    @test isapprox(Bf, Btarget/reinterpret(one(N0f16)), atol=1e-10)
    Argba = reinterpretc(RGBA{N0f16}, reinterpret(N0f16, A))
    B = restrict(Argba)
    @test isapprox(reinterpretc(eltype(eltype(B)), B), restrict(A, (2,3))/reinterpret(one(N0f16)), atol=1e-10)
    A = reshape(1:60, 5, 4, 3)
    B = restrict(A, (1,2,3))
    @test cat(   [  2.6015625   8.71875   6.1171875;
                    4.09375    12.875     8.78125;
                    3.5390625  10.59375   7.0546875],
                 [ 10.1015625  23.71875  13.6171875;
                   14.09375    32.875    18.78125;
                   11.0390625  25.59375  14.5546875], dims=3) ≈ B
    # Issue #395
    img1 = colorview(RGB, fill(0.9, 3, 5, 5))
    img2 = colorview(RGB, fill(N0f8(0.9), 3, 5, 5))
    @test isapprox(channelview(restrict(img1)), channelview(restrict(img2)), rtol=0.01)
    # Non-1 indices
    Ao = OffsetArray(A, (-2,1,0))
    @test parent(@inferred(restrict(Ao, 1))) == restrict(A, 1)
    @test parent(@inferred(restrict(Ao, 2))) == restrict(A, 2)
    @test parent(@inferred(restrict(Ao, (1,2)))) == restrict(A, (1,2))
    # Arrays-of-arrays
    a = Vector{Int}[[3,3,3], [2,1,7],[-11,4,2]]
    @test restrict(a) == Vector{Float64}[[2,3.5/2,6.5/2], [-5,4.5/2,5.5/2]]
    # Images issue #652
    img = testimage("cameraman")
    @test eltype(@inferred(restrict(img))) == Gray{Float32}
    img = testimage("mandrill")
    @test eltype(@inferred(restrict(img))) == RGB{Float32}
    @test eltype(@inferred(restrict(Lab.(img)))) == RGB{Float32}
    img = rand(RGBA{N0f8}, 11, 11)
    @test eltype(@inferred(restrict(img))) == RGBA{Float32}
    @test eltype(@inferred(restrict(LabA.(img)))) == ARGB{Float32}
end

@testset "Image resize" begin
    testcolor = (RGB,Gray)
    testtype = (Float32, Float64, N0f8, N0f16)

    @testset "Interface" begin
        function test_imresize_interface(img, outsz, args...; kargs...)
            img2 = @test_broken @inferred imresize(img, args...; kargs...) # FIXME: @inferred failed
            img2 = @test_nowarn imresize(img, args...; kargs...)
            @test size(img2) == outsz
            @test eltype(img2) == eltype(img)
        end
        for C in testcolor, T in testtype
            img = rand(C{T},10,10)

            test_imresize_interface(img, (5,5), (5,5))
            test_imresize_interface(img, (5,5), (1:5,1:5)) # FIXME: @inferred failed
            test_imresize_interface(img, (5,5), 5,5)
            test_imresize_interface(img, (5,5), 1:5,1:5) # FIXME: @inferred failed
            test_imresize_interface(img, (5,5), ratio = 0.5)
            test_imresize_interface(img, (20,20), ratio = 2)
            test_imresize_interface(img, (20,20), ratio = (2, 2))
            test_imresize_interface(img, (20,10), ratio = (2, 1))
            test_imresize_interface(img, (10,10), ())
            test_imresize_interface(img, (5,10), 5)
            test_imresize_interface(img, (5,10), (5,))
            test_imresize_interface(img, (5,10), 1:5) # FIXME: @inferred failed
            test_imresize_interface(img, (5,10), (1:5,)) # FIXME: @inferred failed

            @test_throws MethodError imresize(img,5.0,5.0)
            @test_throws MethodError imresize(img,(5.0,5.0))
            @test_throws MethodError imresize(img,(5, 5.0))
            @test_throws MethodError imresize(img,[5,5])
            @test_throws UndefKeywordError imresize(img)
            @test_throws DimensionMismatch imresize(img,(5,5,5))
            @test_throws ArgumentError imresize(img, ratio = -0.5)
            @test_throws ArgumentError imresize(img, ratio = (-0.5, 1))
            @test_throws DimensionMismatch imresize(img, ratio=(5,5,5))
            @test_throws DimensionMismatch imresize(img,(5,5,1))
        end
    end

    @testset "Numerical" begin
        # RGB is skipped
        for C in (Gray,), T in testtype
            A = [1 0; 0 1] .|> C{T}
            @test imresize(A, (1, 1)) ≈ reshape([0.5], 1, 1) .|> C{T}
            # if only 1 of 4 pixels is nonzero, check that you get 1/4 the
            # value (no matter where it is)
            for i = 1:4
                fill!(A, 0)
                A[i] = 1
                A = C{T}.(A)
                @test imresize(A, (1, 1)) ≈ reshape([0.25], 1, 1) .|> C{T}
            end
            A = [1 0; 0 0] .|> C{T}
            @test imresize(A, (2, 1)) ≈ reshape([0.5, 0], (2, 1)) .|> C{T}
            @test imresize(A, (1, 2)) ≈ reshape([0.5, 0], (1, 2)) .|> C{T}
            A = [1 0 0; 0 0 0; 0 0 0] .|> C{T}
            @test imresize(A, (2, 2)) ≈ [0.75^2 0; 0 0] .|> C{T}

            A = ones(5,5) .|> C{T}
            for l2 = 3:7, l1 = 3:7
                R = imresize(A, (l1, l2))
                @test all(x->x==1, R)
            end
            for l2 = 3:7, l1 = 3:7
                R = imresize(A, (0:l1-1, 0:l2-1))
                @test all(x->x==1, R)
                @test axes(R) == (0:l1-1, 0:l2-1)
            end
            for l1 = 3:7
                R = imresize(A, (0:l1-1,))
                @test all(x->x==1, R)
                @test axes(R) == (0:l1-1, 1:5)
            end

            R = imresize(A, (2:6, 0:4))
            @test all(x->x==1, R)
            @test axes(R) == (2:6, 0:4)
            R = imresize(A, ())
            @test all(x->x==1, R)
            @test axes(R) == axes(A)
            @test !(R === A)
            Ao = OffsetArray(A, -2:2, 0:4)
            R = imresize(Ao, (5,5))
            @test axes(R) === axes(A)
            R = imresize(Ao, axes(Ao))
            @test axes(R) === axes(Ao)
            @test !(R === A)
            img = reshape([0.5])
            R = imresize(img, ())
            @test ndims(R) == 0
            @test !(R === A)
        end
    end
end


text:
@with_kw mutable struct BatchAdjacentBuffer{T1, T2}
    batches::Vector{T1} = []
    z_dists::Vector{T2} = []
    elements::Int = 0 # max capacity
    next_ind::Int = 1 # Location of the next index for the circular buffer
    total_count::Int = 0
    default_z_dist
    reference_buffer
end

isprioritized(b::BatchAdjacentBuffer) = false

extra_columns(b::BatchAdjacentBuffer) = extra_columns(b.reference_buffer)

DataStructures.capacity(b::BatchAdjacentBuffer) = b.elements

Base.length(buffer::BatchAdjacentBuffer) = length(buffer.batches) > 0 ? sum([length(b) for b in buffer.batches]) : 0

function device(b::BatchAdjacentBuffer)
    d = device(b.batches[1])
    for e in b.batches
        @assert device(e) == d
    end
    return d
end

function buffer_like(b::BatchAdjacentBuffer; capacity=capacity(b.reference_buffer[1]), device=device(b.reference_buffer))
    buffer_like(b.reference_buffer, capacity=capacity, device=device)
end

Base.push!(b::BatchAdjacentBuffer, data) = nothing #push!(b, data, nothing)

function Base.push!(b::BatchAdjacentBuffer, data, z_dist)
    # Convert data to an experience buffer if needed
    if !(data isa ExperienceBuffer)
        data = ExperienceBuffer(data)
    end

    # Add the new buffer to the array
    if length(b.batches) < capacity(b)
        push!(b.batches, data)
        push!(b.z_dists, z_dist)
    else
        b.batches[b.next_ind] = data
        b.z_dists[b.next_ind] = z_dist
    end

    b.next_ind = mod1(b.next_ind + 1, b.elements)
end

function push_reservoir!(b::BatchAdjacentBuffer, data, z_dist; weight=nothing)
    # Convert data to an experience buffer if needed
    if !(data isa ExperienceBuffer)
        data = ExperienceBuffer(data)
    end

    # Only add samples according to their weight
    if !isnothing(weight) && rand() > weight
        return
    end

    b.total_count += 1

    # Add the element to the buffer and return if we haven't hit our max
    if length(b.batches) < capacity(b)
        push!(b.batches, data)
        push!(b.z_dists, z_dist)
    else
        # Choose a random number up to count
        j = rand(1:b.total_count)

        # If its within the buffer replace the element
        if j <= capacity(b)
            b.batches[j] = data
            b.z_dists[j] = z_dist
        end
    end
end

function Random.rand!(target::ExperienceBuffer, source::BatchAdjacentBuffer...; i = 1, fracs = ones(length(source))./length(source))
    lens = [length.(source) ...]
    if any(lens .== 0)
        fracs[lens .== 0] .= 0
        fracs ./= sum(fracs)
    end
    batches = split_batches(capacity(target), fracs)
    for (b, B) in zip(source, batches)
        B == 0 && continue
        for i=1:B
            batch_sizes = length.(b.batches)
            d = Categorical(batch_sizes ./ sum(batch_sizes))
            i1 = rand(d)
            i2 = rand(1:length(b.batches[i1]))
            push!(target, b.batches[i1], ids = [i2])
        end
    end
end

function log_episode_averages(buffer::BatchAdjacentBuffer, keys, period)
    (;kwargs...) -> begin
        d = Dict()
        indices = get_last_N_indices(buffer, period)
        for k in keys
            val = 0
            ends = 0
            for i in indices
                val += sum(buffer.batches[i][k])
                ends += sum(buffer.batches[i][:episode_end])
            end
            d[Symbol(string("avg_", k))] = val / ends
        end
        d
    end
end

function get_last_N_indices(b::BatchAdjacentBuffer, N)
    N = min(length(b.batches), N)
    C = capacity(b)
    start_index = mod1(b.next_ind - N, C)
    mod1.(start_index:start_index + N - 1, C)
end


text:
# Autogenerated wrapper script for ImageMagick_jll for x86_64-unknown-freebsd
export identify, imagemagick_convert, libwand, mogrify, montage

using Zlib_jll
using libpng_jll
using JpegTurbo_jll
using Libtiff_jll
JLLWrappers.@generate_wrapper_header("ImageMagick")
JLLWrappers.@declare_library_product(libwand, "libMagickWand-6.Q16.so.7")
JLLWrappers.@declare_executable_product(identify)
JLLWrappers.@declare_executable_product(imagemagick_convert)
JLLWrappers.@declare_executable_product(mogrify)
JLLWrappers.@declare_executable_product(montage)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, libpng_jll, JpegTurbo_jll, Libtiff_jll)
    JLLWrappers.@init_library_product(
        libwand,
        "lib/libMagickWand-6.Q16.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        identify,
        "bin/identify",
    )

    JLLWrappers.@init_executable_product(
        imagemagick_convert,
        "bin/convert",
    )

    JLLWrappers.@init_executable_product(
        mogrify,
        "bin/mogrify",
    )

    JLLWrappers.@init_executable_product(
        montage,
        "bin/montage",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export make_polygons
export path2turn
export link_turns!

function make_polygons(paths)
    turn_dict = Dict(map(path2turn, paths))
    collector = []
    while !isempty(turn_dict)
        push!(collector, link_turns!(turn_dict))
    end
    return collector
end

function path2turn(path)
    (path[1], path[2]) => (path[2], path[3])
end

function link_turns!(turn_dict)
    carrier = pop!(turn_dict, first(first(turn_dict)), "empty")
    polygon_atoms = [first(carrier)]
    polygon_bonds = [carrier]
    while (carrier = pop!(turn_dict, carrier, "empty")) != "empty"
        push!(polygon_atoms, first(carrier))
        push!(polygon_bonds, carrier)
    end
    return polygon_atoms, polygon_bonds
end


text:
srand(100)

if haskey(ENV, "MOCHA_USE_CUDA")
  const test_gpu = true
else
  const test_gpu = false
end
const test_cpu   = true

using Mocha
using Base.Test

if test_cpu
  backend_cpu = CPUBackend()
  init(backend_cpu)
end

if test_gpu
  backend_gpu = GPUBackend()
  init(backend_gpu)
end

# run test in the whole directory, latest modified files
# are run first, this makes waiting time shorter when writing
# or modifying unit-tests
function test_dir(dir)
  map(reverse(Mocha.glob(dir, r".*\.jl$", sort_by=:mtime))) do file
    include("$dir/$file")
  end
end

############################################################
# Solvers
############################################################
include("solvers/test-adam-solver.jl")
include("solvers/test-sgd-solver.jl")

############################################################
# Network
############################################################
include("net/topology.jl")
include("net/test-gradient-simple.jl")

############################################################
# Utilities functions
############################################################
include("utils/ref-count.jl")
include("utils/glob.jl")
include("utils/blas.jl")
include("utils/blob-reshape.jl")

if test_gpu
  include("cuda/padded-copy.jl")
  include("cuda/cuvec.jl")
  include("cuda/mocha.jl")
  include("cuda/cublas.jl")
end

############################################################
# Activation Functions
############################################################
include("neurons/relu.jl")
include("neurons/sigmoid.jl")
include("neurons/tanh.jl")
include("neurons/exponential.jl")

############################################################
# Regularizers
############################################################
include("regularizers/l2.jl")
include("regularizers/l1.jl")

############################################################
# Constraints
############################################################
include("constraints/l2.jl")

############################################################
# Data Transformers
############################################################
include("data-transformers.jl")

############################################################
# Layers
############################################################
test_dir("layers")

if test_gpu
  shutdown(backend_gpu)
end
if test_cpu
  shutdown(backend_cpu)
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iup_plot.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const iup_plot = "iup_plot") : (const iup_plot = "libiup_plot")  # Name of IUP shared lib.

function IupPlotOpen()
  ccall( (:IupPlotOpen, iup_plot), Void, (), )
end
function IupPlot()
    ccall((:IupPlot, iup_plot), Ptr{Ihandle}, (),)
end
function IupPlotBegin(ih::Ptr{Ihandle}, strXdata::Integer)
  ccall( (:IupPlotBegin, iup_plot), Void, (Ptr{Ihandle}, Cint), ih, strXdata)
end
function IupPlotAdd(ih::Ptr{Ihandle}, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotAdd, iup_plot), Void, (Ptr{Ihandle}, Cdouble, Cdouble), ih, x, y)
end
function IupPlotAddStr(ih::Ptr{Ihandle}, x::String, y::Cdouble)
  ccall( (:IupPlotAddStr, iup_plot), Void, (Ptr{Ihandle}, Ptr{UInt8}, Cdouble), ih, x, y)
end
function IupPlotEnd(ih::Ptr{Ihandle})
  ccall( (:IupPlotEnd, iup_plot), Cint, (Ptr{Ihandle},), ih)
end
function IupPlotLoadData(ih::Ptr{Ihandle}, filename::Ptr{UInt8}, strXdata::Integer)
  ccall( (:IupPlotLoadData, iup_plot), Cint, (Ptr{Ihandle}, Ptr{UInt8}, Cint), ih, filename, strXdata)
end
function IupPlotInsertStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{UInt8}, y::Cdouble)
  ccall( (:IupPlotInsertStr, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{UInt8}, Cdouble), ih, ds_index, sample_index, x, y)
end
function IupPlotInsert(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotInsert, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble), ih, ds_index, sample_index, x, y)
end
function IupPlotInsertStrPoints(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotInsertStrPoints, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}, Cint), ih, ds_index, sample_index, x, y, count)
end
function IupPlotInsertPoints(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotInsertPoints, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), ih, ds_index, sample_index, x, y, count)
end
function IupPlotAddPoints(ih::Ptr{Ihandle}, ds_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotAddPoints, iup_plot), Void, (Ptr{Ihandle}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), ih, ds_index, x, y, count)
end
function IupPlotAddStrPoints(ih::Ptr{Ihandle}, ds_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble}, count::Integer)
  ccall( (:IupPlotAddStrPoints, iup_plot), Void, (Ptr{Ihandle}, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}, Cint), ih, ds_index, x, y, count)
end
function IupPlotGetSample(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:IupPlotGetSample, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), ih, ds_index, sample_index, x, y)
end
function IupPlotGetSampleStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{Ptr{UInt8}}, y::Ptr{Cdouble})
  ccall( (:IupPlotGetSampleStr, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{Ptr{UInt8}}, Ptr{Cdouble}), ih, ds_index, sample_index, x, y)
end
function IupPlotGetSampleSelection(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer)
  ccall( (:IupPlotGetSampleSelection, iup_plot), Cint, (Ptr{Ihandle}, Cint, Cint), ih, ds_index, sample_index)
end
function IupPlotSetSample(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Cdouble, y::Cdouble)
  ccall( (:IupPlotSetSample, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Cdouble, Cdouble), ih, ds_index, sample_index, x, y)
end
function IupPlotSetSampleStr(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, x::Ptr{UInt8}, y::Cdouble)
  ccall( (:IupPlotSetSampleStr, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Ptr{UInt8}, Cdouble), ih, ds_index, sample_index, x, y)
end
function IupPlotSetSampleSelection(ih::Ptr{Ihandle}, ds_index::Integer, sample_index::Integer, selected::Integer)
  ccall( (:IupPlotSetSampleSelection, iup_plot), Void, (Ptr{Ihandle}, Cint, Cint, Cint), ih, ds_index, sample_index, selected)
end
function IupPlotTransform(ih::Ptr{Ihandle}, x::Cdouble, y::Cdouble, cnv_x::Ptr{Cdouble}, cnv_y::Ptr{Cdouble})
    ccall( (:IupPlotTransform, iup_plot), Void, (Ptr{Ihandle}, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), ih, x, y, cnv_x, cnv_y)
end
function IupPlotTransformTo(ih::Ptr{Ihandle}, cnv_x::Cdouble, cnv_y::Cdouble, x::Ptr{Cdouble}, y::Ptr{Cdouble})
    ccall( (:IupPlotTransformTo, iup_plot), Void, (Ptr{Ihandle}, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), ih, cnv_x, cnv_y, x, y)
end
# COMMENTED because cdCanvas is not known here
#function IupPlotPaintTo(ih::Ptr{Ihandle}, cnv::Ptr{cdCanvas})
#  ccall( (:IupPlotPaintTo, iup_plot), Void, (Ptr{Ihandle}, Ptr{cdCanvas}), ih, cnv)
#end


text:
struct MiniLogger{IOT1 <: IO, IOT2 <: IO, DFT <: DateFormat} <: AbstractLogger
    io::IOT1
    ioerr::IOT2
    errlevel::LogLevel
    minlevel::LogLevel
    message_limits::Dict{Any,Int}
    flush::Bool
    format::Vector{Token}
    dtformat::DFT
end

getio(io) = io
getio(io::AbstractString) = open(io)

function MiniLogger(; io = stdout, ioerr = stderr, errlevel = Error, minlevel = Info, message_limits = Dict{Any, Int}(), flush = true, format = "{[{datetime}]:func} {message}", dtformat = dateformat"yyyy-mm-dd HH:MM:SS")
    tio = getio(io)
    tioerr = io == ioerr ? tio : getio(ioerr)
    MiniLogger(tio, tioerr, errlevel, minlevel, message_limits, flush, tokenize(format), dtformat)
end

shouldlog(logger::MiniLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::MiniLogger) = logger.minlevel

catch_exceptions(logger::MiniLogger) = true

# Formatting of values in key value pairs
showvalue(io, msg) = show(io, "text/plain", msg)
function showvalue(io, e::Tuple{Exception,Any})
    ex, bt = e
    Base.showerror(io, ex, bt; backtrace = bt!==nothing)
end
showvalue(io, ex::Exception) = Base.showerror(io, ex)
showvalue(io, ex::AbstractVector{Union{Ptr{Nothing}, Base.InterpreterIP}}) = Base.show_backtrace(io, ex)

function showmessage(io, msg)
    msglines = split(chomp(string(msg)), '\n')
    print(io, msglines[1])
    for i in 2:length(msglines)
        print(io, " ", msglines[i])
    end
end
showmessage(io, e::Tuple{Exception,Any}) = showvalue(io, e)
showmessage(io, ex::Exception) = showvalue(io, e)
showmessage(io, ex::AbstractVector{Union{Ptr{Nothing}, Base.InterpreterIP}}) = Base.show_backtrace(io, ex)

tsnow(dtf) = Dates.format(Dates.now(), dtf)

function colorfunc(level::LogLevel, _module, group, id, filepath, line, element)
    level < Info  ? Color(Base.debug_color()) :
    level < Warn  ? Color(Base.info_color())  :
    level < Error ? Color(Base.warn_color())  :
                    Color(Base.error_color())
end

function extractcolor(token::Token, level, _module, group, id, filepath, line)
    if token.c.c == :func
        color = colorfunc(level, _module, group, id, filepath, line, token.val)
        return Color(color.c, color.isbold || token.c.isbold)
    else
        return token.c
    end
end

function printwcolor(iob, val, color)
    if color.c == -1
        if color.isbold
            printstyled(iob, val; bold = color.isbold)
        else
            print(iob, val)
        end
    else
        printstyled(iob, val; color = color.c, bold = color.isbold)
    end
end

function handle_message(logger::MiniLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    io = if level < logger.errlevel
            isopen(logger.io) ? logger.io : stdout
        else
            isopen(logger.ioerr) ? logger.ioerr : stderr
        end

    buf = IOBuffer()
    iob = IOContext(buf, io)

    for token in logger.format
        c = extractcolor(token, level, _module, group, id, filepath, line)
        val = token.val
        if val == "datetime"
            printwcolor(iob, tsnow(logger.dtformat), c)
        elseif val == "message"
            showmessage(iob, message)

            iscomma = false
            for (key, val) in kwargs
                if string(key) == val
                    print(iob, key)
                    iscomma = false
                else
                    iscomma && print(iob, ", ")
                    print(iob, key, " = ")
                    showvalue(iob, val)
                    iscomma = true
                end
            end
        elseif val == "level"
            levelstr = string(level)
            printwcolor(iob, levelstr, c)
        elseif val == "basename"
            printwcolor(iob, basename(filepath), c)
        elseif val == "filepath"
            printwcolor(iob, filepath, c)
        elseif val == "line"
            printwcolor(iob, line, c)
        elseif val == "group"
            printwcolor(iob, group, c)
        elseif val == "module"
            printwcolor(iob, _module, c)
        elseif val == "id"
            printwcolor(iob, id, c)
        else
            printwcolor(iob, val, c)
        end
    end
    print(iob, "\n")
    write(io, take!(buf))
    if logger.flush
        flush(io)
    end
    nothing
end


text:
using OrdinaryDiffEq, DiffEqBase, Test

## Start on Number
f = (u,p,t) -> (2u)
#(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = u0*exp(t)
prob = ODEProblem(f,1/2+(1/4)im,(0.0,1.0))

sol = solve(prob,Tsit5(),dt=1/2^4)

u0 = rand(ComplexF32,5,5,5)
f = (du,u,p,t) -> (du.=2u)
#(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = u0*exp(t)
prob = ODEProblem(f,u0,(0.0,1.0))

sol = solve(prob,Tsit5(),dt=1/2^4)

@test typeof(sol[1]) == Array{ComplexF32,3}


text:
using MatroidChowRings
using Test

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

max_stars_repo_name: FGlazov/ChowRings

max_issues_repo_name: FGlazov/ChowRings

max_forks_repo_name: FGlazov/ChowRings


text:
export helixnodes

hexlixpoint(radius,pitch,angle)=(radius*cos(angle), radius*sin(angle), pitch*angle/(2π))

"""
        helixnodes(radius, pitch, radians, [radiansperpoint=π/4])

Returns an array of `Node`s in a helix.
"""
function helixnodes(radius::Float64, pitch::Float64, radians::Float64, radiansperpoint::Float64=π/4)
    numberofpoints = Int(cld(radians,radiansperpoint))+1
    nodes = Array{Node}(undef,numberofpoints)
    i = 0
    for angle in 0:radiansperpoint:radians
        i+=1
        nodes[i] = Node(hexlixpoint(radius,pitch,angle)...)
    end
    if i<numberofpoints # end point didn't fall on a multiple of radiansperpoint
        i+=1
        nodes[i] = Node(hexlixpoint(radius,pitch,radians)...)
      end
    return nodes
end


text:
# doc.jl - api

"""
    create_journal_entry(::String, ::String, ::String, ::Int64, ::Int64, ::Float64, ::Float64, ::Float64, ::String)::JournalEntry

    It creates a general ledger `JournalEntry.`
.
# Example

```
julia> entry = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:25:29.277, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
```
"""
function create_journal_entry end

"""
    create_account(::String, ::String, ::String)::Account

# Example

```
julia> account = create_account("8000", "Sales", "")
Account("8000", "Sales", "")
```
"""
function create_account end

"""
    create_record_from(entry::JournalEntry)::Record

It creates a general ledger `Record` from a `JournalEntry.`

@see also [`add_to_file`](@ref), [`read_from_file`](@ref)

# Example

```
julia> entry = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:28:14.696, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> record_from = create_record_from(entry)
Record("20200700", 1300, 2020-03-03T13:28:14.696, "APPLI", "A1001", 1210.0, 0.0, "LS")

julia> add_to_file("./test_ledger.txt", [record_from])
```
"""
function create_record_from end

"""
    create_record_to(::JournalEntry)::Record

It creates a general ledger `Record` from a `JournalEntry.`

@see also [`add_to_file`](@ref)

# Example
```
julia> entry = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:28:54.701, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> record_to = create_record_to(entry)
Record("20200700", 8000, 2020-03-03T13:28:54.701, "APPLI", "A1001", 0.0, 1000.0, "LS")

julia> add_to_file("./test_ledger.txt", [record_to])
```
"""
function create_record_to end

"""
    create_record_vat(::JournalEntry, ::Int)::Record

It creates a general ledger `Record` from a `JournalEntry.`

@see also [`add_to_file`](@ref)

# Example

```
julia> entry = create_journal_entry("20200700", "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")
JournalEntry("20200700", 3, 2020-03-03T13:29:41.847, "APPLI", "A1001", 1300, 8000, 1000.0, 0.0, 210.0, "LS")

julia> record_vat = create_record_vat(entry, 4000)
Record("20200700", 4000, 2020-03-03T13:29:41.847, "APPLI", "A1001", 0.0, 210.0, "LS")

julia> add_to_file("./test_ledger.txt", [record_vat])

```
"""
function create_record_vat end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Cbc
using JuMP
using MathOptInterface
using MIPLearn
const MOI = MathOptInterface

function find_lazy(model::Model, cb_data)::Vector{String}
    x = variable_by_name(model, "x")
    y = variable_by_name(model, "y")
    x_val = value(cb_data, x)
    y_val = value(cb_data, y)
    if x_val + y_val > 1 + 1e-6
        return ["con"]
    end
    return []
end

function enforce_lazy(model::Model, cb_data, violation::String)::Nothing
    if violation == "con"
        x = variable_by_name(model, "x")
        y = variable_by_name(model, "y")
        con = @build_constraint(x + y <= 1)
        submit(cb_data, con)
    end
    return
end

function build_model(data)
    model = Model()
    @variable(model, x, Bin)
    @variable(model, y, Bin)
    @objective(model, Max, 2 * x + y)
    @constraint(model, c1, x + y <= 2)
    @lazycb(model, find_lazy, enforce_lazy)
    return model
end

@testset "Lazy callback" begin
    @testset "JuMPInstance" begin
        model = build_model(nothing)
        instance = JuMPInstance(model)
        solver = LearningSolver(Cbc.Optimizer)
        solve!(solver, instance)
        @test value(model[:x]) == 1.0
        @test value(model[:y]) == 0.0
    end

    @testset "FileInstance" begin
        data = nothing
        filename = tempname()
        MIPLearn.save_data(filename, data)
        instance = FileInstance(filename, build_model)
        solver = LearningSolver(Cbc.Optimizer)
        solve!(solver, instance)
        h5 = MIPLearn.Hdf5Sample("$filename.h5")
        @test h5.get_array("mip_var_values") == [1.0, 0.0]
    end
end


text:
"""
Find root of h(x) in the interval [a, b] using binary search
"""
function bisection(h, a, b, tol=1e-5, nmax=10)
    #find root in [a, b]
    n=1
    c=0
    while n<nmax
        c = (a+b)/2
        if abs(h(c))<tol
            return c
        elseif (sign(h(c))==sign(h(a)))
            a=c
        else
            b=c
        end
        n=n+1
    end
    return c
end

"""
Computes the symmetric 95% confidence interval about the median of a
probability distribution given by its PDF using bisection
"""
function confidence(pdf, median, alpha=0.95)
    f = h -> int1D(pdf, median-h, median+h, "2") - alpha
    c = bisection(f, 0, median, 1e-3, 10)
    if abs(c)>median*.9 #this should mean that integration scheme failed to detect spike, indicating that interval is very thin to begin with
        c=.1
    end
    return c
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Convex
using HDF5, JLD

folder = "greedy_vs_random"

#using CenterOfMass
include("CenterOfMass.jl")
include("setup_parameters.jl")

n_trial = 1000
n_measurement = 20

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

@time critical_forces_by_point = map(attachment_ws) do wrench
  get_critical_values(circle_ws, wrench, prior, interior_q)
end

@time normals_matrices = map(critical_forces_by_point) do forces
  normal_matrix(forces[:], sigma^2)
end

rand_com() = r_attach * rand_in_circle()[1:2] # hack
rand_mass() = min_mass + rand() * (max_mass - min_mass)

rand_parameters() = [rand_com(); rand_mass()]

thetas = hcat(map(x->rand_parameters(), 1:n_trial)...)

select_random(prior) = (rand(1:length(attachment_ps)), ones(length(attachment_ps)))

cycle_val = 1
select_cyclic(prior) = begin
  global cycle_val
  val = cycle_val
  cycle_val = (cycle_val-1+3)%length(attachment_ps) + 1
  (val, ones(length(attachment_ps)))
end

function select_csqmi(prior)
  n = length(critical_forces_by_point)
  csqmis = map(1:n) do ii
    compute_mutual_information(get_data(prior)[:], normals_matrices[ii])
  end

  indmax(csqmis), csqmis
end

selection_functions = [select_csqmi, select_random, select_cyclic]
nf = length(selection_functions)

type TrialData
  #belief::Histogram
  estimate
  error
  applied_f
  f_hat
  csqmis
  robot_measurement
end

data_array = Array{TrialData}(n_trial, nf, n_measurement)

belief = deepcopy(prior)

for ii = 1:n_trial
  theta = thetas[:,ii]
  com_p = [theta[1:2];0]
  mass = theta[3]

  println("trial $(ii)")

  for jj = 1:nf
    println("function $(jj)")
    selection_function = selection_functions[jj]
    belief = deepcopy(prior)

    for kk = 1:n_measurement
      index, csqmis = selection_function(belief)
      applied_p = attachment_ps[index]

      applied_f = critical_force_from_points(circle_ps, com_p, applied_p, mass)

      f_hat = applied_f + sigma * randn()

      critical_forces = critical_forces_by_point[index]

      update_prior!(belief, critical_forces, f_hat, sigma)

      estimate = weighted_average(belief)
      error = normalized_error(estimate - theta)

      data = TrialData(deepcopy(estimate), error, applied_f,
        f_hat, csqmis, index)
      data_array[ii,jj,kk] = data
    end
  end
end

if ~isdir(folder)
  mkdir(folder)
end
@save "$(folder)/data_new" data_array thetas
Void

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
insts = OrderedDict()
insts["minimal"] = [
    ((1, 2, 2, 2),),
    ]
insts["fast"] = [
    ((2, 2, 3, 2),),
    ((3, 3, 3, 3),),
    ((3, 3, 5, 4),),
    ((5, 2, 5, 3),),
    ((1, 30, 2, 30),),
    ((10, 1, 3, 1),),
    ]
insts["various"] = [
    ((3, 3, 3, 3),),
    ((3, 3, 5, 4),),
    ((5, 2, 5, 3),),
    ((2, 30, 3, 30),),
    ((2, 15, 6, 15),),
    ((2, 30, 2, 40),),
    ((8, 3, 3, 3),),
    ]
return (PolyEnvelopeJuMP, insts)


text:
# Spell and einstellgrund encodings in the AMP data, eg the `hv_epi_uni_roh` file.
#
# Available as a submodule, use as eg `AMP.employee`.

module AMP

import ..@enum_and_Dict

@enum_and_Dict(Spell,
               spell_dict,
               "D2" => DLU,
               "BE" => civil_servant,
               "LE" => apprenticeship,
               "AA" => employee,
               "FU" => FUB,
               "FD" => service_contract,
               "SO" => other_employment,
               "LW" => farmer,
               "S1" => self_employed,
               "AL" => unemployed_mixed, # with and without benefits
               "AO" => unemployed_nobenefits,
               "W1" => maternity_active,
               "W2" => maternity_inactive,
               "EO" => parental_leave_inactive,
               "ED" => parental_leave_active,
               "KO" => child_care_allowance_inactive,
               "KG" => child_care_allowance_active,
               "PZ" => military,
               "RE" => retired,
               "AU" => education,
               "SG" => other_insured_nonemployment,
               "G1" => minor_employment,
               "SV" => other_insured_time,
               "LL" => data_gap,
               "GT" => birth,
               "TO" => death,
               "KD" => no_data,
               "66" => transition_allowance,
               "68" => rehabilitation)

end # module


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=3.25)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
# Note, this routine is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_discontinuous_well_balancedness` below.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
solver = DGSEM(polydeg=4, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the TreeMesh and setup a periodic mesh

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solver

tspan = (0.0, 100.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.

# alternative version of the initial conditinon used to setup a truly discontinuous
# bottom topography function for this academic testcase of well-balancedness.
# The errors from the analysis callback are not important but the error for this lake at rest test case
# `∑|H0-(h+b)|` should be around machine roundoff
# In contrast to the usual signature of initial conditions, this one get passed the
# `element_id` explicitly. In particular, this initial conditions works as intended
# only for the TreeMesh2D with initial_refinement_level=2.
function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  b = 0.0

  # Setup a discontinuous bottom topography using the element id number
  if element_id == 7
    b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2])
  end

  return prim2cons(SVector(H, v1, v2, b), equations)
end

# point to the data we want to augment
u = Trixi.wrap_array(ode.u0, semi)
# reset the initial condition
for element in eachelement(semi.solver, semi.cache)
  for j in eachnode(semi.solver), i in eachnode(semi.solver)
    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, j, element)
    u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations)
    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
  end
end

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=3.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using CthulhuVision.Image
using CthulhuVision.Math
using CthulhuVision.Light
using CthulhuVision.Materials
using CthulhuVision.Rendering
using CthulhuVision.Camera
using CthulhuVision.Triangles
using CthulhuVision.Scenes
using CthulhuVision.Random

################
# Define image #
################

width = 512
height = 512
image = PPM(Dimension(width, height))

#################
# SceneSettings #
#################

ambientemission = RGB(0.5f0, 0.7f0, 1.0f0)
settings = SceneSettings(ambientemission)

##########
# Camera #
##########

aspect = Float32(image.dimension.width / image.dimension.height)
vfov = 40.0f0

lookfrom = Vector3(0f0, 0f0, -25f0)
lookat = Vector3(0f0, 0f0, 0f0)
vup = Vector3(0f0, 1f0, 0f0)
focusdist = 10.0f0
aperture = 0.0f0

camera = FovCamera(lookfrom, lookat, vup, vfov, aspect, aperture, focusdist)

###################
# Construct scene #
###################

function constructscene(angle::Float32) :: Scene
    r = rotation(angle, Vector3(1f0, 0f0, 0f0))

    red = lambertian(RGB(1f0, 0f0, 0f0))
    triangle = Triangle(Vector3(0f0, 0f0, 0f0), Vector3(1f0, 0f0, 0f0), Vector3(0f0, 1f0, 0f0), red)

    t = transform(triangle, r)

    Scene([t], settings)
end

##################
# Perform render #
##################

rendersettings = RenderSettings(10)

lastframe = 1

for frame = 0:lastframe
    angle = Float32(2f0 * π) * 3f0 * Float32(frame) / Float32(lastframe)
    scene = constructscene(angle)
    render(image, camera, scene, rendersettings)
    framestr = lpad(frame, 4, '0')
    saveimage(image, examplefilename("spinningtriangle_$(framestr)"))
end


text:
using NLSolvers
using Test
using StaticArrays
#using Optim
#using LineSearches
using Printf
using LinearAlgebra: norm, I
import Random
Random.seed!(41234)
include("optimize/problems.jl")
include("optimize/interface.jl")
include("optimize/preconditioning.jl")
include("optimize/complex.jl")
include("nlsolve/problems.jl")
include("nlsolve/interface.jl")


text:
using DielectricUpdateTechique
using Test

@testset "DielectricUpdateTechique.jl" begin
    # Write your tests here.
    include("testRemainingUtils.jl")
    include("testGetConstants.jl")
    include("testCellToYeeDielectric.jl")
    include("testComputeUpdateMaps.jl")
    include("testOperators.jl")
    include("testInit.jl")
    include("testSolver.jl")
end


text:
using Documenter
using FoldsThreads

makedocs(
    sitename = "FoldsThreads",
    format = Documenter.HTML(),
    modules = [FoldsThreads]
)

deploydocs(; repo = "github.com/JuliaFolds/FoldsThreads.jl", push_preview = true)


text:
using Pope: LJH
using Base.Test
using ReferenceMicrocalFiles

@testset "matter simulator" begin

src = ReferenceMicrocalFiles.dict["good_mnka_mystery"].filename
dest = tempdir()
timeout_s=0.0003
maxchannels = 4
fastforward = 1

sim = @schedule Pope.mattersim(src, dest, timeout_s, fastforward, maxchannels)
sleep(0.1)
filename0, writingbool0 = LJH.matter_writing_status()
wait(sim)
filename1, writingbool1 = LJH.matter_writing_status()

# these two use to pass, but stopped working for reasons I don't understand
# when I made seemingly unrelated changes
# @test writingbool0
# @test !writingbool1
@test filename0==filename1
@test isfile(filename0)
end


text:
using LinearAlgebra: lu
# FIXME: This doesn't really work if x is anything but a Vector or a scalar
function fdjacobian(F, x::Number, t)
    ftx = F(t, x)

    # The 100 below is heuristic
    dx = (x .+ (x==0))./100
    dFdx = (F(t,x+dx)-ftx)./dx

    return dFdx
end

function fdjacobian(F, x, t)
    ftx = F(t, x)
    lx = max(length(x),1)
    dFdx = zeros(eltype(x), lx, lx)
    for j = 1:lx
        # The 100 below is heuristic
        dx = zeros(eltype(x), lx)
        dx[j] = (x[j] .+ (x[j]==0))./100
        dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]
    end
    return dFdx
end

# ODE23S  Solve stiff systems based on a modified Rosenbrock triple
# (also used by MATLABS ODE23s); see Sec. 4.1 in
#
# [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22
#
# supports keywords: points = :all | :specified (using dense output)
#                    jacobian = G(t,y)::Function | nothing (FD)
#                    mass = M::Matrix | nothing (I)
function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
                                                jacobian=nothing,
                                                mass=nothing,
                                                points=:all,
                                                norm=Base.norm)

    # select method for computing the Jacobian
    if typeof(jacobian) == Function
        jac = jacobian
    else
        # fallback finite-difference
        jac = (t, y)->fdjacobian(F, y, t)
    end
    if mass == nothing
        M = I
    else
        M = mass
    end

    # constants
    d = 1/(2 + sqrt(2))
    e32 = 6 + sqrt(2)

    # initialization
    t = tspan[1]
    tfinal = tspan[end]
    tdir = sign(tfinal - t)

    hmax = abs(tfinal - t)/10
    hmin = abs(tfinal - t)/1e9
    h = tdir*abs(tfinal - t)/100  # initial step size

    y = y0
    tout = Array{typeof(t)}(undef, 1)
    tout[1] = t         # first output time
    yout = Array{typeof(y0)}(undef, 1)
    yout[1] = copy(y)         # first output solution

    F0 = F(t,y)

    J = jac(t,y)    # get Jacobian of F wrt y

    while abs(t) < abs(tfinal) && hmin < abs(h)
        if abs(t-tfinal) < abs(h)
            h = tfinal - t
        end

        if size(J,1) == 1
            W = M - h*d*J
        else
            # (see Sec. 5 in [SR97])
            W = lu( M - h*d*J )
        end

        # approximate time-derivative of F
        T = h*d*(F(t + h/100, y) - F0)/(h/100)

        # modified Rosenbrock formula
        k1 = W\(F0 + T)
        F1 = F(t + 0.5*h, y + 0.5*h*k1)
        k2 = W\(F1 - M*k1) + k1
        ynew = y + h*k2
        F2 = F(t + h, ynew)
        k3 = W\(F2 - e32*(M*k2 - F1) - 2*(M*k1 - F0) + T )

        err = (h/6)*norm(k1 - 2*k2 + k3) # error estimate
        delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error

        # check if new solution is acceptable
        if  err <= delta

            if points==:specified || points==:all
                # only points in tspan are requested
                # -> find relevant points in (t,t+h]
                for toi in tspan[(tspan.>t) .& (tspan.<=t+h)]
                    # rescale to (0,1]
                    s = (toi-t)/h

                    # use interpolation formula to get solutions at t=toi
                    push!(tout, toi)
                    push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d)))
                end
            end
            if (points==:all) && (tout[end]!=t+h)
                # add the intermediate points
                push!(tout, t + h)
                push!(yout, ynew)
            end

            # update solution
            t = t + h
            y = ynew

            F0 = F2         # use FSAL property
            J = jac(t,y)    # get Jacobian of F wrt y
                            # for new solution
        end

        # update of the step size
        h = tdir*min( hmax, abs(h)*0.8*(delta/err)^(1/3) )
    end

    return tout, yout
end


text:
module DWDDynam

using Catlab.Theories
using Catlab.WiringDiagrams.DirectedWiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
import Catlab.WiringDiagrams: oapply, input_ports, output_ports

import ..UWDDynam: nstates, eval_dynamics, euler_approx, AbstractInterface

using OrdinaryDiffEq, DynamicalSystems, DelayDiffEq
import OrdinaryDiffEq: ODEProblem
import DelayDiffEq: DDEProblem
import DynamicalSystems: DiscreteDynamicalSystem
using Plots

export AbstractMachine, ContinuousMachine, DiscreteMachine, DelayMachine,
nstates, ninputs, noutputs, eval_dynamics, readout, euler_approx

using Base.Iterators
import Base: show, eltype

### Interface
abstract type AbstractDirectedInterface{T} <: AbstractInterface{T} end

struct DirectedInterface{T} <: AbstractDirectedInterface{T}
    input_ports::Vector
    output_ports::Vector
end
DirectedInterface{T}(ninputs::Int, noutputs::Int) where T =
    DirectedInterface{T}(1:ninputs, 1:noutputs)

struct DirectedVectorInterface{T,N} <: AbstractDirectedInterface{T}
    input_ports::Vector
    output_ports::Vector
end
DirectedVectorInterface{T,N}(ninputs::Int, noutputs::Int) where {T,N} =
    DirectedVectorInterface{T,N}(1:ninputs, 1:noutputs)

input_ports(interface::AbstractDirectedInterface) = interface.input_ports
output_ports(interface::AbstractDirectedInterface) = interface.output_ports
ninputs(interface::AbstractDirectedInterface) = length(input_ports(interface))
noutputs(interface::AbstractDirectedInterface) = length(output_ports(interface))

ndims(::DirectedVectorInterface{T, N}) where {T,N} = N

unit(::Type{I}, dims) where {T, I<:AbstractDirectedInterface{T}} = zeros(T, dims)
unit(::Type{DirectedVectorInterface{T,N}}, dims) where {T, N} = fill(zeros(T, N), dims)

### Dynamics
abstract type AbstractDirectedSystem{T} end

struct ContinuousDirectedSystem{T} <: AbstractDirectedSystem{T}
    nstates::Int
    dynamics::Function
    readout::Function
end

struct DiscreteDirectedSystem{T} <: AbstractDirectedSystem{T}
    nstates::Int
    dynamics::Function
    readout::Function
end

struct DelayDirectedSystem{T} <: AbstractDirectedSystem{T}
    nstates::Int
    dynamics::Function
    readout::Function
end

nstates(dynam::AbstractDirectedSystem) = dynam.nstates
dynamics(dynam::AbstractDirectedSystem) = dynam.dynamics
readout(dynam::AbstractDirectedSystem) = dynam.readout

"""

A directed open dynamical system operating on information fo type `T`.
A machine  `m` has type signature  (`m.ninputs`, `m.outputs`).
"""
abstract type AbstractMachine{T, InterfaceType, SystemType} end

interface(m::AbstractMachine) = m.interface
system(m::AbstractMachine) = m.system

input_ports(m::AbstractMachine) = input_ports(interface(m))
output_ports(m::AbstractMachine) = output_ports(interface(m))
ninputs(m::AbstractMachine) = ninputs(interface(m))
noutputs(m::AbstractMachine) = noutputs(interface(m))
nstates(m::AbstractMachine) = nstates(system(m))
dynamics(m::AbstractMachine) = dynamics(system(m))
readout(m::AbstractMachine) = readout(system(m))

struct Machine{T,I,S} <: AbstractMachine{T,I,S}
    interface::I
    system::S
end

"""    ContinuousMachine{T}(ninputs, nstates, noutputs, f, r)

An directed open continuous system. The dynamics function `f` defines an ODE ``\\dot u(t) = f(u(t),x(t),p,t)`` where ``u`` is the state and ``x`` captures the exogenous variables.

The readout function may depend on the state, parameters, and time, so it must be of the form ``r(u,p,t)``.
"""
const ContinuousMachine{T,I} = Machine{T, I, ContinuousDirectedSystem{T}}

ContinuousMachine{T}(interface::I, system::ContinuousDirectedSystem{T}) where {T, I <: AbstractDirectedInterface} =
    ContinuousMachine{T, I}(interface, system)

ContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout) where T =
    ContinuousMachine{T}(DirectedInterface{T}(ninputs, noutputs), ContinuousDirectedSystem{T}(nstates, dynamics, readout))

ContinuousMachine{T}(ninputs::Int, nstates::Int, dynamics) where T  =
    ContinuousMachine{T}(ninputs, nstates, nstates, dynamics, (u,p,t) -> u)

ContinuousMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout) where {T,I <: AbstractDirectedInterface} =
    ContinuousMachine{T,I}(I(ninputs, noutputs), ContinuousDirectedSystem{T}(nstates, dynamics, readout))

ContinuousMachine{T, N}(ninputs, nstates, noutputs, dynhamics, readout) where {T,N} =
    ContinuousMachine{T, DirectedVectorInterface{T, N}}(ninputs, nstates, noutputs, dynhamics, readout)

ContinuousMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I}  =
    ContinuousMachine{T,I}(ninputs, nstates, nstates, dynamics, (u,p,t) -> u)

"""    DelayMachine{T}(ninputs, nstates, noutputs, f, r)

A delay open continuous system. The dynamics function `f` defines an ODE ``\\dot u(t) = f(u(t), x(t), h(p,t), p, t)`` where
``u`` is the states, ``x`` captures the exogenous variables, and ``h`` is a history function

The readout function may depend on the state, history, parameters, and time, so it has a signature ``r(u,h,p,t)``.
"""
const DelayMachine{T,I} = Machine{T, I, DelayDirectedSystem{T}}

DelayMachine{T}(interface::I, system::DelayDirectedSystem{T}) where {T, I<:AbstractDirectedInterface} =
    DelayMachine{T, I}(interface, system)

DelayMachine{T}(ninputs, nstates, noutputs, dynamics, readout) where T =
    DelayMachine{T}(DirectedInterface{T}(ninputs, noutputs), DelayDirectedSystem{T}(nstates, dynamics, readout))

DelayMachine{T}(ninputs::Int, nstates::Int, dynamics) where T  =
    DelayMachine{T}(ninputs, nstates, nstates, dynamics, (u,p,t) -> u)

DelayMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout) where {T,I<:AbstractDirectedInterface} =
    DelayMachine{T,I}(I(ninputs, noutputs), DelayDirectedSystem{T}(nstates, dynamics, readout))

DelayMachine{T, N}(ninputs, nstates, noutputs, dynhamics, readout) where {T,N} =
    DelayMachine{T, DirectedVectorInterface{T, N}}(ninputs, nstates, noutputs, dynhamics, readout)

DelayMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I}  =
    DelayMachine{T,I}(ninputs, nstates, nstates, dynamics, (u,h,p,t) -> u)

"""    DiscreteMachine{T}(ninputs, nstates, noutputs, f, r)

A directed open discrete dynamical system. The dynamics function `f` defines a discrete update rule ``u_{n+1} = f(u_n, x_n, p, t)`` where ``u_n`` is the state and ``x_n`` is the value of the exogenous variables at the ``n``th time step.

The readout function may depend on the state, parameters, and time step, so it must be of the form ``r(u_n,p,n)``.
"""
const DiscreteMachine{T,I} = Machine{T, I, DiscreteDirectedSystem{T}}

DiscreteMachine{T}(interface::I, system::DiscreteDirectedSystem{T}) where {T, I<:AbstractDirectedInterface} =
    DiscreteMachine{T, I}(interface, system)

DiscreteMachine{T}(ninputs, nstates, noutputs, dynamics, readout) where T =
    DiscreteMachine{T}(DirectedInterface{T}(ninputs, noutputs), DiscreteDirectedSystem{T}(nstates, dynamics, readout))

DiscreteMachine{T}(ninputs::Int, nstates::Int, dynamics) where T  =
    DiscreteMachine{T}(ninputs, nstates, nstates, dynamics, (u,p,t) -> u)

DiscreteMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout) where {T,I<:AbstractDirectedInterface} =
    DiscreteMachine{T,I}(I(ninputs, noutputs), DiscreteDirectedSystem{T}(nstates, dynamics, readout))

DiscreteMachine{T, N}(ninputs, nstates, noutputs, dynhamics, readout) where {T,N} =
    DiscreteMachine{T, VectorInterface{T, N}}(ninputs, nstates, noutputs, dynhamics, readout)

DiscreteMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I}  =
    DiscreteMachine{T,I}(ninputs, nstates, nstates, dynamics, (u,p,t) -> u)

show(io::IO, vf::ContinuousMachine) = print("ContinuousMachine(ℝ^$(nstates(vf)) × ℝ^$(ninputs(vf)) → ℝ^$(nstates(vf)))")
show(io::IO, vf::DelayMachine) = print("DelayMachine(ℝ^$(nstates(vf)) × ℝ^$(ninputs(vf)) → ℝ^$(nstates(vf)))")
show(io::IO, vf::DiscreteMachine) = print("DiscreteMachine(ℝ^$(nstates(vf)) × ℝ^$(ninputs(vf)) → ℝ^$(nstates(vf)))")

eltype(::AbstractMachine{T}) where T = T

readout(f::DelayMachine, u::AbstractVector, h = nothing, p = nothing, t = 0) = readout(f)(u, h, p, t)
readout(f::AbstractMachine, u::AbstractVector, p = nothing, t = 0) = readout(f)(u, p, t)
readout(f::AbstractMachine, u::FinDomFunction, args...) = readout(f, collect(u), args...)

"""    eval_dynamics(m::AbstractMachine, u::AbstractVector, xs:AbstractVector, p, t)

Evaluates the dynamics of the machine `m` at state `u`, parameters `p`, and time `t`. The exogenous variables are set by `xs` which may either be a collection of functions ``x(t)`` or a collection of constant values.

The length of `xs` must equal the number of inputs to `m`.
"""
eval_dynamics(f::DelayMachine, u::AbstractVector, xs::AbstractVector, h::Function, p=nothing, t=0) = begin
    ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
    dynamics(f)(u, xs, h, p, t)
end

eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector, p=nothing, t=0) = begin
    ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
    dynamics(f)(u, xs, p, t)
end

eval_dynamics(f::AbstractMachine, u::S, xs::T, args...) where {S,T <: Union{FinDomFunction, AbstractVector}} =
    eval_dynamics(f, collect(u), collect(xs), args...)

eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector{T}, p=nothing, t=0) where T <: Function =
    eval_dynamics(f, u, [x(t) for x in xs], p, t)

"""    euler_approx(m::ContinuousMachine, h)

Transforms a continuous machine `m` into a discrete machine via Euler's method with step size `h`. If the dynamics of `m` is given by ``\\dot{u}(t) = f(u(t),x(t),p,t)`` the the dynamics of the new discrete system is given by the update rule ``u_{n+1} = u_n + h f(u_n, x_n, p, t)``.
"""
euler_approx(f::ContinuousMachine{T}, h::Float64) where T = DiscreteMachine{T}(
    ninputs(f), nstates(f), noutputs(f),
    (u, x, p, t) -> u + h*eval_dynamics(f, u, x, p, t),
    readout(f)
)

"""    euler_approx(m::ContinuousMachine)

Transforms a continuous machine `m` into a discrete machine via Euler's method where the step size is introduced as a new parameter, the last in the list of parameters.
"""
euler_approx(f::ContinuousMachine{T}) where T = DiscreteMachine{T}(
    ninputs(f), nstates(f), noutputs(f),
    (u, x, p, t) -> u + p[end]*eval_dynamics(f, u, x, p[1:end-1], t),
    readout(f)
)
euler_approx(fs::Vector{M}, args...) where {M<:ContinuousMachine} =
    map(f->euler_approx(f,args...), fs)

euler_approx(fs::AbstractDict{S, M}, args...) where {S, M<:ContinuousMachine} =
    Dict(name => euler_approx(f, args...) for (name, f) in fs)

# Integration with ODEProblem in OrdinaryDiffEq.jl

"""    ODEProblem(m::ContinuousMachine, xs::Vector, u0::Vector, tspan)

Constructs an ODEProblem from the vector field defined by `(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`.
"""
ODEProblem(m::ContinuousMachine, u0::AbstractVector, xs::AbstractVector, tspan::Tuple{Real, Real}, p=nothing)  =
    ODEProblem((u,p,t) -> eval_dynamics(m, u, xs, p, t), u0, tspan, p)

ODEProblem(m::ContinuousMachine, u0::AbstractVector, x, tspan::Tuple{Real, Real}, p=nothing) =
    ODEProblem(m, u0, collect(repeated(x, ninputs(m))), tspan, p)

ODEProblem(m::ContinuousMachine{T}, u0::AbstractVector, tspan::Tuple{Real, Real}, p=nothing) where T =
    ODEProblem(m, u0, T[], tspan, p)

"""    DDEProblem(m::DelayMachine, u0::Vector, xs::Vector, h::Function, tspan)
"""
DDEProblem(m::DelayMachine, u0::AbstractVector, xs::AbstractVector, h, tspan::Tuple{Real, Real}, p=nothing; kwargs...) =
    DDEProblem((u,h,p,t) -> eval_dynamics(m, u, xs, h, p, t), u0, h, tspan, p; kwargs...)

"""    DiscreteDynamicalSystem(m::DiscreteMachine, xs::Vector, u0::Vector, p)

Constructs an DiscreteDynamicalSystem from the equation of motion defined by
`(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`. Pass `nothing` in place of `p` if your system does not have parameters.
"""
DiscreteDynamicalSystem(m::DiscreteMachine, u0::AbstractVector, x,  p; t0::Int = 0) =
  DiscreteDynamicalSystem(m, u0, collect(repeated(x, ninputs(m))), p; t0=t0)

DiscreteDynamicalSystem(m::DiscreteMachine{T}, u0::AbstractVector, xs::AbstractVector, p; t0::Int = 0) where T = begin
  if nstates(m) == 1
    DiscreteDynamicalSystem1d(m, u0[1], xs, p; t0=t0)
  else
    !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
    DiscreteDynamicalSystem(
        (u,p,t) -> SVector{nstates(m)}(eval_dynamics(m, u, xs, p, t)),
        u0, p; t0=t0)
  end
end

DiscreteDynamicalSystem(m::DiscreteMachine, u0::Real, xs::AbstractVector, p; t0::Int = 0) =
  DiscreteDynamicalSystem1D(m, u0, xs, p; t0=t0)

# if the system is 1D then the state must be represented by a number NOT by a 1D array
DiscreteDynamicalSystem1d(m::DiscreteMachine{T}, u0::Real, xs::AbstractVector, p; t0::Int = 0) where T = begin
  nstates(m) == 1 || error("The machine must have exactly 1 state")
  !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
  DiscreteDynamicalSystem(
      (u,p,t) -> eval_dynamics(m, [u], xs, p, t)[1],
      u0, p; t0=t0)
end

DiscreteDynamicalSystem(m::DiscreteMachine, u0, p; t0::Int = 0) =
    DiscreteDynamicalSystem(m, u0, [], p; t0 = t0)

### Plotting backend
@recipe function f(sol, m::AbstractMachine, p=nothing)
    labels = (String ∘ Symbol).(output_ports(m))
    label --> reshape(labels, 1, length(labels))
    vars --> map(1:noutputs(m)) do i
        ((t, args...) -> (t, readout(m)(collect(args), p, t)[i]), 0:nstates(m)...)
    end
    sol
end

"""    oapply(d::WiringDiagram, ms::Vector)

Implements the operad algebras for directed composition of dynamical systems given a
composition pattern (implemented by a directed wiring diagram `d`)
and primitive systems (implemented by a collection of
machines `ms`).

Each box of the composition pattern `d` is filled by a machine with the
appropriate type signature. Returns the composite machine.
"""
function oapply(d::WiringDiagram, ms::Vector{M}) where {M<:AbstractMachine}
    isempty(wires(d, input_id(d), output_id(d))) || error("d is not a valid composition syntax because it has pass wires")
    nboxes(d) == length(ms)  || error("there are $nboxes(d) boxes but $length(ms) machines")
    for box in 1:nboxes(d)
        fills(ms[box], d, box) || error("$ms[box] does not fill box $box")
    end

    S = coproduct((FinSet∘nstates).(ms))
    Inputs = coproduct((FinSet∘ninputs).(ms))

    return M(input_ports(d),
        length(apex(S)),
        output_ports(d),
        induced_dynamics(d, ms, S, Inputs),
        induced_readout(d, ms, S))
end

"""    oapply(d::WiringDiagram, m::AbstractMachine)

A version of `oapply` where each box of `d` is filled with the machine `m`.
"""
function oapply(d::WiringDiagram, x::AbstractMachine)
    oapply(d, collect(repeated(x, nboxes(d))))
end

"""    oapply(d::WiringDiagram, generators::Dict)

A version of `oapply` where `generators` is a dictionary mapping the name of each box to its corresponding machine.
"""
function oapply(d::WiringDiagram, ms::AbstractDict)
    oapply(d, [ms[box.value] for box in boxes(d)])
end

### Helper functions for `oapply`

function induced_dynamics(d::WiringDiagram, ms::Vector{M}, S, Inputs) where {T,I, M<:AbstractMachine{T,I}}

    function v(u::AbstractVector, xs::AbstractVector, p, t::Real)
        states = destruct(S, u) # a list of the states by box
        readouts = get_readouts(ms, states, p, t)
        readins = unit(I, length(apex(Inputs)))
        fill_readins!(readins, d, Inputs, readouts, xs)

        reduce(vcat, map(enumerate(destruct(Inputs, readins))) do (i,x)
            eval_dynamics(ms[i], states[i], x, p, t)
        end)
    end

end

function induced_dynamics(d::WiringDiagram, ms::Vector{M}, S, Inputs) where {T,I, M<:DelayMachine{T,I}}

    function v(u::AbstractVector, xs::AbstractVector, h::Function, p, t::Real)
        states = destruct(S, u) # a list of the states by box
        hists = destruct(S, h)
        readouts = get_readouts(ms, states, hists, p, t)
        readins = unit(I, length(apex(Inputs)))
        fill_readins!(readins, d, Inputs, readouts, xs)

        reduce(vcat, map(enumerate(destruct(Inputs, readins))) do (i,x)
            eval_dynamics(ms[i], states[i], x, hists[i], p, t)
        end)
    end
end

function induced_readout(d::WiringDiagram, ms::Vector{M}, S) where {T, I, M<:AbstractMachine{T,I}}
    function r(u::AbstractVector, p, t)
        states = destruct(S, u)
        readouts = get_readouts(ms, states, p, t)
        outputs = unit(I, length(output_ports(d)))
        fill_outputs!(outputs, d, readouts)
    end
end

function induced_readout(d::WiringDiagram, ms::Vector{M}, S) where {T, I, M<:DelayMachine{T,I}}
    function r(u::AbstractVector, h, p, t)
        states = destruct(S, u)
        hists = destruct(S, h)
        readouts = get_readouts(ms, states, hists, p, t)
        outputs = unit(I, length(output_ports(d)))
        fill_outputs!(outputs, d, readouts)
    end
end

"""    fills(m::AbstractMachine, d::WiringDiagram, b::Int)

Checks if `m` is of the correct signature to fill box `b` of the  wiring diagram `d`.
"""
function fills(m::AbstractMachine, d::WiringDiagram, b::Int)
    b <= nboxes(d) || error("Trying to fill box $b, when $d has fewer than $b boxes")
    b = box_ids(d)[b]
    return ninputs(m) == length(input_ports(d,b)) && noutputs(m) == length(output_ports(d,b))
end

destruct(C::Colimit, xs::FinDomFunction) = map(1:length(C)) do i
    collect(compose(legs(C)[i], xs))
end
destruct(C::Colimit, xs::AbstractVector) = destruct(C, FinDomFunction(xs))

destruct(C::Colimit, h::Function) = map(1:length(C)) do i
    (p,t) -> destruct(C, h(p,t))[i]
end

get_readouts(ms::AbstractArray{M}, states, args...) where {M <: AbstractMachine} = map(enumerate(ms)) do (i, m)
    readout(m, states[i], args...)
end

get_readouts(ms::AbstractArray{M}, states, hists, p, t) where {M<:DelayMachine} = map(enumerate(ms)) do (i, m)
    readout(m, states[i], hists[i], p, t)
end

function fill_readins!(readins, d::WiringDiagram, Inputs::Colimit, readouts, xs)

    for w in wires(d, :Wire)
        readins[legs(Inputs)[w.target.box](w.target.port)] += readouts[w.source.box][w.source.port]
    end
    for w in wires(d, :InWire)
        readins[legs(Inputs)[w.target.box](w.target.port)] += xs[w.source.port]
    end

    return readins
end

function fill_outputs!(outs, d::WiringDiagram, readouts)

    for w in wires(d, :OutWire)
        outs[w.target.port] += readouts[w.source.box][w.source.port]
    end

    return outs
end

end #module


text:
type SAGAsd <: StepData
	A
	y
	m::Int64
	numChunks::Int64
	batchSize::Int64
	getStep::Function
	stepString::AbstractString
	function SAGAsd(numVars::Int64,y, stepString::AbstractString, numDp::Int64,numChunks::Int64, batchSize::Int64)
		A = sum(y)
		new(A, y, numDp, numChunks,batchSize, SAGAComputation, stepString)
	end
end

function SAGA(numVars::Int64, y, numDP::Int64,numChunks::Int64,batchSize::Int64)
	SAGAsd(numVars,y, "SAGA.$numChunks.$batchSize",  numDP,numChunks,batchSize)
end

function SAGAComputation(x, k, gnum, sd::SAGAsd, problem::Problem; outputLevel =0 )

	outputLevel >0 && println("SAGAComputation started")

	B=zeros(size(x))

	sumy=zeros(size(x))

	indices = sample(1:sd.numChunks,sd.batchSize;replace=false)
	nummm=0
	for i in indices

		B +=sd.y[i]

		sd.y[i]=zeros(size(x))

		for j in (i-1)*round(Int,floor(sd.m/sd.numChunks))+1:i*round(Int,floor(sd.m/sd.numChunks))
			outputLevel >1 && println(" Inner loop, j = $j")
			(f,sampleG) = (problem.getSampleFunctionAt(j))(x)
			sd.y[i]+=sampleG
			sumy += sampleG
			nummm +=1
			gnum += 1
		end

		#step 3
	end
	outputLevel >0 && println("Acutal step determination")

	step =((1/sd.m)*sd.A) - (1/nummm)*(B - sumy )

	outputLevel >0 && println("Updating current sum")
	sd.A=sd.A-B+sumy
	outputLevel >0 && println("SAGAComputation finished")
	(step,gnum)

end


text:
export mkchans

function mkchans(chans::Union{Integer, UnitRange, Array{Int64,1}}, n::Int64)

  if chans == Int64[]
    return Int64.(collect(1:n))
  elseif typeof(chans) <: UnitRange
    return Int64.(collect(chans))
  elseif typeof(chans) <: Integer
    return [Int64(chans)]
  else
    return chans
  end
end


text:
@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Beta},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPBetaOutNPP)

@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Function},
                :inbound_types => (Message, Nothing, Message),
                :name          => SPBetaAMNM)

@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Function},
                :inbound_types => (Message, Message, Nothing),
                :name          => SPBetaBMMN)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Beta},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBBetaOut)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBBetaA)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VBBetaB)

mutable struct SPBetaOutNMM <: SumProductRule{Beta} end
outboundType(::Type{SPBetaOutNMM}) = Message{SampleList}
function isApplicable(::Type{SPBetaOutNMM}, input_types::Vector{Type})
  nothing_inputs = 0
  point_inputs = 0

  for input_type in input_types
      if input_type == Nothing
          nothing_inputs += 1
      elseif input_type << Message{PointMass}
          point_inputs += 1
      end
  end

  return (nothing_inputs == 1) && (point_inputs != 2)
end


text:
if abspath(PROGRAM_FILE) == @__FILE__
  if length(ARGS) < 1
    return
  end
  open(Base.ARGS[1]) do f
    lastnum = 0
    result = 0
    for line in eachline(f)
      n = parse(Int, line)
      if n > lastnum
        result += 1
      end
      lastnum = n
    end
    println(result - 1)
  end
end


text:
# Univariate Normal setup, following Richardson & Green (1997).
# This version uses sampler code that is optimized for univariate normal.
# Note: this has been modified to allow for a fixed prior
module Normal

module NormalModel # submodule for component family definitions
export Theta, Data, likelihood, log_likelihood, prior_sample, prior_sample!, log_prior,
       construct_hyperparameters, update_parameter!, update_hyperparameters!

include("Random.jl")
using .Random

const Theta = Array{Float64,1}  # theta = [mu, sigma]
const Data = Float64

# Normal distribution
const constant = 0.5*log(2.0*pi)
normpdf(x) = exp(-0.5*x*x - constant)
normpdf(x,mu,sigma) = normpdf((x-mu)/sigma)/sigma
log_normpdf(x,mu,sigma) = (r=(x-mu)/sigma;  -0.5*r*r - log(sigma) - constant)
# Gamma distribution
gammapdf(x,a,b) = x^(a-1)*exp(-b*x)*(b^a)/gamma(a)
log_gammapdf(x,a,b) = (a-1)*log(x) - b*x + a*log(b) - lgamma(a)

# Likelihood: Normal(x|mu,sigma^2)
likelihood(x,theta) = normpdf(x,theta[1],theta[2])
likelihood(x,theta,c) = normpdf(x,theta[1,c],theta[2,c])
log_likelihood(x,theta) = log_normpdf(x,theta[1],theta[2])

type Hyperparameters
    m::Float64  # prior mean of mu
    s::Float64  # prior stddev of mu
    a::Float64  # prior shape of sigma
    b::Float64  # prior rate of sigma (starting value)
    g::Float64  # (hyper)prior shape of b
    h::Float64  # (hyper)prior rate of b
    alpha_wt::Float64 # alpha weight (not technically a hyperparameter)
end

function construct_hyperparameters(options)
    x = options.x
    # this is the full data set
    x_full = options.x_full
    # Use values in Green & Richardson (1997) to enable comparison.
    println("Using full data set for hyperparameters")
    m = (minimum(x_full) + maximum(x_full))/2
    R = maximum(x_full) - minimum(x_full)
    s = R
    a = 2.0
    b = 1.0  # b will be updated in Gibbs sampling
    g = 0.2
    h = 10/R^2

    # specify the weight
    alpha_eval = eval(parse(options.alpha_wt_str))
    alpha_wt_fn(n) = Base.invokelatest(alpha_eval, n)
    alpha_wt = alpha_wt_fn(options.n)

    return Hyperparameters(m,s,a,b,g,h, alpha_wt)
end

# Prior (Base distribution)
sample_mu(H) = H.s*randn()+H.m
sample_sigma(H) = sqrt(Random.inverse_gamma(H.a,H.b))
prior_sample(H) = [sample_mu(H), sample_sigma(H)]
prior_sample!(theta,n,H) = (for i=1:n; theta[1,i] = sample_mu(H); theta[2,i] = sample_sigma(H); end)
prior_sample!(theta,H) = (theta[1] = sample_mu(H); theta[2] = sample_sigma(H))
log_prior(theta,m,s,a,b) = log_normpdf(theta[1],m,s) + log_gammapdf(1/theta[2]^2,a,b) + log(2/theta[2]^3)
log_prior(theta,H) = log_prior(theta,H.m,H.s,H.a,H.b)

# Update parameter for each component
function update_parameter!(theta_a, theta_b, x, z, c, H, active)
    n = 0
    sum_xc = 0.0
    for i = 1:length(z)
        if z[i] == c
            sum_xc += x[i]
            n += 1
        end
    end

    # update mean
    ell = 1 / H.s^2
    lambda = 1 / (theta_a[2]^2)
    L = lambda * n * H.alpha_wt + ell ###
    M = (lambda * sum_xc * H.alpha_wt + ell * H.m) / L ###
    if active
        theta_b[1] = randn() / sqrt(L) + M
    end

    # update standard deviation
    alpha = H.a + 0.5*n * H.alpha_wt ###
    r = 0.0
    for i = 1:length(z)
        if z[i]==c
            r += (x[i]-theta_b[1])*(x[i]-theta_b[1])
        end
    end
    beta = H.b + 0.5 * r * H.alpha_wt ###
    # if active; theta_b[2] = 1/sqrt(rand(Gamma(alpha,1/beta))); end
    if active
        theta_b[2] = sqrt(Random.inverse_gamma(alpha,beta))
    end

    return log_prior(theta_b,M,1/sqrt(L),alpha,beta)
end

# Update hyperparameters
function update_hyperparameters!(H,theta,ca,t)
    lambda = 1.0./(theta[2,ca[1:t]].^2)
    alpha = H.g + H.a*t
    beta = H.h + sum(lambda)
    #H.b = rand(Gamma(alpha, 1/beta))
    H.b = Random.gamma(alpha, beta)
end

end # module NormalModel
using .NormalModel

# Include generic code
include("generic.jl")

# Include core sampler code
include("coreNormal.jl")

end # module Normal


text:
import TrajectoryOptimization: get_dt_traj

# Set up problem
model = Dynamics.car
n,m = model.n, model.m
Q = (1e-2)*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)
x0 = [0.;0.;0.]
xf = [0.;1.;0.]

xf = [0,1.0,0]
goal_con = goal_constraint(xf)
circle_con = TO.planar_obstacle_constraint(model.n, model.m, (0,2.5), 0.25)
bnd = BoundConstraint(model.n, model.m, x_min=[-0.5,-0.001,-Inf], x_max=[0.5, 1.001, Inf], u_min=-2, u_max=2)

# Initial Controls
N = 101
U0 = [ones(model.m) for k = 1:N]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)

# Create Problem
n,m = model.n, model.m
prob = Problem(rk3(model), obj, constraints=Constraints([bnd,circle_con],N), N=N, tf=3.0)
prob.constraints[N] += goal_con
initial_controls!(prob, U0)
rollout!(prob)
prob = TO.update_problem(prob, model=model)

# Create DIRCOL Solver
opts = TO.DIRCOLSolverOptions{Float64}(verbose=false)
pcon = prob.constraints
dircol = TO.DIRCOLSolver(prob, opts)

N = prob.N
@test length(dircol.Z.U) == N
part_z = create_partition(n,m,N,N)
NN = N*(n+m)
Z = rand(NN)

# Convert Z to X,U
Z0 = Primals(Z,n,m)
Z0 = Primals(Z,part_z)
X = Z0.X
U = Z0.U
@test size(Z0) == (n,m,N)
@test length(Z0) == NN

Primals(Z,n,m).X[1][1] = 10
@test X[1][1] == 10
@test Z[1] == 10

Z1 = copy(Z)
Z1[1] = 100
Z1 = Primals(Z1,X,U)
@test Z1.Z[1] == 100
@test Z1.X[1][1] == 100
@test X[1][1] == 100

Z2 = Primals(Z, Z0)
@test Z2.Z == Z

# Convert fom X,U to Z
X = [rand(n) for k = 1:N]
U = [rand(m) for k = 1:N]
Z2 = Primals(X,U)
@test Z2.equal == true
@test Primals(prob, true).equal == true
Z_rollout = Primals(prob, true)

Z3 = TO.pack(prob)
@test length(Z3) == NN

@test TO.pack(X,U, part_z) == Z2.Z
@test TO.pack(X,U) == Z2.Z

# Test methods

# Cost
Z = Primals(prob, true)
X,U = Z.X, Z.U
@test Z.X isa TrajectoryOptimization.AbstractVectorTrajectory
@test typeof(Z.X) <: Vector{S} where S <: AbstractVector
initial_controls!(prob, Z.U)
TO.initial_states!(prob, Z.X)
@test cost(prob, dircol) == cost(prob)

# Collocation constraints
p_colloc = n*(N-1)
p_custom = sum(num_constraints(prob))
g_colloc = zeros(p_colloc)
g = zeros(p_colloc + p_custom)
# TO.dynamics!(prob, dircol, X, U)
# traj_points!(prob, dircol, Z_rollout)
TO.collocation_constraints!(g_colloc, prob, dircol, X, U)
@test norm(g_colloc,Inf) < 1e-15

# Normal Constraints
g_custom = view(g, p_colloc.+(1:p_custom))
TO.update_constraints!(g, prob, dircol, X, U)
p = num_constraints(prob)
@test p == dircol.p

# Cost gradient
# grad_f = zeros(NN)
# TO.cost_gradient!(grad_f, prob, X, U, get_dt_traj(prob))

# Z = Primals(Z, part_z)
#
# cost(prob.obj, Z.X, Z.U, get_dt_traj(prob))
#
function eval_f(Z)
    Z = Primals(Z, part_z)
    cost(prob.obj, Z.X, Z.U, get_dt_traj(prob))
end
eval_f(Z.Z)
# @test ForwardDiff.gradient(eval_f, Z.Z) == grad_f

# Collocation Constraint jacobian
nG_colloc = p_colloc*2(n+m)
dircol = DIRCOLSolver(prob, opts, Z)

jac = zeros(p_colloc, NN)
TO.collocation_constraint_jacobian!(jac, prob, dircol, X, U)
jac

function jac_colloc(Z)
    X,U = TO.unpack(Z, part_z)
    rows,cols = zeros(nG_colloc), zeros(nG_colloc)
    jac_struct = spzeros(p_colloc,NN)
    TO.collocation_constraint_jacobian_sparsity!(jac_struct, prob)
    r,c = TO.get_rc(jac_struct)
    jac = zeros(nG_colloc)
    TO.collocation_constraint_jacobian!(jac, prob, dircol, X, U)
    return sparse(r,c,jac)
end
jac_colloc(Z.Z)

function c_colloc(Z)
    X,U = TO.unpack(Z, part_z)
    g_colloc = zeros(eltype(Z), p_colloc)
    Z = Primals(Z, part_z)
    solver = DIRCOLSolver(prob, DIRCOLSolverOptions{eltype(Z.Z)}(), Z)
    TO.collocation_constraints!(g_colloc, prob, solver, X, U)
    return g_colloc
end
c_colloc(Z.Z)

@test ForwardDiff.jacobian(c_colloc, Z.Z) ≈ Array(jac_colloc(Z.Z))

# General constraint jacobian
nG_custom = sum(p[1:N-1])*(n+m) + p[N]*n
nG = nG_colloc + nG_custom
jac = spzeros(p_colloc+p_custom, NN)
TO.constraint_jacobian_sparsity!(jac, prob)
rows,cols = TO.get_rc(jac)

function jac_con(Z)
    X,U = TO.unpack(Z, part_z)

    jac = zeros(nG)
    nG_colloc = p_colloc * 2(n+m)
    jac_colloc = view(jac, 1:nG_colloc)
    solver = DIRCOLSolver(prob, DIRCOLSolverOptions{eltype(Z)}(), Primals(Z, part_z))
    TO.collocation_constraint_jacobian!(jac_colloc, prob, solver, X, U)

    # General constraint jacobians
    jac_custom = view(jac, nG_colloc+1:length(jac))
    TO.constraint_jacobian!(jac_custom, prob, solver, X, U)

    return sparse(rows,cols,jac)
end
jac_con(Z.Z)

function c_all(Z)
    X,U = TO.unpack(Z, part_z)
    g = zeros(eltype(Z), p_colloc+p_custom)
    g_colloc = view(g,1:p_colloc)
    g_custom = view(g,p_colloc+1:length(g))
    solver = DIRCOLSolver(prob, DIRCOLSolverOptions{eltype(Z)}(), Primals(Z, part_z))
    TO.collocation_constraints!(g_colloc, prob, solver, X, U)
    TO.update_constraints!(g_custom, prob, solver, X, U)
    return g
end
g = c_all(Z.Z)
@test ForwardDiff.jacobian(c_all, Z.Z) ≈ jac_con(Z.Z)

# Test Ipopt functions
eval_f2, eval_g, eval_grad_f, eval_jac_g = TO.gen_dircol_functions(prob, dircol)

g2 = zero(g)
# grad_f2 = zero(grad_f)
jac2 = zeros(nG)
row, col = zeros(nG), zeros(nG)
@test eval_f2(Z.Z) == eval_f(Z.Z)
eval_g(Z.Z, g2)
@test g ≈ g2
# eval_grad_f(Z.Z, grad_f2)
# @test grad_f ≈ grad_f2

eval_jac_g(Z.Z, :Structure, row, col, jac2)
eval_jac_g(Z.Z, :Values, row, col, jac2)
@test sparse(row,col,jac2) == jac_con(Z.Z)

# Strip out state and control Bounds
prob0 = copy(prob)
p0 = num_constraints(prob0)
bnds = TO.remove_bounds!(prob0)
p = num_constraints(prob0)
@test p0[1] == 9
@test p[1] == 1
@test p[N] == 1
z_U,z_L,g_U,g_L = TO.get_bounds(prob0,bnds)
Z_U = Primals(z_U, part_z)
Z_L = Primals(z_L, part_z)
@test Z_U.X[1] == Z_L.X[1] == prob0.x0
@test Z_U.X[N] == Z_L.X[N] == xf
@test Z_U.U[N] == Z_U.U[N-1]

# @test_nowarn solve!(prob, opts)
# @test_nowarn TO.write_ipopt_options()

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using Libdl

depsfile = joinpath(dirname(@__FILE__),"deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function my_download(url::AbstractString, filename::AbstractString)
    if Sys.iswindows() && VERSION < v"1.1"
        # from https://github.com/JuliaLang/julia/blob/788b2c77c10c2160f4794a4d4b6b81a95a90940c/base/download.jl#L26
        # otherwise files are not properly downloaded - check sizes
        curl_exe = joinpath(get(ENV, "SYSTEMROOT", "C:\\Windows"), "System32\\curl.exe")
        process = run(`$curl_exe -s -S -g -L -f -o $filename $url`)#, wait=false)
        return filename
    else
        return download(url, filename)
    end
end

function write_depsfile(path)
    f = open(depsfile,"w")
    if Sys.iswindows()
        path = replace(path, "\\" => "\\\\")
    end
    print(f,"""
    const xpressdlpath = \"$(path)\"
    """)
    close(f)
end

function local_installation()

    libname = string(Sys.iswindows() ? "" : "lib", "xprs", ".", Libdl.dlext)
    paths_to_try = String[]

    push!(paths_to_try, "")
    push!(paths_to_try, @__DIR__)

    if haskey(ENV, "XPRESSDIR")
    push!(paths_to_try, joinpath(ENV["XPRESSDIR"], Sys.iswindows() ? "bin" : "lib"))
    end

    global found = false
    for l in paths_to_try
        path = joinpath(l, libname)
        d = Libdl.dlopen_e(path)
        if d != C_NULL
            global found = true
            @info("Found $path")
            write_depsfile(l)
            break
        end
    end

    if !found
        error("""
        Unable to locate Xpress installation.
        Please check your enviroment variable XPRESSDIR.
        Note that Xpress must be obtained separately from fico.com.
        """)
    end
end

function ci_installation()
    files = if Sys.iswindows()
    [
        (ENV["SECRET_XPRS_WIN_8110"], "xprs.dll")
        (ENV["SECRET_XPRL_WIN_8110"], "xprl.dll")
        (ENV["SECRET_XPRA_WIN_8110"], "xpauth.xpr")
    ]
    end
    for (url, file) in files
        local_filename = joinpath(@__DIR__, file)
        my_download(url, local_filename)
    end
    path = joinpath(@__DIR__, files[1][2])
    d = Libdl.dlopen_e(path)
    if d != C_NULL
        write_depsfile(@__DIR__)
    else
        error("Could not open xprs.dll")
    end
end

if haskey(ENV, "XPRESS_JL_SKIP_LIB_CHECK")
    # Skip!
elseif get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
    write_depsfile("julia_registryci_automerge")
elseif get(ENV, "SECRET_XPRS_WIN_8110", "") != ""
    ci_installation()
else
    local_installation()
end


text:
tests = ["texpr",
         "tmode",
         "scalar_ref",
         "scalar_ewise",
         "scalar_special",
         "scalar_reduc",
         "scalar_hybrid",
         "extensions" ]

# extensions

println("Testing Devectorize.jl")
println("---------------------------")

for t in tests
    fp = "$(t).jl"
    println("* running $(fp) ...")
    include(fp)
end


text:
using TypeDB_tutorial
using Documenter

DocMeta.setdocmeta!(TypeDB_tutorial, :DocTestSetup, :(using TypeDB_tutorial); recursive=true)

makedocs(;
    modules=[TypeDB_tutorial],
    authors="Frank Urbach",
    repo="https://GitHub.com/FrankUrbach/TypeDB_tutorial.jl/blob/{commit}{path}#{line}",
    sitename="TypeDB_tutorial.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://FrankUrbach.github.io/TypeDB_tutorial.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="GitHub.com/FrankUrbach/TypeDB_tutorial.jl",
)


text:
infile = "./src/Day12/input.txt"
input = readlines(infile)

import Base.:+

import Base.:*

directions = [(x[1:1], parse(Int,x[2:end]) ) for x in input]

struct Point
    x
    y
end
+(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y)
*(a::Point, n::Int) = Point(n*a.x, n*a.y)

function rotate(facing, dir, num)
    x₁ , y₁ = facing.x, facing.y
    if dir == "R"
        num = -num
    end
    x₂ = x₁*cosd(num) - y₁*sind(num)
    y₂ = x₁*sind(num) + y₁*cosd(num)
    return Point(x₂, y₂)
end

testdir = ["F10", "N3", "F7", "R90", "F11"]
testdirections = [(x[1:1], parse(Int,x[2:end]) ) for x in testdir]

function follow_orders(directions, pos, facing)
    visited = Array{Tuple{Point, Point},1}()
    push!(visited,(pos, facing))
    for (d, num) in directions
        if d == "N"
            pos = pos + Point(0, num)
        elseif d == "S"
            pos = pos + Point(0, -num)
        elseif d == "E"
            pos = pos + Point(num, 0)
        elseif d == "W"
            pos = pos + Point(-num, 0)
        elseif d == "L" || d == "R"
            # rotate left by num degrees
            facing = rotate(facing, d, num)
        elseif d == "F"
            pos = pos + (facing*num)
        end
        push!(visited, (pos, facing))
    end
    return pos, visited
end

final_pos, history = follow_orders(directions, Point(0,0), Point(-1, 0))

final_pos
abs(final_pos.x) + abs(final_pos.y)

function track_waypoint(directions, pos, waypoint)
    visited = Array{Tuple{Point, Point},1}()
    push!(visited,(pos, waypoint))
    for (d, num) in directions
        if d == "N"
            waypoint = waypoint + Point(0, num)
        elseif d == "S"
            waypoint = waypoint + Point(0, -num)
        elseif d == "E"
            waypoint = waypoint + Point(num, 0)
        elseif d == "W"
            waypoint = waypoint + Point(-num, 0)
        elseif d == "L" || d == "R"
            # rotate left by num degrees
            waypoint = rotate(waypoint, d, num)
        elseif d == "F"
            pos = pos + (waypoint*num)
        end
        push!(visited, (pos, waypoint))
    end
    return pos, visited
end

p, hist = track_waypoint(directions, Point(0,0), Point(10, 1))
abs(p.x) + abs(p.y)