text:
import MapImages

"""Threaded version: If we can load all stack file into memory, might be much quicker"""
function run_sbas(
    unw_stack::AbstractArray{<:AbstractFloat},
    outfile::String,
    outdset::String,
    geolist,
    intlist,
    constant_velocity::Bool,
    alpha::Real,
    L1::Bool = false,
    prune_outliers = true,
    sigma = 3,
    prune_fast = false;
    demrsc = nothing,
    reference_station = nothing,
    ref_row = nothing,
    ref_col = nothing,
    window = 5,
)

    # Check if this file/dset already exists
    isfile(outfile) &&
    outdset in names(outfile) &&
    error("$outdset already exists in $outfile")

    # TODO : extract from stackavg
    if !isnothing(reference_station)
        println("Using $reference_station as reference")
        ref_row, ref_col = MapImages.station_rowcol(reference_station, demrsc)
    end
    if !(isnothing(ref_row) && isnothing(ref_col))
        println("Shifting input stack to $ref_row, $ref_col")
        # Using col, row since we
        @time unw_stack .=
            InsarTimeseries.shift_stack(unw_stack, ref_row, ref_col, window = window)
    end
    L1 ? println("Using L1 penalty for fitting") :
    println("Using least squares for fitting")
    nrows, ncols, _ = size(unw_stack)

    outstack = Array{Float32,2}(undef, (nrows, ncols))
    println("Out size: ($nrows, $ncols)")

    @time Threads.@threads for col = 1:ncols
        for row = 1:nrows
            # soln_phase, igram_count, geo_clean, unw_clean
            unw_pixel_raw = unw_stack[row, col, :]
            soln_phase, _, _, _ = InsarTimeseries.calc_soln(
                unw_pixel_raw,
                geolist,
                intlist,
                alpha,
                constant_velocity;
                L1 = L1,
                prune_outliers = prune_outliers,
                sigma = sigma,
                prune_fast = prune_fast,
            )
            outstack[row, col] = InsarTimeseries.P2MM * soln_phase[1]
        end
    end

    println("Writing solution into $outfile : $outdset")
    h5write(outfile, outdset, permutedims(outstack))

    !isnothing(demrsc) && Sario.save_dem_to_h5(outfile, demrsc, overwrite = true)
    Sario.save_geolist_to_h5(outfile, outdset, geolist, overwrite = true)

    return outfile, outdset
end

"""Solve the problem Bv = d for each pixel in the stack"""
function invert_sbas(
    unw_stack::Union{HDF5Dataset,Array{Float32,3}},
    B::Array{Float32,2},
    valid_igram_indices;
    extra_zeros = 0,
    L1 = true,
)
    nrows, ncols, nlayers = size(unw_stack)
    total_pixels = nrows * ncols * nlayers

    num_igrams, num_timediffs = size(B)

    # Output holder
    vstack = Array{Float32,3}(undef, (nrows, ncols, num_timediffs))

    println("Using $(Threads.nthreads()) threads for invert_sbas loop")
    # @inbounds Threads.@threads for j in 1:ncols
    # @inbounds for i in 1:nrows
    step = 1000
    row = 1
    col = 1

    # Speeds up inversion to precompute pseudo inverse for L2 least squares case
    if L1
        # v = Convex.Variable(num_timediffs)
        lu_tuple = factor(Float64.(B))
        # TODO: seems to be pretty high variance for alpha...
        # figure out which params are best/ how to adjust on fly
        # Settings found to be okay for now
        rho, alpha, abstol = 1.0, 1.6, 1e-3
    else
        pB = pinv(B)
    end

    chunk = zeros(Float32, (step, step, length(valid_igram_indices)))
    pix_count = 0
    # clear_mem_every = 10000  # TODO: figure out what causes Convex.jl to grow memory and slow down
    while col <= ncols
        while row <= nrows
            end_row = min(row + step - 1, nrows)
            end_col = min(col + step - 1, ncols)

            # If the chunk is not a full square, make sure we assign
            # only part to the chunk buffer to match broadcast
            col_c = length(col:end_col)
            row_c = length(row:end_row)

            # TODO: check that reading all depth, then subselecting is actually better than
            # For-looping and reading one layer at a time
            println("Reading new $row_c x $col_c chunk: ($row:$end_row, $col:$end_col)")
            @time chunk[1:row_c, 1:col_c, :] .=
                unw_stack[row:end_row, col:end_col, :][:, :, valid_igram_indices]

            last_time = time()
            @inbounds Threads.@threads for j = 1:col_c
                @inbounds for i = 1:row_c
                    if L1
                        # vstack[row+i-1, col+j-1, :] .= invert_pixel(chunk[i, j, :], B, v)
                        vstack[row + i - 1, col + j - 1, :] .= invert_pixel(
                            chunk[i, j, :],
                            B,
                            rho = rho,
                            alpha = alpha,
                            lu_tuple = lu_tuple,
                            abstol = abstol,
                        )
                    else
                        vstack[row + i - 1, col + j - 1, :] .=
                            invert_pixel(chunk[i, j, :], pB, extra_zeros)
                    end
                    pix_count += 1
                    last_time = log_count(
                        pix_count,
                        total_pixels,
                        nlayers,
                        every = 10000,
                        last_time = last_time,
                    )
                    # Because of a memory bug in either ECOS or Convex.jl
                    # if pix_count % clear_mem_every == 0
                    # println("CLEARING MEMORY")
                    # Convex.clearmemory()
                    # GC.gc()
                    # v = Convex.Variable(size(B, 2))
                    # end
                end
            end
            row += step
        end

        row = 1
        col += step
    end

    return vstack
end

# function invert_pixel(pixel::Array{Float32, 1}, B::Array{Float32,2}, v::Convex.Variable)
#     # Note: these are defined here to allow multithreading and not reuse across threads
#     solver = ECOS.ECOSSolver(verbose=0)
#     problem = Convex.minimize(norm(B*v - pixel, 1))
#     Convex.solve!(problem, solver)
#
#     if length(v.value) > 1
#         Float32.(vec(v.value))
#     else
#         Float32.([v.value])
#     end
# end

function log_count(pix_count, total_pixels, nlayers; every = 100_000, last_time = nothing)
    if (pix_count % every) == 0
        flush(stdout)
        pct_done = 100 * pix_count * nlayers / total_pixels
        @printf(
            "Processed %.3g pixels out of %.3g (%.2f%% done).",
            pix_count * nlayers,
            total_pixels,
            pct_done
        )
        if !isnothing(last_time)
            t = time()
            elapsed = t - last_time
            @printf(" %.2f seconds elapsed\n", elapsed)
            return t
        else
            @printf("\n")
        end
    else
        return last_time
    end
end

# function save_line_fit(unreg_fname)
#     geolist = InsarTimeseries.load_geolist_from_h5(unreg_fname)
#     geolist_nums = [( g - geolist[1]).value for g in geolist]
#     f = h5open(unreg_fname)
#     dset = f["stack"]
#     fname_out = replace(unreg_fname, ".h5" => "_linefit.h5")
#     fout = h5open(fname_out, "w")
#     d_create(fout, "stack", datatype(Float32), dataspace( (size(dset, 1), size(dset, 2)) ))
#     dset_out = fout["stack"]
#     for i in 1:size(dset, 1)
#         for j in 1:size(dset, 2)
#             p = polyfit(geolist_nums, vec(dset[i, j, :]), 1)
#             dset_out[i, j] = p(geolist_nums[end])
#         end
#     end
#     close(fout); close(f)
# end
#
# pix_count += 1
# last_time = log_count(pix_count, total_pixels, 1, every=10000, last_time=last_time)

function integrate_velocities(vstack::Array{<:AbstractFloat,3}, timediffs::Array{Int,1})
    nrows, ncols, _ = size(vstack)
    num_geos = length(timediffs) + 1
    phi_stack = zeros(nrows, ncols, num_geos)
    phi_diffs = zeros(num_geos - 1)

    phi_arr = zeros(num_geos)  # Buffer to hold each result
    for j = 1:ncols
        for i = 1:nrows
            varr = view(vstack, i, j, :)
            phi_stack[i, j, :] .= integrate1D!(phi_diffs, phi_arr, varr, timediffs)
        end
    end
    return phi_stack
end

function integrate1D!(
    phi_diffs,
    phi_arr,
    velocities::AbstractArray{<:AbstractFloat,1},
    timediffs::Array{Int,1},
)
    # multiply each column of vel array: each col is a separate solution
    phi_diffs .= velocities .* timediffs

    # Now the final phase results are the cumulative sum of delta phis
    # This is equivalent to replacing missing with 0s (like np.ma.cumsum does)
    phi_arr[2:end] .= cumsum(coalesce.(phi_diffs, 0))
    return phi_arr
end

# """2nd Threaded version: If we can load all stack file into memory, might be much quicker"""
# function run_sbas(unw_stack::AbstractArray{<:AbstractFloat},
#                   outfile::String,
#                   outdset::String,
#                   geolist,
#                   intlist,
#                   constant_velocity::Bool,
#                   alpha::Float32,
#                   L1::Bool=false,
#                   prune::Bool=true)
#
#     B = prepB(geolist, intlist, constant_velocity, alpha)
#     L1 ? println("Using L1 penalty for fitting") : println("Using least squares for fitting")
#     # TODO: do I ever really care about the abstol to change as variable?
#     rho, alpha, abstol = 1.0, 1.6, 1e-3
#     nrows, ncols, _ = size(unw_stack)
#
#     outstack = Array{Float32, 2}(undef, (nrows, ncols))
#     countstack = similar(outstack)
#     println("Out size: ($nrows, $ncols)")
#     pix_count, total_pixels = 0, nrows * ncols
#
#     @time Threads.@threads for col in 1:ncols
#         for row in 1:nrows
#             soln_p2mm, igram_count = calc_soln(view(unw_stack, row, col, :), B, geolist, intlist, rho, alpha, L1, prune)
#             outstack[row, col] = soln_p2mm
#             countstack[row, col] = igram_count
#         end
#     end
#
#     println("Writing solution into $outfile : $outdset")
#     h5open(outfile, "cw") do f
#         f[outdset] = outstack
#         f[_count_dset(outdset)] = countstack
#     end
#     return outfile, outdset
# end


text:
export GSSmoother,smooth
using SparseArrays
using LinearAlgebra
using BenchmarkTools

struct GSSmoother
    invh2::Float64
    bc::NTuple{2,NTuple{2,BoundaryCondition}}
    finest::Bool
    function GSSmoother(invh2,nrows,ncols,bc,finest)
        new(invh2,bc,finest)
    end
end

@inline function sij(denom,i,j,sl,rl,ih2,nx,ny)
    @inbounds denom*(rl[i,j]+ih2*(sl[i,j-1]+sl[i-1,j]+sl[i+1,j]+sl[i,j+1]))
end

using LoopVectorization

@inline function smooth_line(nrm1,j,i1,sl,rl,ih2,denom)
    if nrm1 >= 0
    @fastmath @simd for i=i1:2:nrm1
            sl[i,j]=denom*(rl[i,j]+ih2*(sl[i,j-1]+sl[i-1,j]+sl[i+1,j]+sl[i,j+1]))
        end
    else
        # @show nrm1
        @fastmath @inbounds @simd for i=i1:2:nrm1
            sl[i,j]=denom*(rl[i,j]+ih2*(sl[i,j-1]+sl[i-1,j]+sl[i+1,j]+sl[i,j+1]))
        end
    end
end

@inline function smooth_seq(nrows,ncols,sl,rl,ih2)
    denom=1/(4ih2)
    nrm1=nrows-1

    @inbounds for j=2:2:ncols-1
        smooth_line(nrm1,j,2,sl,rl,ih2,denom)
        smooth_line(nrm1,j+1,3,sl,rl,ih2,denom)
    end

    # @inbounds for j=2:2:ncols-1
    #     smooth_line(nrm1,j,2,sl,rl,ih2,denom)
    #     smooth_line(nrm1,j+1,3,sl,rl,ih2,denom)
    # end
    @inbounds for j=2:2:ncols-1
        smooth_line(nrm1,j,3,sl,rl,ih2,denom)
        smooth_line(nrm1,j+1,2,sl,rl,ih2,denom)
    end

end

@inline function smooth_cart(nrows,ncols,sl,rl,ih2)
    denom=1/(4ih2)
    nrm1=nrows-1
    ncm1=ncols-1
    @avx for j=2:ncm1
         for i=2:nrm1
            sl[i,j]=denom*(rl[i,j]+ih2*(sl[i,j-1]+sl[i-1,j]+sl[i+1,j]+sl[i,j+1]))
        end
    end

    # @inbounds for j=2:2:ncols-1
    #     smooth_line(nrm1,j,2,sl,rl,ih2,denom)
    #     smooth_line(nrm1,j+1,3,sl,rl,ih2,denom)
    # end
    # @inbounds for j=2:2:ncols-1
    #     smooth_line(nrm1,j,3,sl,rl,ih2,denom)
    #     smooth_line(nrm1,j+1,2,sl,rl,ih2,denom)
    # end

end

function smooth(sl,rl,a::GSSmoother)
    ih2=a.invh2
    (nrows,ncols)=size(sl)
    @assert(size(rl)==size(sl))
    denom=1/(4ih2)
    # smooth_cart(nrows,ncols,sl,rl,ih2)
    smooth_seq(nrows,ncols,sl,rl,ih2)
    return
end

ext: jl


text:
##############################################
# - Write tests for plot_qq
##############################################

# Import dependencies required
using GLMdiagnostics
using Test
using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Gadfly
using Statistics
using Compose

# Import simulated data and linear model to test functions
include("test_simulated_model.jl")

@testset "plot_qq" begin

    # Test #1
    # - Test that object created by function is of type 'Plot'
    @test   typeof(plot_qq(modAov)) == Plot

end # End of testing for plot_qq function

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_stiffquadito

Random.seed!(100)
prob = prob_sde_stiffquadito
prob = remake(prob;p=(1e5,2.))
alg = AutoSOSRA2(SKenCarp(), maxstiffstep=5, maxnonstiffstep=2, stiffalgfirst=false)
@test StochasticDiffEq.isadaptive(prob,alg)
@test StochasticDiffEq.isadaptive(alg)
@time sol = solve(prob, alg)
@test typeof(alg.algs[sol.alg_choice[end]]) <: SKenCarp
@test length(unique(sol.alg_choice)) == 2

Random.seed!(100)
prob = prob_sde_stiffquadito
prob = remake(prob;p=(1e5,2.))
@time sol = solve(prob, AutoSOSRI2(ImplicitRKMil(), maxstiffstep=1, maxnonstiffstep=10, stiffalgfirst=false))
@test length(unique(sol.alg_choice)) == 2


text:
@testset "Site counting" begin

    dna_alphabets = (DNAAlphabet{4}, DNAAlphabet{2})
    rna_alphabets = (RNAAlphabet{4}, RNAAlphabet{2})

    @testset "Specific count methods" begin

        function generate_possibilities_tester(::Type{A}) where A<:NucAlphs
            symbols = BioSymbols.alphabet(A)
            arra = Vector{eltype(A)}()
            arrb = Vector{eltype(A)}()
            for i in 1:length(symbols), j in i:length(symbols)
                push!(arra, symbols[i])
                push!(arrb, symbols[j])
            end
            return BioSequence{A}(arra), BioSequence{A}(arrb)
        end

        for alphset in (dna_alphabets, rna_alphabets)
            # Answers to these tests were worked out manually to verify
            # counting is working correctly.
            # seqA and seqB contain all possible observations of sites.

            # 4 bit encoded sequences
            seqA, seqB = generate_possibilities_tester(alphset[1])
            @test count(Certain, seqA, seqB) == count(Certain, seqB, seqA) == 10
            @test count(Gap, seqA, seqB) == count(Gap, seqB, seqA) == 16
            @test count(Ambiguous, seqA, seqB) == count(Ambiguous, seqB, seqA) == 121
            @test count(Match, seqA, seqB) == count(Match, seqB, seqA) == length(BioSymbols.alphabet(alphset[1]))
            @test count(Mismatch, seqA, seqB) == count(Mismatch, seqB, seqA) == (length(seqA) - length(BioSymbols.alphabet(alphset[1])))

            # 2 bit encoded sequences
            seqA, seqB = generate_possibilities_tester(alphset[2])
            @test count(Certain, seqA, seqB) == count(Certain, seqB, seqA) == 10
            @test count(Gap, seqA, seqB) == count(Gap, seqB, seqA) == 0
            @test count(Ambiguous, seqA, seqB) == count(Ambiguous, seqB, seqA) == 0
            @test count(Match, seqA, seqB) == count(Match, seqB, seqA) == length(BioSymbols.alphabet(alphset[2]))
            @test count(Mismatch, seqA, seqB) == count(Mismatch, seqB, seqA) == (length(seqA) - length(BioSymbols.alphabet(alphset[2])))
        end
    end

    @testset "Randomized tests" begin

        # A test counting function which is naive.
        @inline function testcount(::Type{P}, a::BioSequence, b::BioSequence) where P<:BioSequences.Position
            k = 0
            @inbounds for idx in 1:min(lastindex(a), lastindex(b))
                k += issite(P, a, b, idx)
            end
            return k
        end
        issite(::Type{Ambiguous}, a::BioSequence, idx) = isambiguous(a[idx])
        @inline function issite(::Type{Ambiguous}, a::BioSequence, b::BioSequence, idx)
            return issite(Ambiguous, a, idx) | issite(Ambiguous, b, idx)
        end
        issite(::Type{Certain}, a::BioSequence, idx) = iscertain(a[idx])
        @inline function issite(::Type{Certain}, a::BioSequence, b::BioSequence, idx)
            return issite(Certain, a, idx) & issite(Certain, b, idx)
        end
        issite(::Type{Gap}, a::BioSequence, idx) = isgap(a[idx])
        @inline function issite(::Type{Gap}, a::BioSequence, b::BioSequence, idx)
            return issite(Gap, a, idx) | issite(Gap, b, idx)
        end
        issite(::Type{Match}, a::BioSequence, b::BioSequence, idx) = a[idx] == b[idx]
        issite(::Type{Mismatch}, a::BioSequence, b::BioSequence, idx) = a[idx] != b[idx]

        # Randomized tests get performed with a naive counting function
        # which is intuitive and works, but that is nowhere near as quick.

        function testcounting(::Type{S}, a, b) where S<:Site
            @test count(S, a, b) == count(S, b, a) == testcount(S, a, b)
        end

        function testforencs(a::Int, b::Int, subset::Bool)
            for alphabet in (DNAAlphabet, RNAAlphabet)
                for _ in  1:50
                    #println("TESTING SEQUENCES:")
                    seqA = random_seq(alphabet{a}, rand(10:100))
                    seqB = random_seq(alphabet{b}, rand(10:100))
                    #=
                    println("A seq: ", seqA)
                    println("B seq: ", seqB)
                    println("A length: ", length(seqA))
                    println("B length: ", length(seqB))
                    =#
                    sa = seqA
                    sb = seqB
                    if subset
                        intA = random_interval(1, length(seqA))
                        intB = random_interval(1, length(seqB))
                        subA = seqA[intA]
                        subB = seqB[intB]
                        #=
                        println("A subset: ", intA)
                        println("B subset: ", intB)
                        println("A subseq: ", subA)
                        println("B subseq: ", subB)
                        =#
                        sa = subA
                        sb = subB
                    end
                    #println("Testing Mismatch")
                    testcounting(Mismatch, sa, sb)
                    #println("Testing Match")
                    testcounting(Match, sa, sb)
                    #println("Testing Certain")
                    testcounting(Certain, sa, sb)
                    #println("Testing Gap")
                    testcounting(Gap, sa, sb)
                    #println("Testing Ambiguous")
                    testcounting(Ambiguous, sa, sb)
                end
            end
        end

        @testset "Testing 4-bit seq pairs" begin
            @testset "Full random sequences" begin
                testforencs(4, 4, false)
            end
            @testset "Subset random sequences" begin
                testforencs(4, 4, true)
            end
        end

        @testset "Testing 2-bit seq pairs" begin
            @testset "Full random sequences" begin
                testforencs(2, 2, false)
            end
            @testset "Subset random sequences" begin
                testforencs(2, 2, true)
            end
        end

        @testset "Testing mixed bit seq pairs" begin
            @testset "Full random sequences" begin
                testforencs(4, 2, false)
                testforencs(2, 4, false)
            end
            @testset "Subset random sequences" begin
                testforencs(4, 2, true)
                testforencs(2, 4, true)
            end
        end
    end

    @testset "Pairwise methods" begin
        @testset "4-bit encoded sequences" begin
            dnas = [dna"ATCGCCA-", dna"ATCGCCTA", dna"ATCGCCT-", dna"GTCGCCTA"]
            rnas = [rna"AUCGCCA-", rna"AUCGCCUA", rna"AUCGCCU-", rna"GUCGCCUA"]
            answer_mismatch = [0 2 1 3;
                               2 0 1 1;
                               1 1 0 2;
                               3 1 2 0]
            answer_match = [0 6 7 5;
                            6 0 7 7;
                            7 7 0 6;
                            5 7 6 0]
            for i in (dnas, rnas)
                @test count_pairwise(Mismatch, i...) == answer_mismatch
                @test count_pairwise(Match, i...) == answer_match
                @test count_pairwise(Certain, i...) == [0 7 7 7;
                                                        7 0 7 8;
                                                        7 7 0 7;
                                                        7 8 7 0]
                @test count_pairwise(Ambiguous, i...) == [0 0 0 0;
                                                          0 0 0 0;
                                                          0 0 0 0;
                                                          0 0 0 0]
                @test count_pairwise(Gap, i...) == [0 1 1 1;
                                                    1 0 1 0;
                                                    1 1 0 1;
                                                    1 0 1 0]
            end
        end
        @testset "2-bit encoded sequences" begin
            dnas = [BioSequence{DNAAlphabet{2}}("ATCGCCAC"),
                    BioSequence{DNAAlphabet{2}}("ATCGCCTA"),
                    BioSequence{DNAAlphabet{2}}("ATCGCCTT"),
                    BioSequence{DNAAlphabet{2}}("GTCGCCTA")]
            rnas = [BioSequence{RNAAlphabet{2}}("AUCGCCAC"),
                    BioSequence{RNAAlphabet{2}}("AUCGCCUA"),
                    BioSequence{RNAAlphabet{2}}("AUCGCCUU"),
                    BioSequence{RNAAlphabet{2}}("GUCGCCUA")]
            answer_mismatch = [0 2 2 3;
                               2 0 1 1;
                               2 1 0 2;
                               3 1 2 0]
            answer_match = [0 6 6 5;
                            6 0 7 7;
                            6 7 0 6;
                            5 7 6 0]
            for i in (dnas, rnas)
                @test count_pairwise(Mismatch, i...) == answer_mismatch
                @test count_pairwise(Match, i...) == answer_match
                @test count_pairwise(Certain, i...) == [0 8 8 8;
                                                        8 0 8 8;
                                                        8 8 0 8;
                                                        8 8 8 0]
                @test count_pairwise(Ambiguous, i...) == [0 0 0 0;
                                                          0 0 0 0;
                                                          0 0 0 0;
                                                          0 0 0 0]
                @test count_pairwise(Gap, i...) == [0 0 0 0;
                                                    0 0 0 0;
                                                    0 0 0 0;
                                                    0 0 0 0]
            end
        end
    end

    @testset "Windowed methods" begin
        @testset "4-bit encoded sequences" begin
            dnaA = dna"ATCGCCA-M"
            dnaB = dna"ATCGCCTAA"
            rnaA = rna"AUCGCCA-M"
            rnaB = rna"AUCGCCUAA"

            for seqs in ((dnaA, dnaB), (rnaA, rnaB))
                @test count(Certain, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                                 IntervalValue(2, 4, 3),
                                                                 IntervalValue(3, 5, 3),
                                                                 IntervalValue(4, 6, 3),
                                                                 IntervalValue(5, 7, 3),
                                                                 IntervalValue(6, 8, 2),
                                                                 IntervalValue(7, 9, 1)]
                @test count(Ambiguous, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                   IntervalValue(2, 4, 0),
                                                                   IntervalValue(3, 5, 0),
                                                                   IntervalValue(4, 6, 0),
                                                                   IntervalValue(5, 7, 0),
                                                                   IntervalValue(6, 8, 0),
                                                                   IntervalValue(7, 9, 1)]
                @test count(Gap, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                             IntervalValue(2, 4, 0),
                                                             IntervalValue(3, 5, 0),
                                                             IntervalValue(4, 6, 0),
                                                             IntervalValue(5, 7, 0),
                                                             IntervalValue(6, 8, 1),
                                                             IntervalValue(7, 9, 1)]
                @test count(Match, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                               IntervalValue(2, 4, 3),
                                                               IntervalValue(3, 5, 3),
                                                               IntervalValue(4, 6, 3),
                                                               IntervalValue(5, 7, 2),
                                                               IntervalValue(6, 8, 1),
                                                               IntervalValue(7, 9, 0)]
                @test count(Mismatch, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                  IntervalValue(2, 4, 0),
                                                                  IntervalValue(3, 5, 0),
                                                                  IntervalValue(4, 6, 0),
                                                                  IntervalValue(5, 7, 1),
                                                                  IntervalValue(6, 8, 2),
                                                                  IntervalValue(7, 9, 3)]
            end
        end
        @testset "2-bit encoded sequences" begin
            dnaA = BioSequence{DNAAlphabet{2}}("ATCGCCATT")
            dnaB = BioSequence{DNAAlphabet{2}}("ATCGCCTAA")
            rnaA = BioSequence{RNAAlphabet{2}}("AUCGCCAUU")
            rnaB = BioSequence{RNAAlphabet{2}}("AUCGCCUAA")

            for seqs in ((dnaA, dnaB), (rnaA, rnaB))
                @test count(Certain, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                                 IntervalValue(2, 4, 3),
                                                                 IntervalValue(3, 5, 3),
                                                                 IntervalValue(4, 6, 3),
                                                                 IntervalValue(5, 7, 3),
                                                                 IntervalValue(6, 8, 3),
                                                                 IntervalValue(7, 9, 3)]
                @test count(Ambiguous, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                   IntervalValue(2, 4, 0),
                                                                   IntervalValue(3, 5, 0),
                                                                   IntervalValue(4, 6, 0),
                                                                   IntervalValue(5, 7, 0),
                                                                   IntervalValue(6, 8, 0),
                                                                   IntervalValue(7, 9, 0)]
                @test count(Gap, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                             IntervalValue(2, 4, 0),
                                                             IntervalValue(3, 5, 0),
                                                             IntervalValue(4, 6, 0),
                                                             IntervalValue(5, 7, 0),
                                                             IntervalValue(6, 8, 0),
                                                             IntervalValue(7, 9, 0)]
                @test count(Match, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                               IntervalValue(2, 4, 3),
                                                               IntervalValue(3, 5, 3),
                                                               IntervalValue(4, 6, 3),
                                                               IntervalValue(5, 7, 2),
                                                               IntervalValue(6, 8, 1),
                                                               IntervalValue(7, 9, 0)]
                @test count(Mismatch, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                  IntervalValue(2, 4, 0),
                                                                  IntervalValue(3, 5, 0),
                                                                  IntervalValue(4, 6, 0),
                                                                  IntervalValue(5, 7, 1),
                                                                  IntervalValue(6, 8, 2),
                                                                  IntervalValue(7, 9, 3)]
            end
        end
        @testset "Mixed encodings" begin
            dnaA = dna"ATCGCCA-M"
            dnaB = BioSequence{DNAAlphabet{2}}("ATCGCCTAA")
            rnaA = rna"AUCGCCA-M"
            rnaB = BioSequence{RNAAlphabet{2}}("AUCGCCUAA")

            for seqs in ((dnaA, dnaB), (rnaA, rnaB))
                @test count(Certain, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                                 IntervalValue(2, 4, 3),
                                                                 IntervalValue(3, 5, 3),
                                                                 IntervalValue(4, 6, 3),
                                                                 IntervalValue(5, 7, 3),
                                                                 IntervalValue(6, 8, 2),
                                                                 IntervalValue(7, 9, 1)]
                @test count(Ambiguous, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                   IntervalValue(2, 4, 0),
                                                                   IntervalValue(3, 5, 0),
                                                                   IntervalValue(4, 6, 0),
                                                                   IntervalValue(5, 7, 0),
                                                                   IntervalValue(6, 8, 0),
                                                                   IntervalValue(7, 9, 1)]
                @test count(Gap, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                             IntervalValue(2, 4, 0),
                                                             IntervalValue(3, 5, 0),
                                                             IntervalValue(4, 6, 0),
                                                             IntervalValue(5, 7, 0),
                                                             IntervalValue(6, 8, 1),
                                                             IntervalValue(7, 9, 1)]
                @test count(Match, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 3),
                                                               IntervalValue(2, 4, 3),
                                                               IntervalValue(3, 5, 3),
                                                               IntervalValue(4, 6, 3),
                                                               IntervalValue(5, 7, 2),
                                                               IntervalValue(6, 8, 1),
                                                               IntervalValue(7, 9, 0)]
                @test count(Mismatch, seqs[1], seqs[2], 3, 1) == [IntervalValue(1, 3, 0),
                                                                  IntervalValue(2, 4, 0),
                                                                  IntervalValue(3, 5, 0),
                                                                  IntervalValue(4, 6, 0),
                                                                  IntervalValue(5, 7, 1),
                                                                  IntervalValue(6, 8, 2),
                                                                  IntervalValue(7, 9, 3)]
            end
        end
    end
end


text:
# this script basially only handles `ARGS`

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
try
    using VSCodeServer
catch err
    local distributed = Base.PkgId(Base.UUID("8ba89e20-285c-5b6f-9357-94700520ee1b"), "Distributed")
    if haskey(Base.loaded_modules, distributed)
        local Distributed = Base.loaded_modules[distributed]
        if err isa CompositeException
            local ex1 = first(err.exceptions)
            if ex1 isa Distributed.RemoteException && ex1.pid == 1
                rethrow()
            end
        else
            rethrow()
        end
    else
        rethrow()
    end
end
popfirst!(LOAD_PATH)

let
    args = [popfirst!(Base.ARGS) for _ in 1:6]
    # load Revise ?
    if "USE_REVISE=true" in args
        try
            @static if VERSION ≥ v"1.5"
                using Revise
            else
                @eval using Revise
                Revise.async_steal_repl_backend()
            end
        catch err
        end
    end

    atreplinit() do repl
        VSCodeServer.toggle_plot_pane(nothing, (;enable="USE_PLOTPANE=true" in args))
        VSCodeServer.toggle_progress(nothing, (;enable="USE_PROGRESS=true" in args))
    end

    conn_pipeline, telemetry_pipeline = args[1:2]
    VSCodeServer.serve(conn_pipeline; is_dev="DEBUG_MODE=true" in args, crashreporting_pipename=telemetry_pipeline)
end

max_stars_repo_name: novitk/julia-vscode

max_issues_repo_name: novitk/julia-vscode

max_forks_repo_name: novitk/julia-vscode


text:
##Specific problems from Kloden&Platen 1997
abstract type Problem end

struct Problem31 <: Problem
    avec::Vector{Float64}
end

function mgen(prob::Problem31, x0::AbstractVector)
    return length(x0)
end

function a!(prob::Problem31)
    coeff = -prob.avec.^2
    return function (t, x, f)
        f .= coeff.*(x.*(1-x.^2))
    end
end

function b!(prob::Problem31)
    return function (t, x, g)
        g .= diagm(prob.avec.*(1-x.^2))
    end
end

function bbp!(prob::Problem31)
    coeff = -2*prob.avec.^2
    return function (t, x, bbp)
        bbp .= coeff.*(x-x.^3)
    end
end

#Diffeq has a different convention from v4.0
function diffeq_a!(prob::Problem31)
    coeff = -prob.avec.^2
    return function (t, x, p, f)
        f .= coeff.*(x.*(1-x.^2))
    end
end

function diffeq_b!(prob::Problem31)
    return function (t, x, p, g)
        g .= diagm(prob.avec.*(1-x.^2))
    end
end

function solution(prob::Problem31, x0, t, Wvec)
    avec = prob.avec
    return [tanh.(avec.*W+atanh.(x0)) for W in Wvec]
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
_tupleme(x::Float64) = x < 0.0 ? error("Max delay must be positive") : (0.0, x)
_tupleme(x::Tuple{Float64, Float64}) = x[1] >= x[2] ? error("Specify delay as (min, max)") : x

struct EventExtents{T <: DiseaseStateSequence}
  exposure::Union{Nothing, Tuple{Float64, Float64}}
  infection::Union{Nothing, Tuple{Float64, Float64}}
  removal::Union{Nothing, Tuple{Float64, Float64}}
  EventExtents{SEIR}(e, i, r) = new{SEIR}(_tupleme(e), _tupleme(i), _tupleme(r))
  EventExtents{SEI}(e, i)     = new{SEI}(_tupleme(e), _tupleme(i), nothing)
  EventExtents{SIR}(i, r)     = new{SIR}(nothing, _tupleme(i), _tupleme(r))
  EventExtents{SI}(i)         = new{SI}(nothing, _tupleme(i), nothing)
end

function Base.show(io::IO, x::EventExtents{T}) where T <: DiseaseStateSequence
  print(io, "$T model event extents")
end


text:
function inner_craft(its::Vector,lin=3,col=3)
	# 载入
	mat=Matrix{String}(undef,lin,col) # 辨识符
	filln=0x0 # 非空数
	for i in 1:lin
		for j in 1:col
			k=i*(col-1)+j
			mat[i,j]=its[k]
			if !isa(mat[i,j],EI)
				filln+=1
			end
		end
	end
# 左上角(mov_i,mov_j)
#=
	fl=false;mov_i=1;mov_j=1
	while mov_i<lin
		for j in 1:col
			if !isa(mat[mov_lin,j],EI)
				fl=true
				break
			end
		end
		if fl break end
		mov_i+=1
	end
	fl=false
	while mov_i<col
		for i in mov_j:lin
			if !isa(mat[i,mov_j],EI)
				fl=true
				break
			end
		end
		if fl break end
		mov_j+=1
	end
=#
	pattern=Vector{UInt8}()
	das=Vector{Item}(undef,filln)
	ids=Vector{String}(undef,filln)
	nas=Vector{String}(undef,filln)
	cou=1
	for i in 1:lin
		v=0
		for j in 1:col
			if !isa(mat[i,j],EI)
				das[cou]=mat[i,j]
				ids[cou]=id(mat[i,j])
				nas[cou]=mat[i,j].name.name
				v+=1<<(j-1)
				cou+=1
			end
		end
		if v!=0
			push!(pattern,v)
		end
	end
	res=Vector{Pair{Item,UInt8}}()
	# 检查普通合成表
	# [[3,1],["a1","a2","a3"],"a",2]
	for c in data[:普通合成表]
		if c[1]==pattern&&c[2]==ids
			push!(res,Pair(guess_struct(c[3]),length(c)==4 ? c[4] : 1))
		end
	end
	# 检查类合成表
	# [[3,1],[]]
	for c in data[:类合成表]
		if c[1]!=pattern
			continue
		end
		fl=false
		tarr=Vector{Union{UInt8,Nothing}}(nothing,4)
		for k in 1:length(c[2])>>1
			reqid::UInt8=c[2][k<<1-1]
			req=c[2][k<<1]
			if req!=nas[k]
				fl=true
				break
			elseif reqid!=0x0
				if tarr[reqid]===nothing
					tarr[reqid]=das[k].t
				elseif tarr[reqid]!=das[k].t
					fl=true
					break
				end
			end
		end
		if !fl
			str=c[3]
			s=replace(str,"\$"=>tarr[1])
			push!(res,Pair(guess_struct(s),length(c)==3 ? 0x1 : c[4]))
		end
	end
	return res
end


text:
#==========================================================================================#

# UNIT WEIGHTS

struct UnitWeights{S<:Real, T<:Real} <: AbstractWeights{S, T, V where V<:Vector{T}}

    el::T
    sum::S

    UnitWeights{S, T}(el::Type{<:Real}, s::Real) where {S, T} = new(one(el), s)
end

UnitWeights(::T, s::S)       where {S, T} = UnitWeights{S, T}(T, s)
UnitWeights(::Type{T}, s::S) where {S, T} = UnitWeights{S, T}(T, s)

Base.eltype(wv::UnitWeights{S, T})  where {S, T} = T
Base.values(wv::UnitWeights{S, T})  where {S, T} = fill(wv.el, length(wv))
Base.sum(wv::UnitWeights{S, T})     where {S, T} = wv.sum
Base.isempty(wv::UnitWeights{S, T}) where {S, T} = iszero(wv.sum)
Base.length(wv::UnitWeights{S, T})  where {S, T} = Int(wv.sum)
Base.size(wv::UnitWeights{S, T})    where {S, T} = Tuple(length(wv))

function Base.getindex(wv::UnitWeights{S, T}, i::Int) where {S, T}
    (i > 0) & (i <= length(wv)) ? wv.el : Base.throw_boundserror(wv, i)
end

function Base.getindex(wv::UnitWeights{S, T}, i::AbstractRange{<:Int}) where {S, T}
    if (minimum(i) > 0) & (maximum(i) <= length(wv))
        return fill(wv.el, length(i))
    else
        Base.throw_boundserror(wv, i)
    end
end

function Base.getindex(wv::UnitWeights{S, T}, i::AbstractArray{<:Int}) where {S, T}
    if (minimum(i) > 0) & (maximum(i) <= length(wv))
        return fill(wv.el, size(i))
    else
        Base.throw_boundserror(wv, i)
    end
end

Base.getindex(wv::UnitWeights{S, T}, ::Colon) where {S, T} = fill(wv.el, length(wv))

sum(v::AbstractArray, ::UnitWeights)  = sum(v)
mean(v::AbstractArray, ::UnitWeights) = mean(v)

isequal(x::UnitWeights, y::UnitWeights) = (isequal(x.sum, y.sum) && isequal(x.el, y.el))


text:
#=
from "Error-free transformations in real and complex floating point arithmetics"
by Stef Graillat and Valérie Ménissier-Morain
2007 International Symposium on Nonlinear Theory and its Applications
NOLTA'07, Vancouver, Canada, September 16-19, 2007
=#
@inline function two_sum(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_sum(a.re, b.re)
    hi2, lo2 = two_sum(a.im, b.im)
    hi = Complex{T}(hi1, hi2)
    lo = Complex{T}(lo1, lo2)
    return hi, lo
end

@inline function one_sum(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_sum(a.re, b.re)
    hi2, lo2 = two_sum(a.im, b.im)
    hi = Complex{T}(hi1, hi2)
    return hi
end

two_hilo_sum(a::Complex{T}, b::Complex{T}) where {T<:Real} = two_sum(a, b)
one_hilo_sum(a::Complex{T}, b::Complex{T}) where {T<:Real} = one_sum(a, b)

@inline function two_diff(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_diff(a.re, b.re)
    hi2, lo2 = two_diff(a.im, b.im)
    hi = Complex{T}(hi1, hi2)
    lo = Complex{T}(lo1, lo2)
    return hi, lo
end

@inline function one_diff(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_diff(a.re, b.re)
    hi2, lo2 = two_diff(a.im, b.im)
    hi = Complex{T}(hi1, hi2)
    return hi
end

two_hilo_diff(a::Complex{T}, b::Complex{T}) where {T<:Real} = two_diff(a, b)
one_hilo_diff(a::Complex{T}, b::Complex{T}) where {T<:Real} = one_diff(a, b)

#=
from "Error-free transformations in real and complex floating point arithmetics"
by Stef Graillat and Valérie Ménissier-Morain
2007 International Symposium on Nonlinear Theory and its Applications
NOLTA'07, Vancouver, Canada, September 16-19, 2007
N.B. modified to order q,r,s
=#
@inline function four_prod(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_prod(a.re, b.re)
    hi2, lo2 = two_prod(a.im, b.im)
    hi3, lo3 = two_prod(a.re, b.im)
    hi4, lo4 = two_prod(a.im, b.re)
    hi5, lo5 = two_diff(hi1, hi2)
    hi6, lo6 = two_sum(hi3, hi4)
    p = Complex{T}(hi5, hi6)
    q = Complex{T}(lo5, lo6)
    r = Complex{T}(lo1, lo3)
    s = Complex{T}(-lo2, lo4)
    r, s = two_sum(r, s)
    q, r = two_sum(q, r)
    return p, q, r, s
end

@inline function three_prod(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_prod(a.re, b.re)
    hi2, lo2 = two_prod(a.im, b.im)
    hi3, lo3 = two_prod(a.re, b.im)
    hi4, lo4 = two_prod(a.im, b.re)
    hi5, lo5 = two_diff(hi1, hi2)
    hi6, lo6 = two_sum(hi3, hi4)
    p = Complex{T}(hi5, hi6)
    q = Complex{T}(lo5, lo6)
    r = Complex{T}(lo1, lo3)
    s = Complex{T}(-lo2, lo4)
    q, r = two_sum(q, r, s)
    return p, q, r
end

@inline function two_prod(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_prod(a.re, b.re)
    hi2, lo2 = two_prod(a.im, b.im)
    hi3, lo3 = two_prod(a.re, b.im)
    hi4, lo4 = two_prod(a.im, b.re)
    hi5, lo5 = two_diff(hi1, hi2)
    hi6, lo6 = two_sum(hi3, hi4)
    p = Complex{T}(hi5, hi6)
    re = one_sum(lo5, lo1, lo2)
    im = one_sum(lo6, lo3, lo4)
    q = Complex{T}(re, im)
    return p, q
end

@inline function one_prod(a::Complex{T}, b::Complex{T}) where {T<:Real}
    hi1, lo1 = two_prod(a.re, b.re)
    hi2, lo2 = two_prod(a.im, b.im)
    hi3, lo3 = two_prod(a.re, b.im)
    hi4, lo4 = two_prod(a.im, b.re)
    hi5, lo5 = two_diff(hi1, hi2)
    hi6 = hi3 + hi4
    p = Complex{T}(hi5, hi6)
    return p
end


text:
abstract type AbstractState end

struct EmptyState <: AbstractState end

"""
    GeneralState{T <: Number}
A mutable structure to store the general state of a solver. All methods applied to `GeneralState` are forwarded to `S <: AbstractSolver`.
```
mutable struct GeneralState{T <: Number} <: AbstractState
    configuration::Configuration{T}
    cons_costs::Dictionary{Int, Float64}
    last_improvement::Int
    tabu::Dictionary{Int, Int}
    vars_costs::Dictionary{Int, Float64}
end
```
"""
mutable struct _State{T} <: AbstractState
    configuration::Configuration{T}
    cons_costs::Dictionary{Int, Float64}
    fluct::Fluct
    icn_computations::Matrix{Float64}
    optimizing::Bool
    last_improvement::Int
    vars_costs::Dictionary{Int, Float64}
end

@forward _State.configuration get_values, get_error, get_value, compute_cost!, set_values!
@forward _State.configuration set_value!, set_sat!

const State = Union{EmptyState, _State}

state() = EmptyState()
function state(m::_Model, pool = pool(); opt = false)
    X = Matrix{Float64}(undef, m.max_vars[], CompositionalNetworks.max_icn_length())
    lc, lv = length_cons(m) > 0, length_vars(m) > 0
    config = Configuration(m, X)
    cons = lc ? zeros(Float64, get_constraints(m)) : Dictionary{Int,Float64}()
    last_improvement = 0
    vars = lv ? zeros(Float64, get_variables(m)) : Dictionary{Int,Float64}()
    fluct = Fluct(cons, vars)
    return _State(config, cons, fluct, X, opt, last_improvement, vars)
end

"""
    _cons_costs(s::S) where S <: Union{_State, AbstractSolver}
Access the constraints costs.
"""
_cons_costs(s::_State) = s.cons_costs

"""
    _vars_costs(s::S) where S <: Union{_State, AbstractSolver}
Access the variables costs.
"""
_vars_costs(s::_State) = s.vars_costs

"""
    _vars_costs(s::S) where S <: Union{_State, AbstractSolver}
Access the variables costs.
"""
_values(s::_State) = get_values(s)

"""
    _optimizing(s::S) where S <: Union{_State, AbstractSolver}
Check if `s` is in an optimizing state.
"""
_optimizing(s::_State) = s.optimizing

"""
    _cons_costs!(s::S, costs) where S <: Union{_State, AbstractSolver}
Set the constraints costs.
"""
_cons_costs!(s::_State, costs) = s.cons_costs = costs

"""
    _vars_costs!(s::S, costs) where S <: Union{_State, AbstractSolver}
Set the variables costs.
"""
_vars_costs!(s::_State, costs) = s.vars_costs = costs

"""
    _values!(s::S, values) where S <: Union{_State, AbstractSolver}
Set the variables values.
"""
_values!(s::_State{T}, values) where T <: Number = set_values!(s, values)

"""
    _optimizing!(s::S) where S <: Union{_State, AbstractSolver}
Set the solver `optimizing` status to `true`.
"""
_optimizing!(s::_State) = s.optimizing = true

"""
    _satisfying!(s::S) where S <: Union{_State, AbstractSolver}
Set the solver `optimizing` status to `false`.
"""
_satisfying!(s::_State) = s.optimizing = false

"""
    _cons_cost(s::S, c) where S <: Union{_State, AbstractSolver}
Return the cost of constraint `c`.
"""
_cons_cost(s::_State, c) = get!(_cons_costs(s), c, 0.0)

"""
    _var_cost(s::S, x) where S <: Union{_State, AbstractSolver}
Return the cost of variable `x`.
"""
_var_cost(s::_State, x) = get!(_vars_costs(s), x, 0.0)

"""
    _value(s::S, x) where S <: Union{_State, AbstractSolver}
Return the value of variable `x`.
"""
_value(s::_State, x) = _values(s)[x]

"""
    _cons_cost!(s::S, c, cost) where S <: Union{_State, AbstractSolver}
Set the `cost` of constraint `c`.
"""
_cons_cost!(s::_State, c, cost) = _cons_costs(s)[c] = cost

"""
    _var_cost!(s::S, x, cost) where S <: Union{_State, AbstractSolver}
Set the `cost` of variable `x`.
"""
_var_cost!(s::_State, x, cost) = _vars_costs(s)[x] = cost

"""
    _value!(s::S, x, val) where S <: Union{_State, AbstractSolver}
Set the value of variable `x` to `val`.
"""
_value!(s::_State, x, val) = _values(s)[x] = val

"""
    _set!(s::S, x, val) where S <: Union{_State, AbstractSolver}
Set the value of variable `x` to `val`.
"""
_set!(s::_State, x, val) = set!(_values(s), x, val)

"""
    _set!(s::S, x, y) where S <: Union{_State, AbstractSolver}
Swap the values of variables `x` and `y`.
"""
function _swap_value!(s::_State, x, y)
    aux = _value(s, x)
    _value!(s, x, _value(s, y))
    _value!(s, y, aux)
end

_last_improvement(s::_State) = s.last_improvement
_inc_last_improvement!(s::_State) = s.last_improvement += 1
_reset_last_improvement!(s::_State) = s.last_improvement = 0

has_solution(s::_State) = is_solution(s.configuration)

function set_error!(s::_State, err)
    sat = err ≈ 0.0
    set_sat!(s, sat)
    !sat && set_value!(s, err)
end

get_error(::EmptyState) = Inf


text:
module Vulnerability
using ExternalBallistics
using OpticSim
using LinearAlgebra
# Write your package code here.
include("trajectories.jl")
include("targets.jl")
include("shotlines.jl")

export raytracing, target, shotlines, targetPosition

end


text:
n =  7 # rand(1:7)
A = matrixdepot("invol", n)
# A*A ≈ eye(n)
@test A*A ≈ Matrix(1.0I, n, n) atol = 1e-5
println("'invol' passed test...")


text:
# args: dictfile outfile < input

using HDF5,JLD # ,DataStructures
using KUparser: Sentence, Position, DepRel

# Store postag/deprel strings as positive integer values
postag_hash = Dict{UTF8String,DepRel}()
deprel_hash = Dict{UTF8String,DepRel}()

# First argument is a dictionary file, could be an existing jld file
dict = load(ARGS[1])
for i=1:length(dict["postag"]); postag_hash[dict["postag"][i]] = i; end
for i=1:length(dict["deprel"]); deprel_hash[dict["deprel"][i]] = i; end

# ROOT deprel is stored as the special value 0
const ROOT="ROOT"
deprel_hash[ROOT] = 0

# The input is in conll07 format followed by wvec+cvec and read from stdin:
# ID,FORM,LEMMA,CPOSTAG,POSTAG,FEATS,HEAD,DEPREL,PHEAD,PDEPREL

corpus  = Sentence[]
s = nothing
mydims = 0
for l in eachline(STDIN)
    if l == "\n"
        @assert s != nothing
        push!(corpus, s)
        s = nothing
        continue
    end
    (id, form, lemma, cpostag, postag, feats, head, deprel, phead, pdeprel, wvec) = split(chomp(l), '\t')
    form = convert(UTF8String, form)
    postag = postag_hash[postag]
    deprel = deprel_hash[deprel]
    head = convert(Position, int(head))
    wvec = map(float32, split(wvec, ' '))
    if s == nothing
        mydims = length(wvec)
        s = Sentence()
        s.form = [form]
        s.postag = [postag]
        s.head = [head]
        s.deprel = [deprel]
        s.wvec = wvec''
    else
        @assert mydims == length(wvec)
        push!(s.form, form)
        push!(s.postag, postag)
        push!(s.head, head)
        push!(s.deprel, deprel)
        s.wvec = hcat(s.wvec, wvec)
    end
end

@assert s == nothing
@assert length(corpus) > 0
save(ARGS[2], "corpus", corpus, "postag", dict["postag"], "deprel", dict["deprel"])


text:
using SparseArrays

import Random: shuffle!

struct SkipGramData{A<:AbstractArray{<:Tuple{<:Integer, <:Integer, <:Integer}, 1}, B<:AbstractArray{<:Context{<:Integer, <:AbstractArray{<:Integer,1}}, 1}, C<:Integer, D<:AbstractDict{<:Integer, <:AbstractFloat}} <: TextData
    ordering::A # to implement with iteration interface
    xs::B # the contexts and words to iterate over
    n::C # the size of the vocabulary
    freq_table::D
end

function Base.push!(data::SkipGramData, c::Context)
    push!(data.xs, c)
    position = length(data.xs)
    num_contexts = length(c.xs)
    for i in 1:num_contexts
        words = length(iterate(c, i)[1][2])
        for j in 1:words
            push!(data.ordering, (position, i, j))
        end
    end
end

SkipGramData(n::Int, freq_table::Dict) = SkipGramData(Tuple{Int, Int, Int}[],Context{Int, Array{Int,1}}[],n, freq_table)

function shuffle!(data::SkipGramData)
    shuffle!(data.ordering)
    return nothing
end

function Base.length(data::SkipGramData)
    return length(data.ordering)
end

function Base.iterate(data::SkipGramData, state = 1)
    if state < 1 || state > length(data) + 1
        error("Invalid iterator state parameter")
    elseif state == length(data) + 1
        return nothing
    else
        (contextid, context_state, word_id) = data.ordering[state]
        (context, _) = iterate(data.xs[contextid], context_state)
        x = SparseVector(data.n, Int[context[1]], [1])
        y = iterate(context[2], word_id)[1]
        return ((x, y), state + 1)
    end
end


text:
# this is for the metric computation
"""
    MajorityVoteClassifier(x, y)

Construct it using training inputs x and outputs y.
"""
struct MajorityVoteClassifier
    x::Vector
    y::Vector
end

function predict(m::MajorityVoteClassifier, x::Number)
    yx = m.y[m.x .== x]
    if length(yx) == 0
        return sample(unique(m.y), 1)[1]
    end
    cm = countmap(yx)
    collect(keys(cm))[argmax(collect(values(cm)))]
end
predict(m::MajorityVoteClassifier, x::Vector) = map(_x->predict(m,_x), x)
error_rate(y_true::Vector, y_pred::Vector) = mean(y_pred .!= y_true)
precision(y_true::Vector, y_pred::Vector) = mean(y_pred .== y_true)

function split_pairs(all_pairs, train_ratio=0.5; seed=nothing)
    tr_inds, val_inds, tst_inds =
        HierarchicalAD.train_val_test_inds(1:length(all_pairs), (train_ratio,0,1-train_ratio); seed=seed)
    tr_x = [_x[1] for _x in all_pairs[tr_inds]]
    tr_y = [_x[2] for _x in all_pairs[tr_inds]]
    tst_x = [_x[1] for _x in all_pairs[tst_inds]]
    tst_y = [_x[2] for _x in all_pairs[tst_inds]]
    (tr_x, tr_y), (tst_x, tst_y)
end

function precision_using_mvc(tr_x, tr_y, tst_x, tst_y)
    mvc = MajorityVoteClassifier(tr_x, tr_y);
    y_pred = predict(mvc, tst_x)
    precision(y_pred, tst_y)
end

# this is for disentanglement per dim
function get_least_varied_dim(x::AbstractArray{T,2}, sds) where T
    x = x ./ sds
    v = var(x, dims=2)
    argmin(v)[1], v
end

function least_varied_dim_pairs(z::AbstractArray{T,2}, y, factors;
	samples_per_factor=100, batchsize=2000) where T
	# first compute sd per dimension
    sds = map(std, eachrow(z))

    # now loop over factors and latents
    all_pairs = []
    vars = []
    for k in 1:length(factors)
        factor = factors[k]
        fvals = unique(y[!,factor])

        k_pairs = []
        for i in 1:samples_per_factor
            fk = sample(fvals)
            zk = sample_tensor(z[:,y[!,factor] .== fk], batchsize)
            ik, vs = get_least_varied_dim(zk, sds)
            push!(k_pairs, (ik, k))
            push!(vars, vs)
        end
        push!(all_pairs, k_pairs)
    end
    # now this contains the pairs with (x=least varied latent index, y=factor index)
    all_pairs = vcat(all_pairs...)
    all_pairs, vars
end

"""
	 disentanglement_per_dim(
	 	z::AbstractArray{T,2}, y, factors;
        samples_per_factor=100, batchsize=2000
        )
"""
function disentanglement_per_dim(z::AbstractArray{T,2}, y, factors; kwargs...) where T
    # get the training pairs
    all_pairs, vars = least_varied_dim_pairs(z, y, factors; kwargs...)

    # here, train and predict using majority voting classifier
	(tr_x, tr_y), (tst_x, tst_y) = split_pairs(all_pairs, 0.5)
	precision_using_mvc(tr_x, tr_y, tst_x, tst_y)
end
disentanglement_per_dim(z::AbstractArray{T,2}, y, factors, M::Int; kwargs...) where T =
	mean(map(_->disentanglement_per_dim(z, y, factors; kwargs...), 1:M))
disentanglement_per_dim(z::AbstractArray{T,2}, y::Nothing, args...; kwargs...) where T = nothing

# this is for disentanglement per latent
function get_least_varied_latent(zs::Vector, sds_per_latent) where T
    dcovs = []
    for latenti in 1:length(zs)
        z = zs[latenti]
        nz = z ./ sds_per_latent[latenti]
        push!(dcovs, det(cov(nz, dims=2)))
    end
    argmin(dcovs)[1], dcovs
end

function least_varied_latent_pairs(zs::Vector, y, factors; samples_per_factor=100, batchsize=2000)
	# first compute sd per dimension
    sds_per_latent = map(z->map(std, eachrow(z)),zs)

    # now loop over factors and latents
    all_pairs = []
    covs = []
    for k in 1:length(factors)
        factor = factors[k]
        fvals = unique(y[!,factor])

        k_pairs = []
        for i in 1:samples_per_factor
        	fk = sample(fvals)
            zks = []
            for latenti in 1:length(zs)
                push!(zks, sample_tensor(zs[latenti][:,y[!,factor] .== fk], batchsize))
            end
            ik, dcovs = get_least_varied_latent(zks, sds_per_latent)
            push!(k_pairs, (ik, k))
            push!(covs, dcovs)
        end
        push!(all_pairs, k_pairs)
    end
    # now this contains the pairs with (x=least varied latent index, y=factor index)
    all_pairs = vcat(all_pairs...)
    all_pairs, covs
end

"""
	disentanglement_per_latent(zs::Vector, y, factors; samples_per_factor=100, batchsize=2000)
"""
function disentanglement_per_latent(zs::Vector, y, factors; kwargs...)
    # get the training pairs
    all_pairs, covs = least_varied_latent_pairs(zs, y, factors; kwargs...)

    # here, train and predict using majority voting classifier
    (tr_x, tr_y), (tst_x, tst_y) = split_pairs(all_pairs, 0.5)
    precision_using_mvc(tr_x, tr_y, tst_x, tst_y)
end
disentanglement_per_latent(zs::Vector, y, factors, M::Int; kwargs...) =
	mean(map(_->disentanglement_per_latent(zs, y, factors; kwargs...), 1:M))
disentanglement_per_latent(zs::Vector, y::Nothing, args...; kwargs...) = nothing

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
using DelimitedFiles

function solve(input)
    raw = readdlm(input,'\n')
    coords =vcat(map(toInt,raw)...)
    @show coords
    num_points = length(coords[:,1])
    @show num_points
    x_min, x_max = extrema(coords[:,1])
    y_min, y_max = extrema(coords[:,2])
    @show (x_min, x_max, y_min, y_max)
    finite_indices = Set(collect(1:num_points))

    safe=zeros(Int16,num_points)
    best=zeros(Int16,num_points)

    for x in x_min:x_max
        for y in y_min:y_max
            distances = sum(abs.([x y] .- coords),dims=2)
            #@show (x,y,distances)
            #@show sum(distances)
            mindis = minimum(distances)
            allmins = findall(d->d==mindis,distances)
            if (length(allmins)==1)
                safe[allmins[1]]+=1
                if x == x_min || x==x_max || y==y_min || y==y_min
                    delete!(finite_indices,allmins[1][1])
                end

            end
            if sum(distances) < 10000
                best[allmins[1]]+=1
            end
        end
    end
    @show safe
    @show best
    @show finite_indices
    part1 = maximum(safe[collect(finite_indices)])
    part2 = sum(best)
    @show part1
    @show part2
end

function toInt(s)
    m = match(r"(?<x>\d+), (?<y>\d+)", s)
    x = tryparse(Int16,m[:x])
    y = tryparse(Int16,m[:y])
    return [x y]
end

solve("input.txt")


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVR

using Test, RLinearAlgebra, Random

Random.seed!(1010)

@testset "LSVR Hop Random Cyclic -- Procedural" begin

    # Verify appropraite super type
    @test supertype(LinSysVecRowHopRandCyclic) == LinSysVecRowSampler

    # Test whether cycles are changed after appropriate hop period
    for k = 1:10
        @test let
            A = rand(10,3)
            b = rand(10)
            x = rand(3)

            cyc = LinSysVecRowHopRandCyclic(k)

            α, β = RLinearAlgebra.sample(cyc, A, b, x, 1)

            order = copy(cyc.order)

            flag = true
            for j = 2:(cyc.hop * 10)
                α, β = RLinearAlgebra.sample(cyc, A, b, x, j)

                # Ordering should not change
                flag = flag & (cyc.order == order)
            end

            α, β = RLinearAlgebra.sample(cyc, A, b, x, cyc.hop * 10 + 1)

            # Ordering should change
            flag = flag & (cyc.order != order)

            flag
        end
    end
end

end