text:
function quadrotor_benchmarks!(suite::BenchmarkGroup)

    suite["line"] = BenchmarkGroup()
    # suite["maze"] = BenchmarkGroup()

    T = Float64

    # options
    max_con_viol = 1.0e-8
    verbose=false

    opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
        iterations=300)

    opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
        opts_uncon=opts_ilqr,
        iterations=40,
        cost_tolerance=1.0e-5,
        cost_tolerance_intermediate=1.0e-4,
        constraint_tolerance=max_con_viol,
        penalty_scaling=10.,
        penalty_initial=1.)

    opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
        feasibility_tolerance=max_con_viol,
        solve_type=:feasible)

    opts_altro = ALTROSolverOptions{T}(verbose=verbose,
        opts_al=opts_al,
        R_inf=1.0e-8,
        resolve_feasible_problem=false,
        opts_pn=opts_pn,
        projected_newton=true,
        projected_newton_tolerance=1.0e-3)

    opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=Ipopt.Optimizer(),
        opts=Dict(:max_iter=>10000),
        feasibility_tolerance=max_con_viol)

    opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=SNOPT7.Optimizer(),
        feasibility_tolerance=max_con_viol,
        opts=Dict(:Iterations_limit=>500000,
            :Major_iterations_limit=>1000))
    run_benchmarks!(suite, Problems.quadrotor, [opts_ilqr, opts_al, opts_altro, opts_ipopt, opts_snopt])

    # Quadrotor in Maze
    opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
        iterations=300)

    opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
        opts_uncon=opts_ilqr,
        iterations=40,
        cost_tolerance=1.0e-5,
        cost_tolerance_intermediate=1.0e-4,
        constraint_tolerance=max_con_viol,
        penalty_scaling=10.,
        penalty_initial=1.)

    opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
        feasibility_tolerance=max_con_viol,
        solve_type=:feasible)

    opts_altro = ALTROSolverOptions{T}(verbose=verbose,
        opts_al=opts_al,
        R_inf=1.0e-8,
        resolve_feasible_problem=false,
        opts_pn=opts_pn,
        projected_newton=true,
        projected_newton_tolerance=1.0e-4)

    opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=Ipopt.Optimizer(),
        opts=Dict(:max_iter=>10000),
        feasibility_tolerance=1.0e-3)

    opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
        nlp=SNOPT7.Optimizer(),
        feasibility_tolerance=1.0e-3,
        opts=Dict(:Iterations_limit=>500000,
            :Major_iterations_limit=>1000))
    # run_benchmarks!(suite["maze"], Problems.quadrotor_maze, [opts_altro,])

end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
module LockandKeyLookups

    include("Lookups.jl")
    export LockandKeyLookup

end # module

ext: jl


text:
using LambdaCalculus
using Test

const N = LambdaCalculus.Named
const D = LambdaCalculus.DeBruijn
const LN = LambdaCalculus.LocallyNameless

const Log = Base.CoreLogging
Log.global_logger(Log.SimpleLogger(stderr, Log.Debug))

# MACROS
@testset "Macros" begin
    @testset "Named" begin
        id = N.@λ x -> x
        l = N.@λ $id(a)
        r = N.@λ (y -> y)(a)
        @test l ≃ r

        @test freevars(N.@λ x -> y) == Set([:y])
        @test (N.@λ x) ≄ (N.@lambda y)

        # interpolation of named terms is _not_ hygienic
        t = N.@lambda x -> y
        @test (N.@lambda y -> $t) ≄ (N.@lambda z1 -> z2 -> y)
        @test (N.@lambda y -> $t) ≃ (N.@lambda z1 -> z2 -> z1)
    end

    @testset "DeBruijn" begin
        id = D.@λ x -> x
        l = D.@λ $id(a)
        r = D.@λ (y -> y)(a)
        @test l ≃ r

        @test freevars(D.@λ x -> y) == Set([1])
        @test (D.@λ [x, y] x) ≄ (D.@lambda [x, y] y)

        # interpolation of de Bruijn terms _is_ hygienic (inserted terms are shifted)
        t = D.@lambda x -> y
        @test (D.@lambda y -> $t) ≃ (D.@lambda z1 -> z2 -> y)
    end

    @testset "LocallyNameless" begin
        id = LN.@λ x -> x
        l = LN.@λ $id(a)
        r = LN.@λ (y -> y)(a)
        @test l ≃ r

        @test freevars(LN.@λ x -> y) == Set([:y])
        @test (LN.@λ x) ≄ (LN.@lambda y)

        # interpolation of locally nameless terms _is_ hygienic (inserted terms are shifted)
        t = LN.@lambda x -> y
        @test (LN.@lambda y -> $t) ≃ (LN.@lambda z1 -> z2 -> y)
    end
end

@testset "Conversion" begin
    Γ = NamingContext([:x, :y, :z])

    named_lambdas = [(N.@lambda x -> x),
                     (N.@lambda (x -> x(x))(x -> x(x))),
                     (N.@lambda x -> y -> x),
                     (N.@lambda x -> y)]

    debruijn_lambdas = [(D.@lambda x -> x),
                        (D.@lambda (x -> x(x))(x -> x(x))),
                        (D.@lambda x -> y -> x),
                        (D.@lambda [x, y, z] x -> y)]

    ln_lambdas = [(LN.@lambda x -> x),
                  (LN.@lambda (x -> x(x))(x -> x(x))),
                  (LN.@lambda x -> y -> x),
                  (LN.@lambda x -> y)]

    # this used to cause a shadowing problem:
    @test convert(D.Term, N.@lambda((x -> (x -> x)))) ≃ D.@lambda x -> (y -> y)

    # conversion back to self
    for t in named_lambdas
        @test convert(N.Term, convert(D.Term, t, Γ), Γ) ≃ t
        @test convert(N.Term, convert(LN.Term, t, Γ), Γ) ≃ t
    end

    for t in debruijn_lambdas
        @test convert(D.Term, convert(N.Term, t, Γ), Γ) ≃ t
        @test convert(D.Term, convert(LN.Term, t, Γ), Γ) ≃ t
    end

    for t in ln_lambdas
        @test convert(LN.Term, convert(D.Term, t, Γ), Γ) ≃ t
        @test convert(LN.Term, convert(N.Term, t, Γ), Γ) ≃ t
    end

    # conversion to other formats
    for (tn, td, tl) in zip(named_lambdas, debruijn_lambdas, ln_lambdas)
        @test convert(D.Term, tn, Γ) ≃ td
        @test convert(LN.Term, tn, Γ) ≃ tl

        @test convert(N.Term, td, Γ) ≃ tn
        @test convert(LN.Term, td, Γ) ≃ tl

        @test convert(N.Term, tl, Γ) ≃ tn
        @test convert(D.Term, tl, Γ) ≃ td
    end
end

@testset "Evaluation" begin
    @testset "DeBruijn" begin
        terms = [D.@lambda((x -> (x -> x))(z -> z)),
                 D.@lambda((x -> x)(z -> (x -> x)(z))),
                 D.@lambda(y -> (f -> x -> f(x))(y)),
                 D.@lambda((f -> x -> f(x))(f -> x -> f(x))),
                 D.@lambda((f -> x -> f(f(x)))(f -> x -> f(f(x)))),
                 D.@lambda((x -> x)(y))]
        results = [D.@lambda(x -> x),
                   D.@lambda(z -> z),
                   D.@lambda(f -> x -> f(x)),
                   D.@lambda(f -> x -> f(x)),
                   D.@lambda(x -> y -> x(x(x(x(y))))),
                   D.@lambda(y)]

        for (t, r) in zip(terms, results)
            @test evaluate(t, 100) ≃ r
        end
    end

    @testset "LocallyNameless" begin
        terms = [LN.@lambda((x -> (x -> x))(z -> z)),
                 LN.@lambda((x -> x)(z -> (x -> x)(z))),
                 LN.@lambda(y -> (f -> x -> f(x))(y)),
                 LN.@lambda((f -> x -> f(x))(f -> x -> f(x))),
                 LN.@lambda((f -> x -> f(f(x)))(f -> x -> f(f(x)))),
                 LN.@lambda((x -> x)(y))]
        results = [LN.@lambda(x -> x),
                   LN.@lambda(z -> z),
                   LN.@lambda(f -> x -> f(x)),
                   LN.@lambda(f -> x -> f(x)),
                   LN.@lambda(x -> y -> x(x(x(x(y))))),
                   LN.@lambda(y)]

        for (t, r) in zip(terms, results)
            @test evaluate(t, 100) ≃ r
        end
    end
end

@testset "Sampling" begin
    for t in rand(D.BoundedTermSampler(20, 30), 100)
        @test 20 ≤ length(t) ≤ 30
    end
end


text:
module FisherySim

using Distributions
using TweedieDistributions
using StatsBase
using LinearAlgebra
using Random
using PDMats
using Arpack
using NeutralLandscapes

import Base:
    rand,
    +, -, *,
    step,
    sum,
    getindex, setindex!,
    size, length, eachindex,
    copy
import StatsBase:
    sample,
    cov
import Distributions:
    location

include("fisherydomain.jl")
export
    AbstractFisheryDomain,
    DiscreteFisheryDomain,
    GriddedFisheryDomain,
    size,
    length,
    eachindex,
    sample

include("covkernels.jl")
export
    AbstractCovarianceKernel,
    ExpCov,
    Matérn32Cov,
    Matern32Cov,
    AR1,
    cov

# include("matrixlognormal.jl")
# export
#     MatrixLogNormal,
#     location

include("domaindistributions.jl")
export
    AbstractDomainDistribution,
    DomainDistribution,
    MultiDomainDistribution,
    BlendedDomainDistribution,
    ClassifiedDomainDistribution,
    domain,
    getindex,
    length

include("habitat.jl")
export
    Habitat,
    getindex,
    length,
    HabitatPreference

include("bathymetry.jl")
export
    BathymetryModel,
    Bathymetry,
    rand

include("pop_dynamics.jl")
export
    PopulationDynamicsModel,
    PopState,
    Schaefer,
    PellaTomlinson,
    StochasticProduction,
    vecstate,
    step,
    sum,
    setindex!,
    copy

include("movement.jl")
export
    MovementModel,
    eqdist,
    MovementRate

include("targeting.jl")
export
    AbstractTargetingBehavior,
    RandomTargeting,
    FixedTargeting,
    AbstractPreferentialTargeting,
    PreferentialTargeting,
    DynamicPreferentialTargeting,
    target,
    reset!

include("catchability.jl")
export
    AbstractCatchability,
    Catchability,
    DensityDependentCatchability,
    HabitatCatchability,
    *

include("vessels.jl")
export
    Vessel,
    Catch,
    CPUE,
    +, -,
    fish!,
    getindex,
    Fleet,
    vessels

include("simulation.jl")
export
    simulate

end # module


text:
struct OpConnection{A<:AbstractArray,B<:AbstractArray,C<:AbstractArray}
    mel::A
    changes::Vector{StateChanges{B,C}}
end

function OpConnection(op::AbsLinearOperator)
    t = conn_type(op)
    v1, v2, v3 = t.parameters
    return OpConnection(v1(), Vector{StateChanges{v2,v3}}())
end

function OpConnection{A,B,C}() where {A,B,C}
    return OpConnection(A(), Vector{typeof(StateChanges(B(),C()))}())
end

function Base.resize!(c::OpConnection, n)
    resize!(c.mel,     n)
    resize!(c.changes, n)
    return c
end

Base.eltype(c::OpConnection{A,B,C}) where {A,B,C} = (eltype(A), eltype(B), eltype(C))
Base.length(c::OpConnection) = length(c.mel)
Base.size(c::OpConnection) = size(c.mel)

Base.:(==)(a::OpConnection, b::OpConnection) = (a.to_change == b.to_change &&
                                              a.new_values == b.new_values)

function Base.push!(c::OpConnection, (m_els, to_change, new_values))
    push!(c.mel, m_els)
    push!(c.changes, StateChanges(to_change, new_values))
    return c
end

function Base.push!(c::OpConnection, (m_els, cngs)::Tuple{Number, StateChanges})
    push!(c.mel, m_els)
    push!(c.changes, cngs)
    return c
end

function Base.append!(c::OpConnection, (m_els, to_change, new_values))
    append!(c.mel, m_els)
    cngs = Vector{StateChanges{eltype(to_change), eltype(new_values)}}()
    for (tcng, nwvls) = zip(to_change, new_values)
        push!(cngs, StateChanges(tcng, nwvls))
    end
    append!(c.changes, cngs)
    return c
end

function Base.append!(c::OpConnection, (m_els, cngs)::Tuple{Vector{<:Number},Vector{StateChanges}})
    append!(c.mel, m_els)
    append!(c.changes, cngs)
    return c
end

function Base.append!(c::OpConnection, new::OpConnection)
    append!(c.mel, new.mel)
    append!(c.changes, new.changes)
    return c
end

Base.getindex(c::OpConnection, i) = (c.mel[i], c.changes[i])

function Base.iterate(iter::OpConnection, state=1)
    state > length(iter) && return nothing
    return (iter[state], state + 1)
end

function add!(c::OpConnection, mel, cngs::StateChanges)
    i = findfirst(isequal(cngs), c.changes)
    isnothing(i) ? push!(c, (mel, cngs)) : c.mel[i] += mel
    return c
end

function add!(c::OpConnection, o::OpConnection)
    for (mel, cngs)=o
        add!(c, mel, cngs)
    end
end

Base.conj!(op::OpConnection) = conj!(op.mel)

function clear_duplicates(c::OpConnection)
    u_cngs = unique(c.changes)
    u_mel  = similar(c.mel, 0)
    for el=u_cngs
        ids = findall(isequal(el), c.changes)
        push!(u_mel, sum(c.mel[ids]))
    end
    return OpConnection(u_mel, u_cngs)
end


text:
# Chance Sampled Monte Carlo Counterfactual Regret Minimization
# Based on Johanson et al 2012 Algorithm 1
# http://modelai.gettysburg.edu/2013/cfr/cfr.pdf page 12 is also useful
# https://justinsermeno.com/posts/cfr/ is also useful for vanilla CFR

# h is represented by just a string, not an entire node structure
# this lets us avoid building the entire game tree and instead only have
# memory structures in the size of the info sets

# TODO: Dict is approx 10x faster than findfirst...restructure code?

# add packages
using Pkg
Pkg.activate(pwd())
using Revise, Random, StatsBase

const N = 1:2 # set of players
const cards = ["J","Q","K"] # cards in the deck
const hands = ["JQ", "JK", "QJ", "QK", "KJ", "KQ"]
const acts = ["p","b"]
const ni_stage = [1, 3, 6, 3] # ni_stage[j] is num of infosets at stage j
const ni = sum(ni_stage) # total number of infosets
const na = vcat(6, fill(2, ni - 1)) # na[j] is number of actions at info set j
const AI = vcat([hands], fill(acts, ni - 1)) # AI[i][j] is the jth action at infoset i
const leafplays = ["pp", "pbp", "pbb", "bp", "bb"]
const leafnodes = reshape([string(i, j) for j in leafplays, i in hands], 30)
const leafvals = [-1.0, -1, -2, 1, -2,
            -1, -1, -2, 1, -2,
            1, -1, 2, 1, 2,
            -1, -1, -2, 1, -2,
            1, -1, 2, 1, 2,
            1, -1, 2, 1, 2]
const Uz = [leafvals, -leafvals]
const Iset = ["",
    "J", "Q", "K",
    "Jp", "Qp", "Kp", "Jb", "Qb", "Kb",
    "Jpb", "Qpb", "Kpb"]

actions(h::String) = length(h) < 1 ? hands : ["p", "b"] # node actions

function player(h::String)
    len = length(h)
    if len < 2
        return 3
    else
        return isodd(length(h) - 2) ? 2 : 1
    end
end

function infoset(h::String)
    depth = length(h)
    if depth < 2
        return 1 # null infoset
    elseif depth == 2
        return findfirst(cards .== h[1:1])[1] + ni_stage[1]
    elseif depth == 3
        bet = h[3:3] == "b" ? 3 : 0
        return findfirst(cards .== h[2:2])[1] + bet + sum(ni_stage[1:2])
    else
        return findfirst(cards .== h[1:1])[1] + sum(ni_stage[1:3])
    end
end

function terminal(h::String)
    depth = length(h)
    if depth > 4
        return true
    elseif depth < 4
        return false
    elseif h[3:4] == "pb"
        return false
    else
        return true
    end
end

function chance(h::String, a::String)
    1 / length(actions(h)) # assumes constant f_c(a|h)
end

function leafutility(h::String, i::Int64, πi::Float64, πo::Float64)
    ind = findfirst(leafnodes .== h)
    u = Uz[i][ind] * πo
    # println("h:$h, πo:$πo, u:$u ")
    return u
end

function matchregret(I::Int64, r::Vector{Vector{Float64}})
    σ = Vector{Float64}(undef, na[I])
    for a in 1:na[I]
        denom = sum(max(r[I][b], 0.0) for b in 1:na[I])
        if denom > 0.0
            σ[a] = max(r[I][a], 0.0) / denom
        else
            σ[a] = 1.0 / na[I]
        end
    end
    σ
end

function updateutility!(i::Int64, I::Int64, h::String, σ::Vector{Float64},
        u::Float64, πi::Float64, πo::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}})
    m = Vector{Float64}(undef, na[I])
    nextactions = actions(h)
    for a in 1:na[I]
        ha = string(h, nextactions[a])
        if player(h) == i
            πip = σ[a] * πi
            up = updatetree!(ha, i, πip, πo, r, s)
            m[a] = up
            u = u + σ[a] * up
        else
            πop = σ[a] * πo
            up = updatetree!(ha, i, πi, πop, r, s)
            u = u + up
        end
    end
    m, u
end

function updateregret!(i::Int64, I::Int64, h::String, σ::Vector{Float64},
        u::Float64, m::Vector{Float64}, πi::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}})
    if player(h) == i
        for a in 1:na[I]
            r[I][a] = r[I][a] + m[a] - u
            s[I][a] = s[I][a] + πi * σ[a]
        end
    end
end

function updatetree!(h::String, i::Int64, πi::Float64, πo::Float64,
    r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}})
    # πo is short for \vec{π}_{-i}
    (terminal(h)) && (return leafutility(h, i, πi, πo))
    if player(h) == 3
        as = actions(h)
        chance_weights = pweights([chance(h, a) for a in as])
        a = as[sample(chance_weights)]
        ha = string(h, a)
        return updatetree!(ha, i, πi, πo, r, s)
    end
    I = infoset(h)
    u = 0.0
    σ = matchregret(I, r)
    m, u = updateutility!(i, I, h, σ, u, πi, πo, r, s)
    updateregret!(i, I, h, σ, u, m, πi, r, s)
    u
end

function cfr(T::Int64)
    r = [zeros(na[i]) for i in 1:ni] # regret
    s = [zeros(na[i]) for i in 1:ni] # cumulative strategies
    u = 0.0
    u1 = Vector{Float64}(undef, T)
    println("Starting $T iterations...")
    for t in 1:T
        for i in N
            u = updatetree!("", i, 1.0, 1.0, r, s) # utility
            (i == 1) && (u1[t] = u)
        end
        t % 1_000 == 0 && print("\rIteration $t complete. Utility: $u")
    end
    denoms = [sum(s[i][:]) for i in 1:ni]
    σ = [s[i][:] ./ denoms[i] for i in 1:ni]
    u1, r, s, σ #-u for player 1's utility since we end on player 2
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
function LE(t; count = false)#this is Label Encoder
    n = 1
    le = length(t)
    target = ones(1, le)
    s = [v for v in Set(t)]
    c = zeros(1, length(s))
    for i in 1:le
        for ii in 1:length(s)
            if t[i] == s[ii]
                target[i] = ii
                c[ii] += 1
                break
            end
        end
    end
    d = Dict()
    for i in 1:length(s)
        d[s[i]] = i
    end
    if count
        return target, d, c
    else
        return target, d
    end
end
function CTOH(ts) #This is Convert to one hot
    try
        ts = broadcast(Int, ts)
    catch
        throw("Perhaps the array of objective variables you passed has a non-zero decimal point in the name of the class (as in 1.1). Please change it to an integer type.")
    end
    if 0 in ts
        ts .+= 1
    end
    c = maximum(ts)
    oh = zeros(size(ts)[2], c)
    i = 1
    for t in ts
        oh[i, t] = 1.0
        i += 1
    end
    return oh
end


text:
using CodecZlib
using FASTX
using TranscodingStreams
using FuzzySearch
using DelimitedFiles
#using TickTock

function write_record(stream::IO,record::FASTX.FASTQ.Record,rng::UnitRange{Int64},rc::Bool)
    # pline=[0x0a,0x2b,0x0a]
    if rc
        revcomplement_record!(record,rng)
    end
    if length(rng)>0
        seq_range=rng.+first(record.sequence).-1
        qual_range=rng.+first(record.quality).-1
        id_range=1:(first(record.sequence).-1)
        unsafe_write(stream,pointer(record.data,first(id_range)),length(id_range))
        unsafe_write(stream,pointer(record.data,first(seq_range)),length(seq_range))
        # the code below is slower
        # for i=id_range
        #     write(stream,record.data[i])
        # end
        # for i=seq_range
        #     write(stream,record.data[i])
        # end
        write(stream,0x0a,0x2b,0x0a)
        unsafe_write(stream,pointer(record.data,first(qual_range)),length(qual_range))
        # for i=qual_range
        #     write(stream,record.data[i])
        # end
        write(stream,0x0a)
    elseif length(rng)>0
        write(stream,record.data)
        write(stream,0x0a)
    else
        id_range=1:(first(record.sequence).-1)
        unsafe_write(stream,pointer(record.data,first(id_range)),length(id_range))
        write(stream,0x4e,0x0a,0x2b,0x0a,0x21)
        write(stream,0x0a);
    end
end

function revcomplement_record!(record::FASTX.FASTQ.Record,rng::UnitRange{Int64})
    seq_range=rng.+first(record.sequence).-1
    qual_range=rng.+first(record.quality).-1
    for i=seq_range
        if !(record.data[i]==0x4e)
            (record.data[i] & 0x02)>0 ? (record.data[i] ⊻= 0x04) : (record.data[i] ⊻= 0x15)
        end
    end
    reverse!(record.data,seq_range.start,seq_range.stop)
    reverse!(record.data,qual_range.start,qual_range.stop)
end

function suffix_mismatch_(data_short::Vector{UInt8},data_long::Vector{UInt8},rng_short::UnitRange{Int64},rng_long::UnitRange{Int64},maxerr::Int64,maxlu::Int64)
    nerr=0
    nmatch=0
    lu=0
    i=rng_short.start
    imax=rng_short.stop.+1
    j=rng_long.stop
    jmin=rng_short.start.-1
    while i<imax && j>jmin && (maxerr<0 || nerr<=maxerr) && (maxlu<0 || lu<maxlu)
        lu+=1
        x=data_short[i] | 0x20 #lower case it
        y=data_long[j] | 0x20
        if (x!==0x6e)
            (x & 0x02)>0 ? (x ⊻= 0x04) : (x ⊻= 0x15) # complement
            if x!==y
                nerr+=1
            else
                nmatch+=1
            end
        # else
            # nmatch+=1 # do not match N
        end
        j-=1
        i+=1
    end
    return nerr, nmatch
end

function is_suffix_(read1::FASTX.FASTQ.Record,read2::FASTX.FASTQ.Record,maxerr::Int64,minmatch::Int64,maxlu::Int64)
    #RNA1 is long, DNA1 is short
    #RNA2 is short, DNA2 is long
    (nerr,nmatch)=suffix_mismatch_(read1.data,read2.data,read1.sequence,read2.sequence,maxerr,maxlu)
    flag=(nerr<=maxerr && nmatch>=minmatch)
    return flag
end

function revcomp_stream(streamIN::IO,streamOUT::IO, verbose::Bool )
    reader = FASTQ.Reader(streamIN)
    record = FASTQ.Record()
    nread=0::Int64
    verbose_period=100000

    cumtime=1.0
    #tick()
    while !eof(reader)
        nread+=1
        if verbose && mod(nread,verbose_period)==0
            #cumtime+=tock()
            print(stderr,"Reverse complemented ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
            #tick()
        end
        read!(reader, record)
        revcomplement_record!(record,1:length(record.sequence))
        write(streamOUT,record.data)
        write(streamOUT,0x0a)
    end
    #cumtime+=tock()
    print(stderr,"DONE. Reverse complemented ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
    return nread
end

function revcompFQ(fqIN::String,fqOUT::String, verbose::Bool=true )
    if length(fqOUT)==0
        streamOUT=stdout
    else
        if endswith(fqOUT, ".gz")
            codec = GzipCompressor()
            streamOUT=TranscodingStream(codec, open(fqOUT,"w"))
        else
            codec = Noop()
            TranscodingStream(codec, open(fqOUT,"w"))
        end
    end
    streamIN=makestream_in(fqIN)
    revcomp_stream(streamIN,streamOUT,verbose)
    !(streamOUT==stdout) ? close(streamOUT) : true
    close(streamIN)
end

function readlength_summary_stream_(stream::IO,blen::Int64,maxreadl::Int64, verbose::Bool)
    M_F0=zeros(Int64,maxreadl,maxreadl)
    M_0R=zeros(Int64,maxreadl,maxreadl)
    M_F02=zeros(Int64,maxreadl,maxreadl)
    M_0R2=zeros(Int64,maxreadl,maxreadl)
    lR=0::Int64
    lD=0::Int64
    lR2=0::Int64
    lD2=0::Int64

    nread=0::Int64
    verbose_period=100000

    cumtime=1.0
    #tick()
    while !eof(stream)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,verbose_period)==0
            #cumtime+=tock()
            print(stderr,"Processed ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
            #tick()
        end

        s=readline(stream)
        x=parse.(split(s,'\t'))
        # print(x)
        if x[1]==1 && x[2]==2 #OF
            lR=x[6]-1
            lD=x[8]-(x[6]+blen)+1
            lR2=lR
            lD2=lD+x[7]
            M_F0[lR+1,lD+1]+=1
            M_F02[lR2+1,lD2+1]+=1
        elseif x[2]==1 && x[1]==2 #F0
            lR=x[5]-1
            lD=x[7]-(x[5]+blen)+1
            lR2=lR
            lD2=lD+x[8]
            M_F0[lR+1,lD+1]+=1
            M_F02[lR2+1,lD2+1]+=1
        elseif x[1]==3 && x[2]==1 #OR
            lD=x[5]-1
            lR=x[7]-(x[5]+blen)+1
            lD2=lD
            lR2=lR+x[8]
            M_0R[lR+1,lD+1]+=1
            M_0R2[lR2+1,lD2+1]+=1
        elseif x[2]==3 && x[1]==1 #OR
            lD=x[6]-1
            lR=x[8]-(x[6]+blen)+1
            lD2=lD
            lR2=lR+x[7]
            M_0R[lR+1,lD+1]+=1
            M_0R2[lR2+1,lD2+1]+=1
        end
    end
    return M_F02, M_0R2, M_F0, M_0R
end

function readlength_SE_summary_stream_(stream::IO,maxreadl::Int64, verbose::Bool)
    M=zeros(Int64,maxreadl,maxreadl)

    nread=0::Int64
    verbose_period=100000

    cumtime=1.0
    #tick()
    while !eof(stream)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,verbose_period)==0
            #cumtime+=tock()
            print(stderr,"Processed ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
            #tick()
        end

        s=readline(stream)
        x=parse.(split(s,' '))
        M[x[1]+1,x[2]+1]+=1
        # print(x)

    end
    return M
end

# function readlength_summary_stream_FAST_(stream::IO,blen::Int64,maxreadl::Int64, verbose::Bool)
#     M_F0=zeros(Int64,maxreadl,maxreadl)
#     M_0R=zeros(Int64,maxreadl,maxreadl)
#     M_F02=zeros(Int64,maxreadl,maxreadl)
#     M_0R2=zeros(Int64,maxreadl,maxreadl)
#     lR=0::Int64
#     lD=0::Int64
#     lR2=0::Int64
#     lD2=0::Int64
#
#     nread=0::Int64
#     verbose_period=100000
#
#     cumtime=0.0
#     tick()
#     while !eof(stream)
#         nread+=1
#         # print(nread,"\n")
#         if verbose && mod(nread,verbose_period)==0
#             #cumtime+=tock()
#             print(stderr,"Processed ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
#             tick()
#         end
#
#         x=read(stream,UInt8)
#         s=readline(stream)
#         x=parse.(split(s,'\t'))
#         # print(x)
#         if x[1]==1 && x[2]==2 #OF
#             lR=x[6]-1
#             lD=x[8]-(x[6]+blen)+1
#             lR2=lR
#             lD2=lD+x[7]
#             M_F0[lR+1,lD+1]+=1
#             M_F02[lR2+1,lD2+1]+=1
#         elseif x[2]==1 && x[1]==2 #F0
#             lR=x[5]-1
#             lD=x[7]-(x[5]+blen)+1
#             lR2=lR
#             lD2=lD+x[8]
#             M_F0[lR+1,lD+1]+=1
#             M_F02[lR2+1,lD2+1]+=1
#         elseif x[1]==3 && x[2]==1 #OR
#             lD=x[5]-1
#             lR=x[7]-(x[5]+blen)+1
#             lD2=lD
#             lR2=lR+x[8]
#             M_0R[lR+1,lD+1]+=1
#             M_0R2[lR2+1,lD2+1]+=1
#         elseif x[2]==3 && x[1]==1 #OR
#             lD=x[6]-1
#             lR=x[8]-(x[6]+blen)+1
#             lD2=lD
#             lR2=lR+x[7]
#             M_0R[lR+1,lD+1]+=1
#             M_0R2[lR2+1,lD2+1]+=1
#         end
#     end
#     return M_F02, M_0R2, M_F0, M_0R
# end

function readlength_summary(fname::String,fout::String,blen::Int64,maxreadl::Int64,verbose::Bool=true)
    in=makestream_in(fname)
    stream=open(fout,"w")
    M_F02, M_0R2, M_F0, M_0R = readlength_summary_stream_(in,blen,maxreadl,verbose)
    println(stream, "\n>>RNA/DNA length distribution")
    print(stream,"MF02"*"\t")
    Base.print(stream,M_F02,true)
    println(stream,"")
    print(stream,"M0R2"*"\t")
    Base.print(stream,M_0R2,true)
    println(stream,"")
    print(stream,"MF0"*"\t")
    Base.print(stream,M_F0,true)
    println(stream,"")
    print(stream,"M0R"*"\t")
    Base.print(stream,M_0R,true)
    println(stream,"")
    close(stream)
    close(in)
    return M_F02, M_0R2, M_F0, M_0R
end

function readlength_SE_summary(fname::String,fout::String,maxreadl::Int64,verbose::Bool=true)
    in=makestream_in(fname)
    stream=open(fout,"w")
    M = readlength_SE_summary_stream_(in,maxreadl,verbose)
    println(stream, "\n>>RNA/DNA length distribution")
    print(stream,"M"*"\t")
    Base.print(stream,M,true)
    println(stream,"")
    close(stream)
    close(in)
    return M
end

function splice_mates_stream_(streamR1::IO,streamR2::IO,streamD1::IO,streamD2::IO,stream_splicedR::IO,stream_splicedD::IO,stream_unsplicedR1::IO,stream_unsplicedR2::IO,stream_unsplicedD1::IO,stream_unsplicedD2::IO, err_match_lu_R::Tuple{Int64,Int64,Int64},err_match_lu_D::Tuple{Int64,Int64,Int64}, maxreadl::Int64, checkpairs::Bool, verbose::Bool )
    readerR1 = FASTQ.Reader(streamR1)
    recordR1 = FASTQ.Record()
    readerR2 = FASTQ.Reader(streamR2)
    recordR2 = FASTQ.Record()
    readerD1 = FASTQ.Reader(streamD1)
    recordD1 = FASTQ.Record()
    readerD2 = FASTQ.Reader(streamD2)
    recordD2 = FASTQ.Record()

    rd_len=zeros(Int64,maxreadl,maxreadl)
    rd_len_unspliced1=zeros(Int64,maxreadl,maxreadl)
    rd_len_unspliced2=zeros(Int64,maxreadl,maxreadl)

    nread=0::Int64
    nspliced=0::Int64
    verbose_period=100000

    cumtime=1.0
    #tick()
    while !eof(readerR1)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,verbose_period)==0
            #cumtime+=tock()
            print(stderr,"Spliced ", nread/(1000000.0), " million reads [", nread/cumtime/1000000*60, " million reads/min]\n")
            #tick()
        end
        read!(readerR1, recordR1)
        read!(readerR2, recordR2)
        read!(readerD1, recordD1)
        read!(readerD2, recordD2)

        if checkpairs
        end

        if (is_suffix_(recordR2,recordR1,err_match_lu_R[1],err_match_lu_R[2],err_match_lu_R[3]) && is_suffix_(recordD1,recordD2,err_match_lu_D[1],err_match_lu_D[2],err_match_lu_D[3]))
            nspliced+=1
            write(stream_splicedR,recordR1.data)
            write(stream_splicedR,0x0a)

            revcomplement_record!(recordD2,1:length(recordD2.sequence))
            write(stream_splicedD,recordD2.data)
            write(stream_splicedD,0x0a)

            rd_len[min(maxreadl+1,length(recordR1.sequence)+1),min(maxreadl+1,length(recordD2.sequence)+1)]+=1
        else
            write(stream_unsplicedR1,recordR1.data)
            write(stream_unsplicedR1,0x0a)
            write(stream_unsplicedR2,recordR2.data)
            write(stream_unsplicedR2,0x0a)
            write(stream_unsplicedD1,recordD1.data)
            write(stream_unsplicedD1,0x0a)
            write(stream_unsplicedD2,recordD2.data)
            write(stream_unsplicedD2,0x0a)
            rd_len_unspliced1[min(maxreadl+1,length(recordR1.sequence)+1),min(maxreadl+1,length(recordD1.sequence)+1)]+=1
            rd_len_unspliced2[min(maxreadl+1,length(recordR1.sequence)+1),min(maxreadl+1,length(recordD1.sequence)+1)]+=1

        end
    end
    return (nspliced,nread,rd_len,rd_len_unspliced1,rd_len_unspliced2)
end

function summarize_splicing_(read_counts, stream::IO)
    println(stream, ">>Total counts")
    println(stream,read_counts[2])
    println(stream, ">>Spliced counts")
    println(stream,read_counts[1])
    println(stream, "\n>>RNA/DNA length distribution")
    print(stream,"spliced"*"\t")
    Base.print(stream,read_counts[3])
    println(stream,"")
    print(stream,"unspliced_1"*"\t")
    Base.print(stream,read_counts[4])
    println(stream,"")
    print(stream,"unspliced_2"*"\t")
    Base.print(stream,read_counts[5])
    println(stream,"")

    # println("")
end

function splice_mates(fq_prefix::String,fq_spliced_prefix::String,fq_unspliced_prefix::String,err_match_lu_R_string::String,err_match_lu_D_string::String, maxreadl::Int64, checkpairs::Bool=false, verbose::Bool=true)

    streamSUMMARY=makestream_out(fq_spliced_prefix,".summary",false,".txt")

    streamR1=makestream_in(fq_prefix*".rna.1.fastq.gz")
    streamR2=makestream_in(fq_prefix*".rna.2.fastq.gz")
    streamD1=makestream_in(fq_prefix*".dna.1.fastq.gz")
    streamD2=makestream_in(fq_prefix*".dna.2.fastq.gz")
    stream_splicedR=makestream_out(fq_spliced_prefix,".rna",true,".fastq")
    stream_splicedD=makestream_out(fq_spliced_prefix,".dna",true,".fastq")
    stream_unsplicedR1=makestream_out(fq_unspliced_prefix,".rna.1",true,".fastq")
    stream_unsplicedR2=makestream_out(fq_unspliced_prefix,".rna.2",true,".fastq")
    stream_unsplicedD1=makestream_out(fq_unspliced_prefix,".dna.1",true,".fastq")
    stream_unsplicedD2=makestream_out(fq_unspliced_prefix,".dna.2",true,".fastq")

    err_match_lu_R=tuple(map(x->parse(Int64,x),split(err_match_lu_R_string,":"))...)
    err_match_lu_D=tuple(map(x->parse(Int64,x),split(err_match_lu_D_string,":"))...)

    read_counts = splice_mates_stream_(streamR1,streamR2,streamD1,streamD2,stream_splicedR,stream_splicedD,stream_unsplicedR1,stream_unsplicedR2,stream_unsplicedD1,stream_unsplicedD2,err_match_lu_R,err_match_lu_D,maxreadl,checkpairs, verbose)

    # println(streamSUMMARY, "total,spliced")
    # print(streamSUMMARY,nread)
    # print(streamSUMMARY,",")
    # println(streamSUMMARY, nspliced)

    summarize_splicing_(read_counts,streamSUMMARY)

    close(streamSUMMARY)
    close(stream_splicedR)
    close(stream_splicedD)
    close(stream_unsplicedR1)
    close(stream_unsplicedR2)
    close(stream_unsplicedD1)
    close(stream_unsplicedD2)
    close(streamR1)
    close(streamR2)
    close(streamD1)
    close(streamD2)

    return read_counts

end

function write_record_skipempty(stream::IO,record::FASTX.FASTQ.Record,rng::UnitRange{Int64},rc::Bool)
    # pline=[0x0a,0x2b,0x0a]
    if length(rng)>0
        seq_range=rng.+first(record.sequence).-1
        qual_range=rng.+first(record.quality).-1
        id_range=1:(first(record.sequence)-1)
        if rc
            revcomplement_record!(record,rng)
        end
        unsafe_write(stream,pointer(record.data,first(id_range)),length(id_range))
        unsafe_write(stream,pointer(record.data,first(seq_range)),length(seq_range))
        # the code below is slower
        # for i=id_range
        #     write(stream,record.data[i])
        # end
        # for i=seq_range
        #     write(stream,record.data[i])
        # end
        write(stream,0x0a,0x2b,0x0a)
        unsafe_write(stream,pointer(record.data,first(qual_range)),length(qual_range))
        # for i=qual_range
        #     write(stream,record.data[i])
        # end
        write(stream,0x0a)
    elseif length(rng)>0
        if rc
            revcomplement_record!(record,rng)
        end
        write(stream,record.data)
        write(stream,0x0a)
    end
end

function debridgeSE_stream_(bridgeF::String, bridgeR::String, streamX::IO, streamSTATS::Union{IO,Bool},streamsOUT::Union{Array{Array{IO,1},1},Bool}, verbose::Bool,maxreadlen::Int64,revcomp_bRreads::Bool)
    wF=Vector{UInt8}(bridgeF)
    wR=Vector{UInt8}(bridgeR)
    wF_length=length(wF)
    wR_length=length(wR)

    read_counts=zeros(Int64,4,1)
    bridge_pos_mean_X=zeros(Int64,4,1)
    maxl=maxreadlen-min(wF_length,wR_length).+1 #+1 because 1st element is length 0
    rd_len=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)]) #rna length / dna length when both, otherwise just read length

    readerX = FASTQ.Reader(streamX)
    recordX = FASTQ.Record()
    #flag1,flag2,N1,N2,position1,position2,l1,l2
    read_stats=zeros(Int64,1,4)
    writeDebridged=!(streamsOUT==false)
    writePositions=!(streamSTATS==false)
    if writeDebridged
        stream_array_DNA=streamsOUT[1]
        stream_array_RNA=streamsOUT[2]

        flip_read=convert(Vector{Bool},[0;0;1;0]) .& revcomp_bRreads
    end

    nread=0::Int64
    cumtime=1.0
    #tick()
    while !eof(readerX)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,100000)==0
            #cumtime+=tock()
            print(stderr,"Processed ", nread/1000000.0, " million reads [", round(nread/cumtime), " reads/s]\n")
            #tick()
        end
        read!(readerX, recordX)
        n1X=first(recordX.sequence)
        bridge_pos_X_F=brute_force_search_wildcard(wF,recordX.data,n1X,last(recordX.sequence))
        bridge_pos_X_R=brute_force_search_wildcard(wR,recordX.data,n1X,last(recordX.sequence))

        readlength_X=length(recordX.sequence)
        flagX=0x01+0x01*(bridge_pos_X_F[2]>0)+0x02*(bridge_pos_X_R[2]>0)
        bridge_pos_X_F[2]+bridge_pos_X_R[2]>1 ? flagX=0x04 : true

        read_counts[flagX]+=1
        bridge_pos_mean_X[flagX]+=bridge_pos_X_F[1]+bridge_pos_X_R[1]

        if writePositions
            read_stats[1]=flagX
            read_stats[2]=bridge_pos_X_F[2]+bridge_pos_X_R[2]
            read_stats[3]=max(bridge_pos_X_F[1],bridge_pos_X_R[1])
            read_stats[4]=readlength_X
            writedlm(streamSTATS,read_stats)
        end

        dnaX=1:readlength_X
        rnaX=0:-1
        if flagX==0x02 #unique bridgeF
            rnaX=1:(bridge_pos_X_F[1]-1)
            dnaX=(bridge_pos_X_F[1]+wF_length):readlength_X
            rd_len[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        elseif flagX==0x03 #unique bridgeR
            dnaX=1:(bridge_pos_X_R[1]-1)
            rnaX=(bridge_pos_X_R[1]+wR_length):readlength_X
            rd_len[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        else
            rd_len[flagX][min(maxreadlen,readlength_X)+1]+=1
        end

        if writeDebridged

            write_record(stream_array_DNA[flagX],recordX,dnaX,flip_read[flagX])
            write_record(stream_array_RNA[flagX],recordX,rnaX,flip_read[flagX])
        end

    end
    bridge_pos_mean_X_ok=bridge_pos_mean_X./read_counts
    return read_counts, bridge_pos_mean_X_ok, rd_len
end

function debridgeSE_stream_fuzzy_(bridgeF::String, bridgeR::String, streamX::IO, maxerr::Int64, streamSTATS::Union{IO,Bool}, streamsOUT::Union{Array{Array{IO,1},1},Bool}, verbose::Bool, maxreadlen::Int64,revcomp_bRreads::Bool) #note: using that method with maxerr=0 gives same result as debridgePE_stream_ but is 1.5x slower
    wF=Vector{UInt8}(bridgeF)
    wR=Vector{UInt8}(bridgeR)
    wF_length=length(wF)
    wR_length=length(wR)

    read_counts=zeros(Int64,4,1)
    bridge_pos_mean_X=zeros(Int64,4,1)
    maxl=maxreadlen-min(wF_length,wR_length)+1 #+1 because 1st element is length 0
    rd_len=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)]) #rna length / dna length when both, otherwise just read length

    readerX = FASTQ.Reader(streamX)
    recordX = FASTQ.Record()

    #flag1,flag2,N1,N2,position1,position2,l1,l2
    read_stats=zeros(Int64,1,6)
    writeDebridged=!(streamsOUT==false)
    writePositions=!(streamSTATS==false)
    if writeDebridged
        stream_array_DNA=streamsOUT[1]
        stream_array_RNA=streamsOUT[2]
        flip_read=convert(Vector{Bool},[0;0;1;0]) .& revcomp_bRreads
    end
    nread=0
    cumtime=1.0
    #tick()
    while !eof(readerX)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,100000)==0
           #cumtime+=tock()
            print(stderr,"Processed ", nread/1000000.0, " million reads [", round(nread/cumtime), " reads/s]\n")
            #tick()
        end
        read!(readerX, recordX)
        n1X=first(recordX.sequence)
        bridge_pos_X_F=brute_force_search_bestfuzzy(wF,recordX.data,n1X,last(recordX.sequence),maxerr)
        bridge_pos_X_R=brute_force_search_bestfuzzy(wR,recordX.data,n1X,last(recordX.sequence),maxerr)

        bridge_pos_X=bestof2searches(bridge_pos_X_F,bridge_pos_X_R)

        readlength_X=length(recordX.sequence)

        flagX=bridge_pos_X[1]

        read_counts[flagX]+=1
        bridge_pos_mean_X[flagX]+=bridge_pos_X_F[1]+bridge_pos_X_R[1]

        if writePositions
            read_stats[1]=flagX
            read_stats[2]=bridge_pos_X[3] #total including confusing ones
            read_stats[3]=bridge_pos_X[1]
            read_stats[4]=readlength_X
            read_stats[5]=bridge_pos_X[2]
            read_stats[6]=bridge_pos_X[4]
            writedlm(streamSTATS,read_stats)
        end

        dnaX=1:readlength_X
        rnaX=0:-1

        if flagX==0x02 #unique bridgeF
            rnaX=1:(bridge_pos_X_F[1]-1)
            dnaX=(bridge_pos_X_F[1]+wF_length):readlength_X
            rd_len[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        elseif flagX==0x03 #unique bridgeR
            dnaX=1:(bridge_pos_X_R[1]-1)
            rnaX=(bridge_pos_X_R[1]+wR_length):readlength_X
            rd_len[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        else
            rd_len[flagX][min(maxreadlen,readlength_X)+1]+=1
        end

        if writeDebridged
            write_record_skipempty(stream_array_DNA[flagX],recordX,dnaX,flip_read[flagX])
            write_record_skipempty(stream_array_RNA[flagX],recordX,rnaX,flip_read[flagX])
        end
    end
    bridge_pos_mean_X_ok=bridge_pos_mean_X./read_counts
    return read_counts, bridge_pos_mean_X_ok, rd_len
end

function debridgePE_stream_(bridgeF::String, bridgeR::String, streamX::IO, streamY::IO, streamSTATS::Union{IO,Bool},streamsOUT::Union{Array{Array{IO,2},1},Bool}, verbose::Bool, maxreadlen::Int64, revcomp_bRreads::Bool)
    wF=Vector{UInt8}(bridgeF)
    wR=Vector{UInt8}(bridgeR)
    wF_length=length(wF)
    wR_length=length(wR)

    read_counts=zeros(Int64,4,4)
    bridge_pos_mean_X=zeros(Int64,4,4)
    bridge_pos_mean_Y=zeros(Int64,4,4)
    maxl=maxreadlen-min(wF_length,wR_length)+1 #+1 because 1st element is length 0
    rd_lenX=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)])
    rd_lenY=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)])

    readerX = FASTQ.Reader(streamX)
    readerY = FASTQ.Reader(streamY)
    recordX = FASTQ.Record()
    recordY = FASTQ.Record()

    #flag1,flag2,N1,N2,position1,position2,l1,l2
    read_stats=zeros(Int64,1,8)
    writeDebridged=!(streamsOUT==false)
    writePositions=!(streamSTATS==false)
    if writeDebridged
        stream_array_DNA_PE1=streamsOUT[1]
        stream_array_DNA_PE2=streamsOUT[2]
        stream_array_RNA_PE1=streamsOUT[3]
        stream_array_RNA_PE2=streamsOUT[4]

        # RC in flip are :
        # dnaY in 0R
        # dnaX in R0
        # all in RM or RR or MR
        flip_DNA_PE1=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 1 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_DNA_PE2=convert(Array{Bool,2},[0 0 1 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_RNA_PE1=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_RNA_PE2=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads

    end

    nread=0::Int64
    cumtime=1.0
    #tick()
    while !eof(readerX)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,100000)==0
            #cumtime+=tock()
            print(stderr,"Processed ", nread/1000000.0, " million reads [", round(nread/cumtime), " reads/s]\n")
            #tick()
        end
        read!(readerX, recordX)
        read!(readerY, recordY)
        n1X=first(recordX.sequence)
        n1Y=first(recordY.sequence)
        bridge_pos_X_F=brute_force_search_wildcard(wF,recordX.data,n1X,last(recordX.sequence))
        bridge_pos_X_R=brute_force_search_wildcard(wR,recordX.data,n1X,last(recordX.sequence))
        bridge_pos_Y_F=brute_force_search_wildcard(wF,recordY.data,n1Y,last(recordY.sequence))
        bridge_pos_Y_R=brute_force_search_wildcard(wR,recordY.data,n1Y,last(recordY.sequence))

        readlength_X=length(recordX.sequence)
        readlength_Y=length(recordY.sequence)
        flagX=0x01+0x01*(bridge_pos_X_F[2]>0)+0x02*(bridge_pos_X_R[2]>0)
        flagY=0x01+0x01*(bridge_pos_Y_F[2]>0)+0x02*(bridge_pos_Y_R[2]>0)
        bridge_pos_X_F[2]+bridge_pos_X_R[2]>1 ? flagX=0x04 : true
        bridge_pos_Y_F[2]+bridge_pos_Y_R[2]>1 ? flagY=0x04 : true

        read_counts[flagX,flagY]+=1
        bridge_pos_mean_X[flagX,flagY]+=bridge_pos_X_F[1]+bridge_pos_X_R[1]
        bridge_pos_mean_Y[flagX,flagY]+=bridge_pos_Y_F[1]+bridge_pos_Y_R[1]

        if writePositions
            read_stats[1]=flagX
            read_stats[2]=flagY
            read_stats[3]=bridge_pos_X_F[2]+bridge_pos_X_R[2]
            read_stats[4]=bridge_pos_Y_F[2]+bridge_pos_Y_R[2]
            read_stats[5]=max(bridge_pos_X_F[1],bridge_pos_X_R[1])
            read_stats[6]=max(bridge_pos_Y_F[1],bridge_pos_Y_R[1])
            read_stats[7]=readlength_X
            read_stats[8]=readlength_Y
            writedlm(streamSTATS,read_stats)
        end

        dnaX=1:readlength_X
        dnaY=1:readlength_Y
        rnaX=0:-1
        rnaY=0:-1

        if flagX==0x02 #unique bridgeF
            rnaX=1:(bridge_pos_X_F[1]-1)
            dnaX=(bridge_pos_X_F[1]+wF_length):readlength_X
            rd_lenX[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        elseif flagX==0x03 #unique bridgeR
            dnaX=1:(bridge_pos_X_R[1]-1)
            rnaX=(bridge_pos_X_R[1]+wR_length):readlength_X
            rd_lenX[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        else
            rd_lenX[flagX][min(maxreadlen,readlength_X)+1]+=1
        end

        if flagY==0x02 #unique bridgeF
            rnaY=1:(bridge_pos_Y_F[1]-1)
            dnaY=(bridge_pos_Y_F[1]+wF_length):readlength_Y
            rd_lenY[flagY][min(maxl,length(rnaY)+1),min(maxl,length(dnaY)+1)]+=1
        elseif flagY==0x03 #unique bridgeR
            dnaY=1:(bridge_pos_Y_R[1]-1)
            rnaY=(bridge_pos_Y_R[1]+wR_length):readlength_Y
            rd_lenY[flagY][min(maxl,length(rnaY)+1),min(maxl,length(dnaY)+1)]+=1
        else
            rd_lenY[flagY][min(maxreadlen,readlength_Y)+1]+=1
        end

        if writeDebridged
            write_record(stream_array_DNA_PE1[flagX,flagY],recordX,dnaX,flip_DNA_PE1[flagX,flagY])
            write_record(stream_array_DNA_PE2[flagX,flagY],recordY,dnaY,flip_DNA_PE2[flagX,flagY])
            write_record(stream_array_RNA_PE1[flagX,flagY],recordX,rnaX,flip_RNA_PE1[flagX,flagY])
            write_record(stream_array_RNA_PE2[flagX,flagY],recordY,rnaY,flip_RNA_PE2[flagX,flagY])
        end
    end
    bridge_pos_mean_X_ok=bridge_pos_mean_X./read_counts
    bridge_pos_mean_Y_ok=bridge_pos_mean_Y./read_counts
    return read_counts, bridge_pos_mean_X_ok, bridge_pos_mean_Y_ok, rd_lenX, rd_lenY
end

function debridgePE_stream_fuzzy_(bridgeF::String, bridgeR::String, streamX::IO, streamY::IO, maxerr::Int64, streamSTATS::Union{IO,Bool}, streamsOUT::Union{Array{Array{IO,2},1},Bool}, verbose::Bool, maxreadlen::Int64, revcomp_bRreads::Bool) #note: using that method with maxerr=0 gives same result as debridgePE_stream_ but is 1.5x slower
    wF=Vector{UInt8}(bridgeF)
    wR=Vector{UInt8}(bridgeR)
    wF_length=length(wF)
    wR_length=length(wR)

    read_counts=zeros(Int64,4,4)
    bridge_pos_mean_X=zeros(Int64,4,4)
    bridge_pos_mean_Y=zeros(Int64,4,4)
    maxl=maxreadlen-min(wF_length,wR_length)+1 #+1 because 1st element is length 0
    rd_lenX=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)])
    rd_lenY=vcat([zeros(Int64,maxreadlen+1,1)],[zeros(Int64,maxl,maxl) for i in range(1,2, step=1)],[zeros(Int64,maxreadlen+1,1)])

    readerX = FASTQ.Reader(streamX)
    readerY = FASTQ.Reader(streamY)
    recordX = FASTQ.Record()
    recordY = FASTQ.Record()

    #flag1,flag2,N1,N2,position1,position2,l1,l2
    read_stats=zeros(Int64,1,12)
    writeDebridged=!(streamsOUT==false)
    writePositions=!(streamSTATS==false)
    if writeDebridged
        stream_array_DNA_PE1=streamsOUT[1]
        stream_array_DNA_PE2=streamsOUT[2]
        stream_array_RNA_PE1=streamsOUT[3]
        stream_array_RNA_PE2=streamsOUT[4]

        # RC in flip are :
        # dnaY in 0R
        # dnaX in R0
        # all in RM or RR or MR
        flip_DNA_PE1=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 1 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_DNA_PE2=convert(Array{Bool,2},[0 0 1 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_RNA_PE1=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
        flip_RNA_PE2=convert(Array{Bool,2},[0 0 0 0 ; 0 0 0 0 ; 0 0 1 1 ; 0 0 1 0]) .& revcomp_bRreads
    end
    nread=0
    cumtime=1.0
    #tick()
    while !eof(readerX)
        nread+=1
        # print(nread,"\n")
        if verbose && mod(nread,100000)==0
            #cumtime+=tock()
            print(stderr,"Processed ", nread/1000000.0, " million reads [", round(nread/cumtime), " reads/s]\n")
            #tick()
        end
        read!(readerX, recordX)
        read!(readerY, recordY)
        n1X=first(recordX.sequence)
        n1Y=first(recordY.sequence)
        bridge_pos_X_F=brute_force_search_bestfuzzy(wF,recordX.data,n1X,last(recordX.sequence),maxerr)
        bridge_pos_X_R=brute_force_search_bestfuzzy(wR,recordX.data,n1X,last(recordX.sequence),maxerr)
        bridge_pos_Y_F=brute_force_search_bestfuzzy(wF,recordY.data,n1Y,last(recordY.sequence),maxerr)
        bridge_pos_Y_R=brute_force_search_bestfuzzy(wR,recordY.data,n1Y,last(recordY.sequence),maxerr)

        bridge_pos_X=bestof2searches(bridge_pos_X_F,bridge_pos_X_R)
        bridge_pos_Y=bestof2searches(bridge_pos_Y_F,bridge_pos_Y_R)

        readlength_X=length(recordX.sequence)
        readlength_Y=length(recordY.sequence)

        flagX=bridge_pos_X[1]
        flagY=bridge_pos_Y[1]

        read_counts[flagX,flagY]+=1
        bridge_pos_mean_X[flagX,flagY]+=bridge_pos_X_F[1]+bridge_pos_X_R[1]
        bridge_pos_mean_Y[flagX,flagY]+=bridge_pos_Y_F[1]+bridge_pos_Y_R[1]

        if writePositions
            read_stats[1]=flagX
            read_stats[2]=flagY
            read_stats[3]=bridge_pos_X[3] #total including confusing ones
            read_stats[4]=bridge_pos_Y[3]
            read_stats[5]=bridge_pos_X[1]
            read_stats[6]=bridge_pos_Y[1]
            read_stats[7]=readlength_X
            read_stats[8]=readlength_Y
            read_stats[9]=bridge_pos_X[2]
            read_stats[10]=bridge_pos_Y[2]
            read_stats[11]=bridge_pos_X[4]
            read_stats[12]=bridge_pos_Y[4]
            writedlm(streamSTATS,read_stats)
        end

        dnaX=1:readlength_X
        dnaY=1:readlength_Y
        rnaX=0:-1
        rnaY=0:-1

        if flagX==0x02 #unique bridgeF
            rnaX=1:(bridge_pos_X_F[1]-1)
            dnaX=(bridge_pos_X_F[1]+wF_length):readlength_X
            rd_lenX[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        elseif flagX==0x03 #unique bridgeR
            dnaX=1:(bridge_pos_X_R[1]-1)
            rnaX=(bridge_pos_X_R[1]+wR_length):readlength_X
            rd_lenX[flagX][min(maxl,length(rnaX)+1),min(maxl,length(dnaX)+1)]+=1
        else
            rd_lenX[flagX][min(maxreadlen,readlength_X)+1]+=1
        end

        if flagY==0x02 #unique bridgeF
            rnaY=1:(bridge_pos_Y_F[1]-1)
            dnaY=(bridge_pos_Y_F[1]+wF_length):readlength_Y
            rd_lenY[flagY][min(maxl,length(rnaY)+1),min(maxl,length(dnaY)+1)]+=1
        elseif flagY==0x03 #unique bridgeR
            dnaY=1:(bridge_pos_Y_R[1]-1)
            rnaY=(bridge_pos_Y_R[1]+wR_length):readlength_Y
            rd_lenY[flagY][min(maxl,length(rnaY)+1),min(maxl,length(dnaY)+1)]+=1
        else
            rd_lenY[flagY][min(maxreadl,readlength_Y)+1]+=1
        end

        if writeDebridged
            # write_record(stream_array_DNA_PE1[flagX,flagY],recordX,dnaX)
            # write_record(stream_array_DNA_PE2[flagX,flagY],recordY,dnaY)
            # write_record(stream_array_RNA_PE1[flagX,flagY],recordX,rnaX)
            # write_record(stream_array_RNA_PE2[flagX,flagY],recordY,rnaY)
            write_record(stream_array_DNA_PE1[flagX,flagY],recordX,dnaX,flip_DNA_PE1[flagX,flagY])
            write_record(stream_array_DNA_PE2[flagX,flagY],recordY,dnaY,flip_DNA_PE2[flagX,flagY])
            write_record(stream_array_RNA_PE1[flagX,flagY],recordX,rnaX,flip_RNA_PE1[flagX,flagY])
            write_record(stream_array_RNA_PE2[flagX,flagY],recordY,rnaY,flip_RNA_PE2[flagX,flagY])
        end
    end
    bridge_pos_mean_X_ok=bridge_pos_mean_X./read_counts
    bridge_pos_mean_Y_ok=bridge_pos_mean_Y./read_counts
    return read_counts, bridge_pos_mean_X_ok, bridge_pos_mean_Y_ok, rd_lenX, rd_lenY
end

function makestream_out(fileprefix::String,filename::String, compress::Bool=true, ext::String=".fastq")
    if isempty(fileprefix)
        return stdout
    else
        if compress
            codec=GzipCompressor()
            return TranscodingStream(codec, open(fileprefix * filename * ext * ".gz","w"))
        else
            codec=Noop()
            return TranscodingStream(codec, open(fileprefix * filename * ext,"w"))
        end
    end
end

function makestream_in(filepath::String)
    if endswith(filepath, ".gz")
        codec = GzipDecompressor()
        return TranscodingStream(codec, open(filepath))
    else
        codec = Noop()
        return TranscodingStream(codec, open(filepath))
    end
end

function summarize_debridgingPE_(read_counts, stream::IO)
    println(stream, ">>Total counts")
    println(stream, sum(read_counts[1]))
    println(stream,">>Count matrix")
    Base.print(stream,read_counts[1])
    println(stream, "\n>>Count %")
    Base.print(stream,read_counts[1]*100.0/sum(read_counts[1]))
    println(stream, "\n>>Bridge mean position read 1")
    Base.print(stream,read_counts[2])
    println(stream, "\n>>Bridge mean position read 2")
    Base.print(stream,read_counts[3])
    println(stream, "\n>>RNA/DNA length distribution read 1/2")
    bridgecodes=["0","F","R","M"]
    for i=1:4
        print(stream,"bridge_" * bridgecodes[i] * "_1\t")
        Base.print(stream,read_counts[4][i])
        print(stream,"\nbridge_" * bridgecodes[i] * "_2\t")
        Base.print(stream,read_counts[5][i])
        println(stream,"")
    end
    # println("")
end

function summarize_debridgingSE_(read_counts, stream::IO)
    println(stream, ">>Total counts")
    println(stream,sum(read_counts[1]))
    println(stream,">>Count matrix")
    Base.print(stream,read_counts[1])
    println(stream, "\n>>Count %")
    Base.print(stream,read_counts[1]*100.0/sum(read_counts[1]))
    println(stream, "\n>>Bridge mean position")
    Base.print(stream,read_counts[2])
    println(stream, "\n>>RNA/DNA length distribution")
    bridgecodes=["0","F","R","M"]
    for i=1:4
        print(stream,"bridge_" * bridgecodes[i] * "\t")
        Base.print(stream,read_counts[3][i])
        println(stream,"")
    end
    # println("")
end

function debridgePE(bridgeF::String, bridgeR::String, fqX::String, fqY::String, fqout_prefix::String, maxerr::Int64=0, writePositions::Int64=2, writeDebridged::Int64=0,  verbose::Bool=true, maxreadl::Int=151, rc::Bool=true)

    streamX=makestream_in(fqX)
    streamY=makestream_in(fqY)

    streamSUMMARY=makestream_out(fqout_prefix,"summary.PE",false,".txt")
    if rc
        println(streamSUMMARY, "Running debridge PE with reverse complement mode ON")
    end
    if writePositions>0
        streamSTATS=makestream_out(fqout_prefix,"positions.PE",writePositions>1,".txt") #2 stands for PE
    else
        streamSTATS=false
    end
    if writeDebridged>0
        compressOut=writeDebridged>1
        #00 reads

        # RC in flip are :
        # dnaY in 0R
        # dnaX in R0
        # all in RM or RR
        NULL1=makestream_out(fqout_prefix,"00.xxx.1",compressOut)
        NULL2=makestream_out(fqout_prefix,"00.xxx.2",compressOut)
        #F0 reads
        F0rna1=makestream_out(fqout_prefix,"F0.rna.1",compressOut)#sense
        F0dna1=makestream_out(fqout_prefix,"F0.dna.1",compressOut)
        F0dna2=makestream_out(fqout_prefix,"F0.dna.2",compressOut)
        #R0 reads (read2 more likely RNA?)

        if rc
            R0rna1=makestream_out(fqout_prefix, "0R.rna.2",compressOut) #comes from the R read, which is mate 2 in flip mode. This RNA is antisense
            R0dna1=makestream_out(fqout_prefix, "0R.dna",compressOut) #THIS WILL GET RC
            R0rna2=makestream_out(fqout_prefix, "0R.rna.1",compressOut) # comes from the zero read which is .1 in flip version. This RNA is sense
        else
            R0rna1=makestream_out(fqout_prefix, "R0.rna.1",compressOut) # comes from the R read, which is mate 1 in original mode. This RNA is antisense.
            R0dna1=makestream_out(fqout_prefix, "R0.dna.1",compressOut)
            R0rna2=makestream_out(fqout_prefix, "R0.rna.2",compressOut) # comes from the zero read which is .2 in the initial no flip version. This RNA is sense
        end

        #FM reads
        FMrna1=makestream_out(fqout_prefix, "FM.rna.1",compressOut) #sense
        FMdna1=makestream_out(fqout_prefix, "FM.dna.1",compressOut)
        FMdna2=makestream_out(fqout_prefix, "FM.xxx.2",compressOut)
        #RM reads
        RMrna1=makestream_out(fqout_prefix, "RM.rna.1",compressOut) #antisense, or gets RC in flip mode so sense
        RMdna1=makestream_out(fqout_prefix, "RM.dna.1",compressOut) #antisense, or gets RC in flip mode so sense
        RMdna2=makestream_out(fqout_prefix, "RM.xxx.2",compressOut) # gets RC in flip mode
        #M0 reads
        M01=makestream_out(fqout_prefix, "M0.xxx.1",compressOut)
        M02=makestream_out(fqout_prefix, "M0.xxx.2",compressOut)
        #MM reads
        MM1=makestream_out(fqout_prefix, "MM.xxx.1",compressOut)
        MM2=makestream_out(fqout_prefix, "MM.xxx.2",compressOut)
        #FF reads
        FFrna1=makestream_out(fqout_prefix, "FF.rna.1",compressOut)#sense
        FFdna1=makestream_out(fqout_prefix, "FF.dna.1",compressOut)
        FFrna2=makestream_out(fqout_prefix, "FF.rna.2",compressOut)#sense
        FFdna2=makestream_out(fqout_prefix, "FF.dna.2",compressOut)
        #RR reads
        RRrna1=makestream_out(fqout_prefix, "RR.rna.1",compressOut) #antisense, get Rc in flip
        RRdna1=makestream_out(fqout_prefix, "RR.dna.1",compressOut) # get RC in flip
        RRrna2=makestream_out(fqout_prefix, "RR.rna.2",compressOut)#antisense #gets RC in flip
        RRdna2=makestream_out(fqout_prefix, "RR.dna.2",compressOut)#gets RC in flip
        #FR reads
        FRrna1=makestream_out(fqout_prefix, "FR.rna.1",compressOut) #sense
        FRdna1=makestream_out(fqout_prefix, "FR.dna.1",compressOut)
        FRrna2=makestream_out(fqout_prefix, "FR.rna.2",compressOut) #antisense
        FRdna2=makestream_out(fqout_prefix, "FR.dna.2",compressOut)

        NULLstream=stdout #

        # pritories for reversing reads: F, R, M, 0
        #[00, 0F*, 0R*&, 0M*; F0, FF, FR, FM; R0&, RF*, RR, RM; M0, MF*, MR*, MM] * we reverse the reads, & 2nd read (after reverse) contains no DNA

        #if rc
        #    R0rna1=makestream_out(fqout_prefix, "0R.rna.2",compressOut) #comes from the R read, which is mate 2 in flip mode. This RNA is antisense
    #        R0dna1=makestream_out(fqout_prefix, "0R.dna",compressOut) #THIS WILL GET RC
#            R0rna2=makestream_out(fqout_prefix, "0R.rna.1",compressOut) # comes from the zero read which is .1 in flip version. This RNA is sense

        stream_array_DNA_PE1=[NULL1 F0dna2 R0rna2 M02; F0dna1 FFdna1 FRdna1 FMdna1; R0dna1 FRdna2 RRdna1 RMdna1; M01 FMdna2 RMdna2 MM1]
        stream_array_DNA_PE2=[NULL2 F0dna1 R0dna1 M01; F0dna2 FFdna2 FRdna2 FMdna2; R0rna2 FRdna1 RRdna2 RMdna2; M02 FMdna1 RMdna1 MM2]
        #[00&&, 0F*&, 0R*, 0M*&&; F0&, FF, FR, FM&; R0, RF*, RR, RM&; M0&&, MF*&, MR*&, MM&&] * we reverse the reads, & 2nd read contains no RNA
        stream_array_RNA_PE1=[NULLstream NULLstream NULLstream NULLstream; F0rna1 FFrna1 FRrna1 FMrna1; R0rna1 FRrna2 RRrna1 RMrna1; NULLstream NULLstream NULLstream NULLstream]
        stream_array_RNA_PE2=[NULLstream F0rna1 R0rna1 NULLstream; NULLstream FFrna2 FRrna2 NULLstream; NULLstream FRrna1 RRrna2 NULLstream; NULLstream FMrna1 RMrna1 NULLstream]

        streamsOUT=[stream_array_DNA_PE1,stream_array_DNA_PE2,stream_array_RNA_PE1,stream_array_RNA_PE2]
    else
        streamsOUT=false
    end
    if maxerr>0
        read_counts=debridgePE_stream_fuzzy_(bridgeF, bridgeR, streamX, streamY, maxerr,streamSTATS,streamsOUT, verbose, maxreadl, rc)
    else
        read_counts=debridgePE_stream_(bridgeF, bridgeR, streamX, streamY,streamSTATS,streamsOUT, verbose, maxreadl, rc)
    end
    if writeDebridged>0
        for i=1:4
            for j=1:4
                !(stream_array_DNA_PE1[i,j]==stdout) ? close(stream_array_DNA_PE1[i,j]) : true
                !(stream_array_DNA_PE2[i,j]==stdout) ? close(stream_array_DNA_PE2[i,j]) : true
                !(stream_array_RNA_PE1[i,j]==stdout) ? close(stream_array_RNA_PE1[i,j]) : true
                !(stream_array_RNA_PE1[i,j]==stdout) ? close(stream_array_RNA_PE1[i,j]) : true
            end
        end
    end
    if writePositions>0
        !(streamSTATS==stdout) ? close(streamSTATS) : true
    end

    close(streamX)
    close(streamY)

    summarize_debridgingPE_(read_counts,streamSUMMARY)
    !(streamSUMMARY==stdout) ? close(streamSUMMARY) : true
    return read_counts

end

function debridgeSE(bridgeF::String, bridgeR::String, fqX::String, fqout_prefix::String, maxerr::Int64=0, writePositions::Int64=1, writeDebridged::Int64=0,  verbose::Bool=true, maxreadl::Int=151, rc::Bool=true)
    streamX=makestream_in(fqX)
    streamSUMMARY=makestream_out(fqout_prefix,"summary.SE",false,".txt")

    if rc
        println(streamSUMMARY, "Running debridge SE with reverse complement mode ON")
    end

    if writePositions>0
        streamSTATS=makestream_out(fqout_prefix,"positions.SE",writePositions>1,".txt")
    else
        streamSTATS=false
    end

    if writeDebridged>0
        compressOut=writeDebridged>1
        #0 reads
        NULL=makestream_out(fqout_prefix,"0.xxx",compressOut)
        #F reads
        Frna=makestream_out(fqout_prefix,"F.rna",compressOut)#sense
        Fdna=makestream_out(fqout_prefix,"F.dna",compressOut)

        #R reads (read2 more likely RNA?)
        Rrna=makestream_out(fqout_prefix, "R.rna",compressOut) # RNA is antisense in nomrmal or sense in FLIP mode
        Rdna=makestream_out(fqout_prefix, "R.dna",compressOut)

        #M reads
        M=makestream_out(fqout_prefix, "M.xxx",compressOut)

        NULLstream=stdout #

        # pritories for reversing reads: F, R, M, 0
        #[00, 0F*, 0R*&, 0M*; F0, FF, FR, FM; R0&, RF*, RR, RM; M0, MF*, MR*, MM] * we reverse the reads, & 2nd read (after reverse) contains no DNA
        stream_array_DNA=[NULL, Fdna, Rdna, M]
        stream_array_RNA=[NULLstream, Frna, Rrna, NULLstream]
        streamsOUT=[stream_array_DNA,stream_array_RNA]
    else
        streamsOUT=false
    end

    if maxerr>0
        read_counts=debridgeSE_stream_fuzzy_(bridgeF, bridgeR, streamX, maxerr,streamSTATS,streamsOUT, verbose, maxreadl, rc)
    else
        read_counts=debridgeSE_stream_(bridgeF, bridgeR, streamX ,streamSTATS,streamsOUT, verbose, maxreadl, rc)
    end
    if writeDebridged>0
        for i=1:4
            !(stream_array_DNA[i]==stdout) ? close(stream_array_DNA[i]) : true
            !(stream_array_RNA[i]==stdout) ? close(stream_array_RNA[i]) : true
        end
    end
    if writePositions>0
        !(streamSTATS==stdout) ? close(streamSTATS) : true
    end

    close(streamX)
    summarize_debridgingSE_(read_counts,streamSUMMARY)
    !(streamSUMMARY==stdout) ? close(streamSUMMARY) : true
    return read_counts

end

#

# bF="AAAAAAACCGGCGTCCAAG"
# bR="CTTGGACGCCGGTTTTTTT"
# fqX="s1head_R1.fastq.gz"
# fqY="s1head_R2.fastq.gz"
# fqOutPrefix="out_"

max_stars_repo_name: straightlab/chartools-dev

max_issues_repo_name: straightlab/chartools-dev

max_forks_repo_name: straightlab/chartools-dev


text:
using Pkg
using PyCall

if lowercase(get(ENV, "CI", "false")) == "true"

        ENV["PYTHON"] = ""
        Pkg.build("PyCall")

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1ReplicaSetStatus <: SwaggerModel
    availableReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: availableReplicas
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1ReplicaSetCondition} } # spec name: conditions
    fullyLabeledReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: fullyLabeledReplicas
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration
    readyReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: readyReplicas
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas

    function IoK8sApiExtensionsV1beta1ReplicaSetStatus(;availableReplicas=nothing, conditions=nothing, fullyLabeledReplicas=nothing, observedGeneration=nothing, readyReplicas=nothing, replicas=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("availableReplicas"), availableReplicas)
        setfield!(o, Symbol("availableReplicas"), availableReplicas)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("fullyLabeledReplicas"), fullyLabeledReplicas)
        setfield!(o, Symbol("fullyLabeledReplicas"), fullyLabeledReplicas)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("readyReplicas"), readyReplicas)
        setfield!(o, Symbol("readyReplicas"), readyReplicas)
        validate_property(IoK8sApiExtensionsV1beta1ReplicaSetStatus, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        o
    end
end # type IoK8sApiExtensionsV1beta1ReplicaSetStatus

const _property_map_IoK8sApiExtensionsV1beta1ReplicaSetStatus = Dict{Symbol,Symbol}(Symbol("availableReplicas")=>Symbol("availableReplicas"), Symbol("conditions")=>Symbol("conditions"), Symbol("fullyLabeledReplicas")=>Symbol("fullyLabeledReplicas"), Symbol("observedGeneration")=>Symbol("observedGeneration"), Symbol("readyReplicas")=>Symbol("readyReplicas"), Symbol("replicas")=>Symbol("replicas"))
const _property_types_IoK8sApiExtensionsV1beta1ReplicaSetStatus = Dict{Symbol,String}(Symbol("availableReplicas")=>"Int32", Symbol("conditions")=>"Vector{IoK8sApiExtensionsV1beta1ReplicaSetCondition}", Symbol("fullyLabeledReplicas")=>"Int32", Symbol("observedGeneration")=>"Int64", Symbol("readyReplicas")=>"Int32", Symbol("replicas")=>"Int32")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1ReplicaSetStatus }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1ReplicaSetStatus))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1ReplicaSetStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1ReplicaSetStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1ReplicaSetStatus }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1ReplicaSetStatus[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1ReplicaSetStatus)
    (getproperty(o, Symbol("replicas")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1ReplicaSetStatus }, name::Symbol, val)
end


text:
immutable SMC <: InferenceAlgorithm
  n_particles :: Int
  resampler :: Function
  resampler_threshold :: Float64
  SMC(n) = new(n, resampleSystematic, 0.5)
end

assume(spl :: Sampler{SMC}, distr :: Distribution)  = rand(current_trace(), distr)

## wrapper for smc: run the sampler, collect results.
function Base.run(spl::Sampler{SMC})

  spl.particles = ParticleContainer{TraceC}(spl.model)
  push!(spl.particles, spl.alg.n_particles)

  while consume(spl.particles) != Val{:done}
    ess = effectiveSampleSize(spl.particles)
    if ess <= spl.alg.resampler_threshold * length(spl.particles)
      resample!(spl.particles)
    end
  end

  res = Chain(spl.particles)

end


text:
using SparseArrays
using LinearAlgebra
using MAT
import Base.exp;
import Base.sin;
import Base.cos;

struct HelmholtzProblem2D
    sigma
    Aders
    Afun
    b
end

struct HelmholtzProblem2D_nlterms
    sigma
    Afun
    fv
    fders
    Cv
    b
end

function compute_Mlincomb(prob::HelmholtzProblem2D,sigma,Y)
    if (prob.sigma != sigma)
        @show sigma
        @show prob.sigma
        error("Wrong sigma");
    end
    z=zeros(eltype(Y),size(Y,1));
    for i=1:size(Y,2)
        z += prob.Aders[i]*Y[:,i];
    end
    return z;
end
function compute_Mder(prob::HelmholtzProblem2D,mu)
    return prob.Afun(mu);
end

"""

    prob=create_helmholtz(n,a,b,c,N,k_eval,beta_eval,h_eval)

Creates an object representing the discretization of

```u''(x)+f1(mu)*(1+mu*k(x))^2*u+f2(mu)*beta(x)*u=h(x)```

using matrices from a finite element discretization.

"""
function create_helmholtz2D(N,A0,A1,A2,A3,A4,b)

    n = length(b);

    ders=Vector{typeof(A0)}(undef,N);
    ders[1] = A0;
    ders[2] = A1 + A4;

    coeff=1;
    for i=4:2:N
       coeff = coeff*(-1);
       ders[i] = coeff*A4;
    end
    ders[4] += 6*A3;
    for i=3:2:N
        ders[i] = spzeros(n,n);
    end
    ders[3] = 4*A2;

    # Create the A-function
    Afun = mu-> (A0 + mu*A1 + (2*mu^2)*A2 + (mu^3)*A3 + sin(mu)*A4);

    rhs=b;  # right-hand side

    myproblem=HelmholtzProblem2D(0,ders,Afun,rhs);

end

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
#####
##### parsing utilities
#####

function Base.tryparse(::Type{PatientID}, raw::AbstractString)
    metadata = split(raw, ' '; keepempty=false)
    length(metadata) == 4 || return nothing
    code_raw, sex_raw, dob_raw, name_raw = metadata
    length(sex_raw) == 1 || return nothing
    code = edf_unknown(code_raw)
    sex = edf_unknown(first, sex_raw)
    dob = edf_unknown(parse_date, dob_raw)
    dob === nothing && return nothing
    name = edf_unknown(name_raw)
    return PatientID(code, sex, dob, name)
end

function Base.tryparse(::Type{RecordingID}, raw::AbstractString)
    startswith(raw, "Startdate") || return nothing
    metadata = split(chop(raw; head=9, tail=0), ' '; keepempty=false)
    length(metadata) == 4 || return nothing
    start_raw, admin_raw, tech_raw, equip_raw = metadata
    startdate = edf_unknown(parse_date, start_raw)
    startdate === nothing && return nothing
    admincode = edf_unknown(admin_raw)
    technician = edf_unknown(tech_raw)
    equipment = edf_unknown(equip_raw)
    return RecordingID(startdate, admincode, technician, equipment)
end

# EDF does not specify what to do if the month or day is not in a valid range
# so if it is not we snap month or day to "01" and try again
function parse_header_date(date_str::AbstractString)
    m = match(r"^(\d{2})\.(\d{2})\.(\d{2}) (\d{2})\.(\d{2})\.(\d{2})$", date_str)
    if m === nothing
        throw(ArgumentError("Malformed date string: expected 'dd.mm.yy HH.MM.SS', " *
                            "got '$date_str'"))
    end
    day, month, year, hour, minute, second = parse.(Int, m.captures)
    if year <= 84
        year += 2000
    else
        year += 1900
    end

    # FIXME: EDF "avoids" the Y2K problem by punting it to the year 2084, after which
    # we ignore the above entirely and use `recording_id.startdate`. We could add a
    # check here on `year(today())`, but that will be dead code for the next 60+ years.

    month = clamp(month, 1, 12)
    day = clamp(day, 1, daysinmonth(year, month))
    hour = clamp(hour, 0, 23)
    minute = clamp(minute, 0, 59)
    second = clamp(second, 0, 59)
    return DateTime(year, month, day, hour, minute, second)
end

parse_float(raw::AbstractString) = something(tryparse(Float32, raw), NaN32)

parse_date(raw::AbstractString) = tryparse(Date, raw, dateformat"d-u-y")

"""
    edf_unknown([f,] field::String)

Check whether the given text is simply "X", which per the EDF+ specification means that the
value is unknown, not applicable, or must be made anonymous, and if so return `missing`,
otherwise return the result of applying `f` to the given text. If unspecified, `f` defaults
to `identity`.
"""
edf_unknown(f, field::AbstractString) = field == "X" ? missing : f(field)
edf_unknown(field::AbstractString) = edf_unknown(identity, field)

#####
##### header reading utilities
#####

function read_file_header(io::IO)
    version = strip(String(Base.read(io, 8)))

    patient_id_raw = strip(String(Base.read(io, 80)))
    patient_id = something(tryparse(PatientID, patient_id_raw), String(patient_id_raw))

    recording_id_raw = strip(String(Base.read(io, 80)))
    recording_id = something(tryparse(RecordingID, recording_id_raw), String(recording_id_raw))

    start_raw = Base.read(io, 8)
    push!(start_raw, UInt8(' '))
    append!(start_raw, Base.read(io, 8))  # Add the time
    # Parsing the date per the given format will validate EDF+ item 2
    start = parse_header_date(String(start_raw))

    # NOTE: These 8 bytes are supposed to define the byte count of the header,
    # which in reality is trivially computable from constants defined by the
    # specification + directly available information in the header already. I'm
    # not sure why the EDF standard requires that it be written out at all; AFAICT
    # it only serves as a potential bug source for readers/writers that might write
    # the incorrect value here. Since we don't actually use this value anywhere in
    # our read/write process, we skip it here.
    skip(io, 8)

    reserved = String(Base.read(io, 44))
    is_contiguous = !startswith(reserved, "EDF+D")
    record_count = parse(Int, String(Base.read(io, 8)))
    seconds_per_record = parse(Float64, String(Base.read(io, 8)))
    signal_count = parse(Int, String(Base.read(io, 4)))
    return FileHeader(version, patient_id, recording_id, start, is_contiguous,
                      record_count, seconds_per_record), signal_count
end

function read_signal_headers(io::IO, signal_count)
    fields = [String(Base.read(io, size)) for signal in 1:signal_count, (_, size) in SIGNAL_HEADER_FIELDS]
    signal_headers = [SignalHeader(strip(fields[i, 1]), strip(fields[i, 2]),
                                   strip(fields[i, 3]), parse_float(fields[i, 4]),
                                   parse_float(fields[i, 5]), parse_float(fields[i, 6]),
                                   parse_float(fields[i, 7]), strip(fields[i, 8]),
                                   parse(Int16, fields[i, 9])) for i in 1:size(fields, 1)]
    skip(io, 32 * signal_count) # reserved
    return signal_headers
end

#####
##### signal reading utilities
#####

# NOTE: The fast-path in `Base.read!` that uses `unsafe_read` will read too much when
# the element type is `Int24`, since it will try to include the alignment padding for
# each value read and will thus read too much. To get around this, we'll fall back to
# a naive implementation when the size of the element type doesn't match its aligned
# size. (See also `write_from`)
function read_to!(io::IO, x::AbstractArray{T}) where {T}
    if sizeof(T) == Base.aligned_sizeof(T)
        Base.read!(io, x)
    else
        @inbounds for i in eachindex(x)
            x[i] = Base.read(io, T)
        end
    end
    return x
end

function read_signals!(file::File)
    for record_index in 1:file.header.record_count, signal in file.signals
        read_signal_record!(file, signal, record_index)
    end
    return nothing
end

function read_signal_record!(file::File, signal::Signal, record_index::Int)
    if isempty(signal.samples)
        resize!(signal.samples, file.header.record_count * signal.header.samples_per_record)
    end
    record_start = 1 + (record_index - 1) * signal.header.samples_per_record
    record_stop = record_index * signal.header.samples_per_record
    read_to!(file.io, view(signal.samples, record_start:record_stop))
    return nothing
end

function read_signal_record!(file::File, signal::AnnotationsSignal, record_index::Int)
    bytes_per_sample = sizeof(sample_type(file))
    io_for_record = IOBuffer(Base.read(file.io, bytes_per_sample * signal.samples_per_record))
    tals_for_record = TimestampedAnnotationList[]
    while !eof(io_for_record) && Base.peek(io_for_record) != 0x00
        push!(tals_for_record, read_tal(io_for_record))
    end
    push!(signal.records, tals_for_record)
    return nothing
end

function read_tal(io::IO)
    sign = read_tal_onset_sign(io)
    bytes = readuntil(io, 0x14)
    timestamp = split(String(bytes), '\x15'; keepempty=false)
    onset_in_seconds = flipsign(parse(Float64, timestamp[1]), sign)
    duration_in_seconds = length(timestamp) == 2 ? parse(Float64, timestamp[2]) : nothing
    annotations = convert(Vector{String}, split(String(readuntil(io, 0x00)), '\x14'; keepempty=true))
    isempty(last(annotations)) && pop!(annotations)
    return TimestampedAnnotationList(onset_in_seconds, duration_in_seconds, annotations)
end

function read_tal_onset_sign(io::IO)
    sign = Base.read(io, UInt8)
    sign === 0x2b && return 1
    sign === 0x2d && return -1
    error("starting byte of a TAL must be '+' or '-'; found $sign")
end

#####
##### API functions
#####

"""
    EDF.File(io::IO)

Return an `EDF.File` instance that wraps the given `io`, as well as EDF-formatted
file, signal, and annotation headers that are read from `io`. This constructor
only reads headers, not the subsequent sample data; to read the subsequent sample
data from `io` into the returned `EDF.File`, call `EDF.read!(file)`.
"""
function File(io::IO)
    file_header, signal_count = read_file_header(io)
    T = sample_type(file_header)
    signals = Union{Signal{T},AnnotationsSignal}[]
    for header in read_signal_headers(io, signal_count)
        if header.label == ANNOTATIONS_SIGNAL_LABEL
            push!(signals, AnnotationsSignal(header))
        else
            push!(signals, Signal{T}(header, T[]))
        end
    end
    file_size = _size(io)
    if file_size > 0
        bytes_left = file_size - position(io)
        total_expected_samples = sum(signals) do signal
            if signal isa Signal
                return signal.header.samples_per_record
            else
                return signal.samples_per_record
            end
        end
        readable_records = div(div(bytes_left, 2), total_expected_samples)
        if file_header.record_count > readable_records
            @warn("Number of data records in file header does not match file size. " *
                  "Skipping $(file_header.record_count - readable_records) truncated " *
                  "data record(s).")
            file_header = FileHeader(file_header.version,
                                     file_header.patient,
                                     file_header.recording,
                                     file_header.start,
                                     file_header.is_contiguous,
                                     readable_records,
                                     file_header.seconds_per_record)
        end
    end
    return File(io, file_header, signals)
end

_size(io::IOStream) = filesize(io)
_size(io::IOBuffer) = io.size
_size(::IO) = -1  # type-stable unknown

"""
    EDF.is_bdf(file)

Return `true` if `file` is a BDF (BioSemi Data Format) file, otherwise `false`.
"""
is_bdf(file::File) = is_bdf(file.header)
is_bdf(header::FileHeader) = header.version == "\xffBIOSEMI"

"""
    EDF.sample_type(file::EDF.File{T})

Return the encoded type `T` of the samples stored in `file`.
"""
sample_type(file::File{T}) where {T} = T
sample_type(header::FileHeader) = is_bdf(header) ? BDF_SAMPLE_TYPE : EDF_SAMPLE_TYPE

"""
    EDF.read!(file::File)

Read all EDF sample and annotation data from `file.io` into `file.signals` and
`file.annotations`, returning `file`. If `eof(file.io)`, return `file` unmodified.
"""
function read!(file::File)
    isopen(file.io) && !eof(file.io) && read_signals!(file)
    return file
end

"""
    EDF.read(io::IO)

Return `EDF.read!(EDF.File(io))`.

See also: [`EDF.File`](@ref), [`EDF.read!`](@ref)
"""
read(io::IO) = read!(File(io))

"""
    EDF.read(path::AbstractString)

Return `open(EDF.read, path)`.
"""
read(path::AbstractString) = open(read, path)


text:
using Test
using SWIG

function simpletest()
  a = 1 ~ unif
  b = 2 ~ a
  sample(b)
end

function testlogpdf()
  x = 1 ~ ω -> Uniform(ω, 0, 1)
  y = 2 ~ ω -> Uniform(ω, 0, 1)
  z = 3 ~ (ω -> Normal(ω, x(ω), 1.0)) <| (x,)
  x_ = 0.2
  y_ = 0.4
  z_ = 0.5
  ω = Ω(Dict((1,) => x_, (2,) => y_, (3,) => z_))
  l = logpdf(rt(x, y, z), ω)
  @test l == Distributions.logpdf(Normal(x_, 1), z_)
end

function testconditioning()
  # model
  x = 1 ~ unif
  y = 2 ~ unif
  z = 3 ~ ω -> normal(ω, x(ω), 1)
  summed = x + y + z
  rt(x, y, z) | observe(normal(ω, summed, 0.1), 3.4)

  # Trace
  x_ = 0.2
  y_ = 0.4
  z_ = 0.5
  ω = Dict(1 => x_, 2 => y_, 3 => z_)
  logpdf(rt, ω)
end

function test()
  uniform(a, b) = ω -> unif(ω) * (b - a) + a

  x = 1 ~ uniform(0, 1)
  y = 2 ~ uniform(0, 1)
  d = 7 ~ y
  z = ω -> (x(ω), x(ω), y(ω), d(ω))
  @show sample(z)
  function a(ω)
    x_ = x(ω)
    d = 3 ~ uniform(0, 4)
    e = 4 ~ uniform(0, 4)
    x_ + d(ω) + e(ω)
  end

  @show sample(rt(x, y, z, a))

  # Conditional Independence
  x = 1 ~ uniform(0, 19)
  measure(ω) = x(ω) + uniform(0, 1)(ω)
  m1 = 3 ~ measure <| (x,)
  m2 = 4 ~ measure <| (x,)
  m3 = 5 ~ measure

  @show sample(rt(m1, m2, m3))

  # Linaer regression
  α = 1 ~ uniform(0, 1)
  β = 2 ~ uniform(0, 1)
  f(x, i) = i ~ (ω -> x*α(ω) + β(ω) + uniform(0.0, 1.0)(ω)) <| (α, β)
  y1 = f(0.3, 3)
  y2 = f(0.3, 4)
  @show sample(rt(y1, y2))

  y1b, y2b = pointwise() do
    f2(x, i) = i ~ (x * α + β + uniform(0.0, 1.0)) <| (α, β)
    y1b = f2(0.3, 3)
    y2b = f2(0.3, 4)
    y1b, y2b
  end

  @show sample(rt(y1, y2, y1b, y2b))

  # A bunch of random variables with a shared parents
  pointwisetest() do
    x = 1 ~ unifQD
    y = x * 100
    function a(ω)
      y(ω) + 5
    end
    function b(ω)
      y(ω) + 10
    end
    c = 2 ~ b <| (y,)
    d = 3 ~ b
    @show sample(rt(y, a, b, c, d))
  end
end


text:
@testset "Count2Reification" begin
@testset "CountEqualTo2Reification: $(fct_type), $(array_size) items, $(T)" for fct_type in ["vector of variables", "vector affine function"], array_size in [2, 3], T in [Int, Float64]
    base_model = if T == Int
        IntReificationEqualToModel{T}()
    elseif T == Float64
        FloatReificationEqualToModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.CountEqualTo2Reification{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.Count{MOI.EqualTo{T}},
    )

    x_count, _ = MOI.add_constrained_variable(model, MOI.Integer())

    if T == Int
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
    elseif T == Float64
        x_array = MOI.add_variables(model, array_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_count, x_array...])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_count, x_array...])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.Count(array_size, MOI.EqualTo(zero(T))))

    @test MOI.is_valid(model, x_count)
    for i in 1:array_size
        @test MOI.is_valid(model, x_array[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Count{MOI.EqualTo{T}}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Count{MOI.EqualTo{T}}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == 1

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == bridge.vars
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == bridge.vars_bin
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}()) == bridge.cons
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == [bridge.con_sum]
    end

    @testset "Set of variables" begin
        @test length(bridge.vars) == array_size
        @test length(bridge.vars_bin) == array_size
        for i in 1:array_size
            @test MOI.is_valid(model, bridge.vars[i])
            @test MOI.is_valid(model, bridge.vars_bin[i])
            @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_bin[i]) == bridge.vars[i]
            @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_bin[i]) == MOI.ZeroOne()
        end
    end

    @testset "Reification" begin
        @test length(bridge.cons) == array_size
        for i in 1:array_size
            @test MOI.is_valid(model, bridge.cons[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons[i])
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.cons[i]) == CP.Reification(MOI.EqualTo(zero(T)))

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == bridge.vars[i]

            t2 = f.terms[2]
            @test t2.output_index == 2
            @test t2.scalar_term.coefficient === one(T)
            @test t2.scalar_term.variable == x_array[i]
        end
    end

    @testset "Sum" begin
        @test MOI.is_valid(model, bridge.con_sum)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_sum)
        @test length(f.terms) == array_size + 1
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_sum) == MOI.EqualTo(zero(T))

        t = f.terms[1]
        @test t.coefficient === -one(T)
        @test t.variable == x_count

        for i in 1:array_size
            t = f.terms[1 + i]
            @test t.coefficient === one(T)
            @test t.variable == bridge.vars[i]
        end
    end
end
end


text:
function multiply_all(;text_data, num_data)
    f = num_data.factor
    a = num_data.num_arr
    t1 = text_data^f
    a1 = a .* f
    f1 = f^2
    return (;text_data=t1, num_data=(factor=f1, num_arr=a1))
end

fns = (;multiply_all, )


text:
abstract type InitialEnergyDistribution{N} end

"""
    (::InitialEnergyDistribution{N})(
        positions::Vector{SVectorF{N}},
        mass_densities::Vector{Float},
        eos::EquationOfState,
    ) -> Vector{Float}

Compute the fluid particle specific energies corresponding to this initial energy
distribution.
"""
function (::InitialEnergyDistribution) end

include("energy/uniform.jl")
include("energy/piecewise_uniform.jl")
include("energy/isothermal.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl