text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04gd(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    NMAX = 10
    MMAX = 10
    LDA = MMAX
    LDTAU = min(MMAX,NMAX)
    LDWORK = 3*MMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    JPVT = Array{BlasInt,1}(undef, MMAX)
    TAU = Array{Float64,1}(undef, LDTAU)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end
    if ( M<0 || M>MMAX )
        @error "Illegal M=$M"
    end

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz = M
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    JPVT[1:_isz] .= parsex.(BlasInt, vs)
    close(f)

# interp call 1

    INFO = SLICOT.mb04gd!(M, N, A, JPVT, TAU)
    @test INFO == 0
    INFO == 0 || return

    if ( INFO!=0 )
        @warn "mb04gd returns info=$INFO"
        return
    end

# interp output 1
    println(io, "JPVT:")
    _nr = M
    show(io, "text/plain", JPVT[1:_nr])
    println(io)

if ( M>=N )
    if  N>1
# interp call 2

        #FOREIGN.dlaset!( 'Lower', N-1, N-1, ZERO, ZERO, A(M-N+2,1), LDA )
        triu!(view(A,M-N+1:M,1:N-1))

    end
else
# interp call 3

    #FOREIGN.dlaset!( 'Full', M, N-M-1, ZERO, ZERO, A, LDA )
    A[1:M,1:N-M-1] .= ZERO

# interp call 4

    #FOREIGN.dlaset!( 'Lower', M, M, ZERO, ZERO, A(1,N-M), LDA )
    triu!(view(A,1:M,N-M:N-1),1)

end # if
# interp output 2
    println(io, "A:")
    _nc = N
    _nr = M
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

end # run_mb04gd()


text:
using GemmKernels.Tiling

################################################################################

@test_if "tiling" @testset "Tiling API" begin
    @testset "Tiles" begin
        @testset "Index" begin
            @test Tile(M = 4, N = 4, K = 4).index == (M = 0, N = 0, K = 0)
        end

        @testset "Projection" begin
            @test Tile(M = 1, N = 2, K = 3).MN  == Tile(M = 1, N = 2)
            @test Tile(M = 1, N = 2, K = 3).NM  == Tile(N = 2, M = 1)
            @test Tile(M = 1, N = 2, K = 3).M   == Tile(M = 1)
            @test Tile(M = 1, N = 2, K = 3).KMN == Tile(K = 3, M = 1, N = 2)
        end

        @testset "Transposition" begin
            @test transpose(Tile(M = 1, N = 2)) == Tile(N = 2, M = 1)
            @test transpose(Tile(M = 1, N = 2, K = 3)) == Tile(K = 3, N = 2, M = 1)
        end

        @testset "Translate base" begin
            tile = translate_base(Tile(M = 10, N = 20), (M = 1, N = 2))
            @test tile.size == (M = 10, N = 20)
            @test tile.base == (M = 1, N = 2)
            @test tile.offset == (M = 0, N = 0)
        end

        @testset "Translate offset" begin
            tile = translate_offset(Tile(M = 10, N = 20), (M = 1, N = 2))
            @test tile.size == (M = 10, N = 20)
            @test tile.base == (M = 0, N = 0)
            @test tile.offset == (M = 1, N = 2)
        end

        @testset "Linearise" begin
            tile = Tile(M = 3, N = 5)
            for i = 0 : 2, j = 0 : 4
                tile_t = translate_offset(tile, (M = i, N = j))
                @test linearise(tile_t.index, (M = 100, N = 200)) == j * 100 + i + 1
                @test linearise(tile_t.NM.index, (N = 200, M = 100)) == i * 200 + j + 1
            end
        end
    end

    @testset "Tile iteration" begin
        @testset "Subdivide" begin
            tile_size = (M = 8, N = 4)
            num_tiles = (M = 2, N = 4)
            tile = Tile(M = num_tiles.M * tile_size.M, N = num_tiles.N * tile_size.N)

            for i = 1 : num_tiles.M * num_tiles.N
                t = subdivide(tile, Tile(tile_size), i, num_tiles.M * num_tiles.N)

                @test t.offset == (M = 0, N = 0)
                @test t.base   == (M = tile_size.M * mod(i - 1, num_tiles.M), N = tile_size.N * fld(i - 1, num_tiles.M))
                @test t.size   == tile_size
            end
        end

        @testset "Parallellise" begin
            tile_size = (M = 8, N = 4)
            num_tiles = (M = 2, N = 8)
            tile = Tile(M = num_tiles.M * tile_size.M, N = num_tiles.N * tile_size.N)

            for i = 1 : (num_tiles.M * num_tiles.N) ÷ 2
                t1, t2 = parallellise(tile, Tile(tile_size), i, (num_tiles.M * num_tiles.N) ÷ 2)

                @test t1.offset == (M = 0, N = 0)
                @test t2.offset == (M = 0, N = 4 * tile_size.N)

                @test t1.base   == (M = tile_size.M * mod(i - 1, num_tiles.M), N = tile_size.N * fld(i - 1, num_tiles.M))
                @test t2.base   == (M = tile_size.M * mod(i - 1, num_tiles.M), N = tile_size.N * fld(i - 1, num_tiles.M))

                @test t1.size == tile_size
                @test t2.size == tile_size
            end
        end
    end
end

################################################################################

ext: jl


text:
using Documenter, AbstractLogic
# using AbstractLogic

# include("C:\\Users\\francis.smart.ctr\\GitDir\\AbstractLogic.jl\\src\\AbstractLogic.jl")

# Installation
# ] dev https://github.com/EconometricsBySimulation/AbstractLogic.jl

# Form DataFrames - using as referense

# Build documentation.
# ====================

makedocs(
    # options
    modules = [AbstractLogic],
    doctest = false,
    clean = true,
    sitename = "AbstractLogic.jl",
    format = Documenter.HTML(),
    pages = Any[
        "Introduction" => "index.md",
        "User Guide" => Any[
            "Getting Started" => "man/getting_started.md",
            "Command Snytax" => "man/command_syntax.md",
            "RELP" => "man/repl.md",
            "Generators" => "man/generators.md",
            "Operators" => "man/operators.md",
            "Super Operators" => "man/superoperators.md",
            "Meta Operators" => "man/metaoperators.md",
            "Wildcards" => "man/wildcards.md",
            "Queries" => "man/queries.md"
        ]
        ,
         "API" => Any[
        #     "Types" => "lib/types.md",
             "Functions" => "lib/functions.md"
        #     "Indexing" => "lib/indexing.md"
         ]
         ,
         "Examples" => Any[
             "Harry Potter" => "examples/harrypotter.md",
             "LSAT 2007" => "examples/LSATAnalyticalReasoning.md"
         ]
    ]

)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    repo = "github.com/EconometricsBySimulation/AbstractLogic.jl.git"
)


text:
import DataFrames
import Statistics

"""
"""
function risk_score_cutoff_values(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        multiply_by::Real = 1.0,
        average_function = Statistics.mean,
        )
    #
    ytrue = Int.(
        singlelabelbinaryytrue(
            labels_df[single_label_name],
            positive_class,
            )
        )
    #
    predictedprobabilitiesalllabels =
        predict_proba(estimator, features_df)
    yscore = Float64.(
        singlelabelbinaryyscore(
            predictedprobabilitiesalllabels[single_label_name],
            positive_class,
            )
        )
    #
    cutoffs, risk_group_prevalences = risk_score_cutoff_values(
            ytrue,
            yscore;
            multiply_by = multiply_by,
            average_function = average_function,
            )
    return cutoffs, risk_group_prevalences
end

"""
"""
function risk_score_cutoff_values(
        ytrue::AbstractVector{<:Integer},
        yscore::AbstractVector{<:AbstractFloat};
        multiply_by::Real = 1.0,
        average_function = Statistics.mean,
        )
    true_negative_rows = findall(
        ytrue .== 0
        )
    true_positive_rows = findall(
        ytrue .== 1
        )
    #
    average_score_true_negatives = average_function(
        yscore[true_negative_rows]
        )
    average_score_true_positives = average_function(
        yscore[true_positive_rows]
        )
    #
    lower_cutoff = multiply_by * average_score_true_negatives
    higher_cutoff = multiply_by * average_score_true_positives
    #
    cutoffs = (lower_cutoff, higher_cutoff,)
    #
    low_risk_group_rows = findall(
        yscore .<= average_score_true_negatives
        )
    medium_risk_group_rows = findall(
        average_score_true_negatives .<=
            yscore .<=
            average_score_true_positives
        )
    high_risk_group_rows = findall(
        average_score_true_positives .<= yscore
        )
    #
    risk_group_prevalences = DataFrames.DataFrame()
    risk_group_prevalences[:Risk_group] = [
        "Low risk",
        "Medium risk",
        "High risk",
        ]
    risk_group_prevalences[:User_supplied_average_function] = [
        average_function( ytrue[low_risk_group_rows] ),
        average_function( ytrue[medium_risk_group_rows] ),
        average_function( ytrue[high_risk_group_rows] ),
        ]
    risk_group_prevalences[:Arithmetic_mean] = [
        Statistics.mean( ytrue[low_risk_group_rows] ),
        Statistics.mean( ytrue[medium_risk_group_rows] ),
        Statistics.mean( ytrue[high_risk_group_rows] ),
        ]
    risk_group_prevalences[:Median] = [
        Statistics.median( ytrue[low_risk_group_rows] ),
        Statistics.median( ytrue[medium_risk_group_rows] ),
        Statistics.median( ytrue[high_risk_group_rows] ),
        ]
    if average_function==Statistics.mean || average_function==Statistics.median
        DataFrames.deletecols!(
            risk_group_prevalences,
            [:User_supplied_average_function],
            )
    end
    return cutoffs, risk_group_prevalences
end


text:
using Statistics

@testset "StructureFactor" begin
    p = Parameters(0.4, 1.4737, 7^3, 1e-5, 393216, Float64)
    s = SimulationSystem(p, 3, Float64)
    potential = PseudoHS(50.0, p.kT)
    initialize!(s)
    move!(10, s, potential; interval = 10)
    gfunc = PairDistributionFunction(s, 2^9)
    move!(3, s, potential, gfunc; interval = 1)
    sqfactor = StructureFactor(s, gfunc)
    compute_structure!(sqfactor, gfunc, s)

    # Should be "almost" 1
    sq_values = @view sqfactor.sq[:, 2]
    half_idx = fld1(length(sq_values), 2)
    @test isapprox(mean(sq_values[half_idx:end]), 1.0; atol=1e-1)
    # Check spacing
    dq = 0.1 * π / s.rc
    @test sqfactor.sq[2, 1] == dq
    # Check the normalizing constant
    const_val = 4.0 * π * s.ρ
    @test sqfactor.norm_const == const_val
    # And just check consistency
    @test sqfactor.nm == gfunc.nm
end


text:
orig_input = readlines("day7.txt")
Init_state = hcat(split.(orig_input,",")...) .|> x -> parse(Int64,x)

un = unique(Init_state)
positions = sort(collect(minimum(un):maximum(un)))
all_dist = zeros(Int64, size(positions,1), 2)

all_dist[:,1] = positions

# for (ind, pos) in enumerate(positions)
#     all_dist[ind,2] = sum(binomial.((abs.(Init_state .- all_dist[ind,1])).+1,2))
# end

# minimum(all_dist[:,2])

minimum([sum(binomial.((abs.(Init_state .- all_dist[ind,1])).+1,2)) for ind in (positions .+1)])


text:
# Apply boundary conditions for FV

function FV_flux_bc_rhs_only(val, facex, Qvec, t=0, dofind = 1, totaldofs = 1)
    if typeof(val) <: Number
        return val;

    elseif typeof(val) == Coefficient && typeof(val.value[1]) == GenFunction
        bvals = zeros(size(facex,2));
        if config.dimension == 1
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],0,0,t);
            end
        elseif config.dimension == 2
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],facex[2,i],0,t);
            end
        else
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],facex[2,i],facex[3,i],t);
            end
        end

    elseif typeof(val) == GenFunction
        bvals = zeros(size(facex,2));
        if config.dimension == 1
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],0,0,t);
            end
        elseif config.dimension == 2
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],facex[2,i],0,t);
            end
        else
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],facex[2,i],facex[3,i],t);
            end
        end
    end

    # Do quadrature over the face
    b = Qvec * bvals;

    return b[1];
end

function FV_flux_bc_rhs_only_simple(val, fid, t=0)
    # if typeof(val) <: Number
    #     return val;

    # elseif typeof(val) == Coefficient && typeof(val.value[1]) == GenFunction
    #     if config.dimension == 1
    #         bval = val.value[1].func(facex[1],0,0,t);
    #     elseif config.dimension == 2
    #         bval = val.value[1].func(facex[1],facex[2],0,t);
    #     else
    #         bval = val.value[1].func(facex[1],facex[2],facex[3],t);
    #     end

    # elseif typeof(val) == GenFunction
    #     if config.dimension == 1
    #         bval = val.func(facex[1],0,0,t);
    #     elseif config.dimension == 2
    #         bval = val.func(facex[1],facex[2],0,t);
    #     else
    #         bval = val.func(facex[1],facex[2],facex[3],t);
    #     end
    # end
    eid = fv_grid.face2element[1,fid];

    return FV_evaluate_bc(val, eid, fid, t);
end

function FV_copy_bdry_vals_to_vector(var, vec, grid, dofs_per_node)
    if typeof(var) <: Array
        dofind = 0;
        for vi=1:length(var)
            for compo=1:length(var[vi].symvar)
                dofind = dofind + 1;
                for bid=1:size(prob.bc_type,2)
                    if prob.bc_type[var[vi].index, bid] == DIRICHLET
                        for i = 1:length(grid.bdryface[bid]) # loop over faces with this BID
                            fid = grid.bdryface[bid][i];
                            eid = grid.face2element[1,fid];
                            vec[(eid-1)*dofs_per_node + dofind] = var[vi].values[compo, eid];
                        end
                    end
                end
            end
        end
    else
        for d=1:dofs_per_node
            dofind = d;
            for bid=1:size(prob.bc_type,2)
                if prob.bc_type[var.index, bid] == DIRICHLET
                    for i = 1:length(grid.bdryface[bid]) # loop over faces with this BID
                        fid = grid.bdryface[bid][i];
                        eid = grid.face2element[1,fid];
                        vec[(eid-1)*dofs_per_node + dofind] = var.values[d, eid];
                    end
                end
            end
        end
    end
end

# This evaluates the BC at a specific node.
# That could mean:
# - the value of constant BCs
# - evaluate a genfunction.func for BCs defined by strings->genfunctions
# - evaluate a function for BCs defined by callback functions
function FV_evaluate_bc(val, eid, fid, t)
    dim = config.dimension;
    if typeof(val) <: Number
        result = val;

    elseif typeof(val) == Coefficient
        facex = fv_info.faceCenters[:,fid];
        result = evaluate_coefficient(val, 1, facex, t, eid, fid);

    elseif typeof(val) == GenFunction
        facex = fv_info.faceCenters[:,fid];
        if dim == 1
            result=val.func(facex[1],0,0,t,eid,fid);
        elseif dim == 2
            result=val.func(facex[1],facex[2],0,t,eid,fid);
        else
            result=val.func(facex[1],facex[2],facex[3],t,eid,fid);
        end

    elseif typeof(val) == CallbackFunction
        facex = fv_info.faceCenters[:,fid];
        #form a dict for the arguments. x,y,z,t are always included
        arg_list = [];
        append!(arg_list, [("x", facex[1]),
                    ("y", dim>1 ? facex[2] : 0),
                    ("z", dim>2 ? facex[3] : 0),
                    ("t", t)]);
        # Add in the requires list
        for r in val.args
            foundit = false;
            # is it an entity?
            if typeof(r) == Variable
                foundit = true;
                if size(r.values,1) == 1
                    push!(arg_list, (string(r.symbol), r.values[1,eid]));
                else
                    push!(arg_list, (string(r.symbol), r.values[:,eid]));
                end
            elseif typeof(r) == Coefficient
                foundit = true;
                cvals = zeros(size(r.value));
                facex = fv_info.faceCenters[:,fid];
                for i=1:length(cvals)
                    cvals[i] = evaluate_coefficient(r, i, facex, t, eid, fid);
                end
                push!(arg_list, (string(r.symbol), cvals));
            elseif typeof(r) == Indexer
                foundit = true;
                push!(arg_list, (string(r.symbol), r.value));
            elseif typeof(r) == String
                # This could also be a variable, coefficient, or special thing like normal,
                if r in ["x","y","z","t"]
                    foundit = true;
                    # These are already included
                elseif r == "normal"
                    foundit = true;
                    push!(arg_list, ("normal", fv_grid.facenormals[:,fid]));
                else
                    for v in variables
                        if string(v.symbol) == r
                            foundit = true;
                            if size(v.values,1) == 1
                                push!(arg_list, (r, v.values[1,eid]));
                            else
                                push!(arg_list, (r, v.values[:,eid]));
                            end
                            break;
                        end
                    end
                    for c in coefficients
                        if string(c.symbol) == r
                            foundit = true;
                            push!(arg_list, (r, FV_evaluate_bc(c, eid, fid, t)));
                            break;
                        end
                    end
                end
            else
                # What else could it be?
            end

            if !foundit
                # Didn't figure this thing out.
                printerr("Unknown requirement for callback function "*string(fun)*" : "*string(requires[i]));
            end
        end

        # Build the dict
        args = Dict(arg_list);

        # call the function
        result = val.func(args);

    end
    return result;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module TestStaticTapirWithSequential

using FoldsTapir
using Folds.Testing: parse_tests, test_with_sequential

rawdata = """
count(isodd(x) for x in 1:10)
extrema((x - 5)^2 for x in 1:10)
mapreduce(x -> x^2, +, 1:10)
mapreduce(*, +, 1:10, 11:20)
mapreduce(*, +, 1:10, 11:20, 21:30)
maximum(0:9)
maximum(9:-1:0)
maximum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
minimum(0:9)
minimum(9:-1:0)
minimum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
prod(1:2:10)
prod(y for x in 1:11 if isodd(x) for y in 1:x:x^2; init = 1)
sum(1:10)
sum(x^2 for x in 1:11)
sum(x^2 for x in 1:11 if isodd(x); init = 0)
"""

tests = parse_tests(rawdata, @__MODULE__)

executors = [StaticTapirEx(), StaticTapirEx(nthreads = 2), StaticTapirEx(nthreads = 16)]

test_with_sequential(tests, executors)

end


text:
export Softmax

"""
Description:

    Softmax mapping between a real variable in1 ∈ R^d and discrete variable out ∈ {0, 1}^d.

    f(out, in1, xi, a) = Cat(out_j | exp(in1_j)/Σ_k exp(in1_k)),
    where log(Σ_k exp(x_k)) is upper-bounded according to (Bouchard, 2007).

Interfaces:

    1. out (discrete)
    2. in1 (real)
    3. xi  (auxiliary variable)
    4. a   (auxiliary variable)

Construction:

    Softmax(out, in1, xi, a)
"""
mutable struct Softmax <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol, Interface}

    function Softmax(out, in1, xi, a; id=generateId(Softmax))
        @ensureVariables(out, in1, xi, a)
        self = new(id, Vector{Interface}(undef, 4), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:xi] = self.interfaces[3] = associate!(Interface(self), xi)
        self.i[:a] = self.interfaces[4] = associate!(Interface(self), a)

        return self
    end
end

slug(::Type{Softmax}) = "σ"

function unsafeBoundMean(marg_in1::ProbabilityDistribution{Multivariate}, marg_xi::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Univariate})
    xi_hat = unsafeMode(marg_xi)
    a_hat = unsafeMode(marg_a)
    x_bar = unsafeMean(marg_in1)
    v_x = unsafeVar(marg_in1)
    d = length(x_bar)

    s = 0.5*(x_bar .- a_hat - xi_hat) +
        logisticLambda.(xi_hat).*(v_x + (x_bar .- a_hat).^2 - xi_hat) +
        log.(1 .+ exp.(xi_hat))

    return a_hat + sum(s)
end

# Average energy functionals
function averageEnergy(::Type{Softmax}, marg_out::ProbabilityDistribution, marg_in1::ProbabilityDistribution{Multivariate}, marg_xi::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Univariate})
    return unsafeBoundMean(marg_in1, marg_xi, marg_a) - unsafeMean(marg_out)'*unsafeMean(marg_in1)
end


text:
function init_continual_world()
	py"""
		from continualworld.envs import get_single_env
		"""
end

function get_clworld_task(name)
	pyenv = py"get_single_env"(name)
	pystate = pycall(pyenv."reset", PyObject)
	state = convert(PyArray, pystate)
	T = typeof(state)
    GymPOMDP(GymEnv{T}(Symbol(name), :v2, pyenv, pystate, state))
end

cw10() = ["hammer-v2",
        "push-wall-v2",
        "faucet-close-v2",
        "push-back-v2",
        "stick-pull-v2",
        "handle-press-side-v2",
        "push-v2",
        "shelf-place-v2",
        "window-close-v2",
        "peg-unplug-side-v2"]


text:
import BinDeps

rm′(f) = (isdir(f) || isfile(f)) && rm(f, recursive = true)

const version = "4.0.4"

folder() = normpath(joinpath(@__FILE__, "../../../deps"))

@static if Sys.isapple()
    uninstall() = map(rm′, filter(x -> !endswith(x, "build.jl"), readdir(folder())))
else
    uninstall() = rm′(joinpath(folder(), "atom"))
end

isinstalled() = Sys.isapple() ?
    isfile(joinpath(folder(), "version")) :
    isdir(joinpath(folder(), "atom"))

function install()
  dir = folder()
  if Sys.isapple()
    _icons = normpath(joinpath(@__FILE__, "../../../res/julia-icns.icns"))
  end
  !isdir(dir) && mkpath(dir)
  uninstall()
  cd(dir) do
    download(x) = run(BinDeps.download_cmd(x, basename(x)))

    download("http://junolab.s3.amazonaws.com/blink/julia.png")

    if Sys.isapple()
      file = "electron-v$version-darwin-x64.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file`)
      rm(file)
      if isdir("Julia.app") # Issue # 178
          rm("Julia.app", recursive=true)
      end
      run(`mv Electron.app Julia.app`)
      run(`mv Julia.app/Contents/MacOS/Electron Julia.app/Contents/MacOS/Julia`)
      run(`sed -i.bak 's/Electron/Julia/' Julia.app/Contents/Info.plist`)
      run(`cp $_icons Julia.app/Contents/Resources/electron.icns`)
      run(`touch Julia.app`)  # Apparently this is necessary to tell the OS to double-check for the new icons.
    end

    if Sys.iswindows()
       arch = Int == Int64 ? "x64" : "ia32"
       file = "electron-v$version-win32-$arch.zip"
       download("https://github.com/electron/electron/releases/download/v$version/$file")
       zipper = joinpath(Base.Sys.BINDIR, "7z.exe")
       if !isfile(zipper)
           #=
           This is likely built with cygwin, which includes its own version of 7z.
           But if we unzip with cmd = Cmd(`$(fwhich("bash.exe")) 7z x $file -oatom`)
           The resulting files would not be windows executable.
           So we want the 7z.exe included with a binary download of Julia.
           The PATH environment variable likely includes to the locally built
           Julia, so instead we look in the default Julia binary location and
           pick the latest version.
           =#
           juliafolders = filter(readdir(ENV["LOCALAPPDATA"])) do f
                                     startswith(f, "Julia-")
              end
           juliaversions = VersionNumber.([replace(f, "Julia-" => "") for f in juliafolders])
           i = findlast(isequal(maximum(juliaversions)), juliaversions)
           zipper = joinpath(ENV["LOCALAPPDATA"], juliafolders[i], "bin", "7z.exe")
           isfile(zipper) || error("could not find $zipper. Try also installing a binary version of Julia.")
        end
        cmd = Cmd([zipper, "x", file, "-oatom"])
        run(cmd)
        rm(file)
    end

    if Sys.islinux()
      arch = Int == Int64 ? "x64" : "ia32"
      file = "electron-v$version-linux-$arch.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file -d atom`)
      rm(file)
    end
  end
end


text:
using Weave
weave("notebooks/*.jmd", out_path="notebooks", doctype="github")


text:
module Archimed

import XMLDict.xml_dict
import EzXML.readxml
import EzXML.root
import Colors.RGBA
import LinearAlgebra.I
import StaticArrays.SMatrix

using MeshIO
using GeometryTypes

# 3D:
export read_opf,extract_opf_ref_meshes,extract_opf_materials
export meshes_from_topology,merge_meshes
export material,extract_opf_shapes,meshes_color

# For debbug:
export attr_type
export mesh_from_topology!

include("3D_related/read_opf.jl")
include("3D_related/meshes.jl")
include("3D_related/materials.jl")
include("3D_related/shapes.jl")
include("3D_related/mesh_colors.jl")
include("3D_related/3D_types.jl")

end # module


text:
using ChainRulesTestUtils
using Documenter

makedocs(;
    modules=[ChainRulesTestUtils],
    format=Documenter.HTML(; prettyurls=false, assets=["assets/chainrules.css"]),
    sitename="ChainRulesTestUtils",
    authors="JuliaDiff contributors",
    pages=[
        "ChainRulesTestUtils" => "index.md",
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
   )

const repo = "github.com/JuliaDiff/ChainRulesTestUtils.jl.git"
deploydocs(; repo=repo, push_preview=true)


text:
export bandwidth_ratio

function fixlength(x,n;fill=0)
  if length(x) > n
    view(x,1:n)
  elseif length(x) < n
    append!(convert(Array{Union{eltype(x),Missing}},x),
            Base.fill(fill,n-length(x)))
  else
    x
  end
end

function dominant_channels(sp;level_threshold=0.9,threshold=0.25,
                           window=500ms,delta=250ms,thresh_ratio=1.0,
                            max_channels=10)
  x = map_windowing(sp,length=window,step=delta,flatten=true) do window
    thresh_level = thresh_ratio*quantile(vec(sp),1 - threshold/size(sp,2))
    # @show thresh_level
    levels = dropdims(mapslices(x -> quantile(x,level_threshold),
                                Array(window),dims=1), dims=1)
    # @show sort(levels,rev=true)[1:min(end,20)]
    # thresh_level = thresh_ratio*
    #   quantile(levels,1 - threshold/length(levels))

    x = fixlength(findall(levels .> thresh_level),max_channels,fill=missing)
    # @show x
    x
  end
end

function bandwidth_ratio(spmask::AbstractMatrix, sp::AbstractMatrix,
                         settings)
  settings = read_settings(settings)
  startHz, stopHz = asHz.(settings.rates.freqbound.freq_limits)
  bandwidth_ratio(spmask, sp[:,startHz .. stopHz]; settings.bandwidth_ratio...)
end

function bandwidth_ratio(spmask, sp; threshold=1.5,
                         window=500ms,
                         full_band_ratio=2,
                         level_threshold=0.9,
                         thresh_ratio=1.0,
                         delta=250ms)
  window = asseconds(window)
  @assert frequencies(spmask) == frequencies(sp) "Frequency axes are not equal."
  @assert times(spmask) == times(sp) "Time axes are not equal."

  fullchan = dominant_channels(sp,window=window*full_band_ratio,
                               delta=delta,threshold=threshold,
                               level_threshold=level_threshold,
                               thresh_ratio=thresh_ratio)

  maskchan = dominant_channels(spmask,window=window,delta=delta,
                               threshold=threshold,
                               level_threshold=level_threshold,
                               thresh_ratio=thresh_ratio)

  dim = axisdim(fullchan,Axis{:time})
  others = setdiff(1:ndims(fullchan),[dim])
  fullband = mapslices(Array(fullchan),dims=others) do x
    maximum(skipmissing(x)) - minimum(skipmissing(x))
  end |> x -> AxisArray(vec(x),Axis{:time}(times(fullchan)))

  maskband = AxisArray(Array{Float64}(undef,size(maskchan,dim)),
                       Axis{:time}(times(maskchan)))
  ratio = AxisArray(Array{Float64}(undef,size(maskchan,dim)),
                    Axis{:time}(times(maskchan)))
  winwidth = window*full_band_ratio
  for (i,t) in enumerate(times(maskchan))
    val = (-window*full_band_ratio*0.51 .. window*full_band_ratio*0.51) + t
    maskwin = view(maskchan,Axis{:time}(t))
    fullwin = view(fullchan,Axis{:time}(val))
    maskcount = skipmissing(intersect(maskwin,fullwin))
    if isempty(maskcount)
      maskband[i] = 0
    else
      maskband[i] = maximum(maskcount) - minimum(maskcount)
    end
    ratio[i] = maskband[i] / maximum(view(fullband,Axis{:time}(val)))
  end

  ratio, fullband, maskband
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
"""
Parallel Factor Analysis (PARAFAC) decomposition
"""
immutable PARAFAC2{T<:Number, N} <: TensorDecomposition{T, N}
    factors::NTuple{N, Matrix{T}}
    D::Vector{Matrix{T}}
    A::Matrix{T}
    props::Dict{Symbol, Any}

    function Base.call{T}(::Type{PARAFAC2},
                   X::Vector{Matrix{T}},
                   F::Matrix{Float64},
                   D::Vector{Matrix{Float64}},
                   A::Matrix{Float64})
        factors = (map(function (Xi, Di)
            U = svd(A_mul_Bt(F .* Di, Xi * A))
            U[3] * (U[1]'F)
        end, X, D)...)
        return new{T, length(factors)}(factors, D, A, Dict{Symbol,Any}())
    end

end

"""
PARAFAC2 model
"""
function parafac2{S<:Matrix}(X::Vector{S},
                             r::Integer;
                             tol::Float64=1e-5,
                             maxiter::Integer=100,
                             compute_error::Bool=false,
                             verbose::Bool=true)
    m = length(X)
    n = size(X[1], 2)
    for i in 2:m
        size(X[i], 2) == n || error("All input matrices should have the same number of rows.")
    end

    F = eye(r)
    D = Matrix{Float64}[ones(1, r) for _ in 1:m]
    A = eigs(sum(map(Xi -> Xi'Xi, X)), nev=r)[2]
    G = Matrix{Float64}[eye(r), eye(r), ones(r, r) * m]
    H = Matrix{Float64}[(size(X[i], 1) > size(X[i], 2) ? qr(X[i])[2]: X[i]) for i in 1:m]
    P = Array(Matrix{Float64}, m)

    niters = 0
    converged = false
    resid = vecnorm(vcat(X...))
    while !converged && niters < maxiter
        map!(function (Hi, Di)
            U = svd((F .* Di) * (Hi * A)')
            U[3]*U[1]'
        end, P, H, D)
        T = cat(3, [P[i]'H[i] for i in 1:m]...)

        B = vcat(D...)
        F = _row_unfold(T, 1) * khatrirao(B, A) / (G[3] .* G[2])
        At_mul_B!(G[1], F, F)
        A = _row_unfold(T, 2) * khatrirao(B, F) / (G[3] .* G[1])
        At_mul_B!(G[2], A, A)
        B = _row_unfold(T, 3) * khatrirao(A, F) / (G[2] .* G[1])
        At_mul_B!(G[3], B, B)

        D = Matrix{Float64}[B[i, :] for i in 1:m]

        resid_old = resid
        resid = sum(k -> sumabs2(H[k] - P[k] * A_mul_Bt(F .* D[k], A)), 1:m)
        converged = abs(resid - resid_old) < tol * resid_old

        niters += 1
    end

    verbose && _iter_status(converged, niters, maxiter)
    res = PARAFAC2(X, F, D, A)
    if compute_error
        _set_rel_residue(res, sqrt(resid) / vecnorm(vcat(X...)))
    end
    return res
end


text:
using IntervalArithmetic
using Test

A = [ 2..4   -2..1
     -1..2    2..4]

b = [-2..2
     -2..2]

@testset "Linear algebra with intervals tests" begin

    @test A * b == [-12..12
                    -12..12]

    # Example from Moore et al., Introduction to Interval Analysis (2009), pg. 88:

    @test A \ b == [-5..5
                    -4..4]

end


text:
push!(LOAD_PATH, "src/")
using ArgParse,XAM, BioAlignments, nanoJulia, DataFrames, CSV

# Arguments parser
function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table! s begin
		"--inputbam", "-i"
			help = "The path to BAM file"
			arg_type = String
			required = true
		"--outputdir", "-o"
			help = "Folder to output QC results"
			required = true
	end

	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	inputbamfile = parsed_args["inputbam"]
    outputdir = parsed_args["outputdir"]
    bamindexfile = string(inputbamfile, ".bai")
    if !isfile(bamindexfile)
        println("\033[1;31m[Error] Please make sure if the index of input BAM file exists!\033[0m")
        exit(1)
    end

	if isdir(outputdir)
		println("\033[1;33m[Warning] The output folder has already existed, so it will be renamed with additional suffix \"_bak\"\033[0m")
		new_outputdir = string(outputdir, "_bak")
		mv(outputdir, new_outputdir)
	end
	mkdir(outputdir)

	length2qualityTextFile = joinpath(outputdir, "readlength_vs_readquality.tsv")

    println("\033[1;32m* Start parsing this BAM file...\033[0m")
	output_table = readBAM(inputbamfile)
	println("\033[1;32m* Writing output file...\033[0m")
	output_table |> CSV.write(length2qualityTextFile)
	totalLength = nanoJulia.totalLen(output_table.length)
	n50 = nanoJulia.calculate_N50(output_table.length, totalLength)
	println("\033[1;32m* Generating summry file and figures...\033[0m")
	generateStatSummary(output_table, totalLength, n50, outputdir)
end

main()


text:
using jInv.InverseSolve
using jInv.LinearSolvers
using Statistics
using jInv.Utils
using Distributed

function calculateDobs(m::Array{Float64, 2}, pForp::Array{RemoteChannel},
	omega::Vector, nfreq::Int64)
	Dp, = getData(vec(m),pForp)
	nfreq = length(omega);

	for k=1:length(Dp)
		wait(Dp[k]);
	end
	Dobs = Array{Array}(undef, length(Dp));
	for k=1:length(Dp)
		Dobs[k] = fetch(Dp[k]);
	end

	return Dobs;
end

export solveForwardProblem
function solveForwardProblem(m::Array{Float64, 2}, pForp::Array{RemoteChannel},
	omega::Vector, nfreq::Int64)
	Dobs = calculateDobs(m, pForp, omega, nfreq)
	Wd = Array{Array}(undef, nfreq);

	for k=1:length(Dobs)
		Wd[k] = 1.0./(abs.(real.(Dobs[k])) .+ 1e-1*mean(abs.(Dobs[k])));
	end

	return Dobs, Wd;
end

export solveForwardProblemExtendedSources
function solveForwardProblemExtendedSources(m::Array{Float64, 2}, pForp::Array{RemoteChannel},
	omega::Vector, nfreq::Int64)
	Dobs = calculateDobs(m, pForp, omega, nfreq)
	Wd = Array{Array}(undef, length(Dobs));

	avgDobs = mean(abs.(map(dobs_j -> mean(abs.(dobs_j)), Dobs[:])));
	for k=1:length(Dobs)
		Wd[k] = ones(size(Dobs[k])) ./ avgDobs ;
	end

	return Dobs, Wd;
end

export solveForwardProblemNoProcs
function solveForwardProblemNoProcs(m::Array{Float64, 2}, pFor::BasicFWIparam,
	omega::Vector, nfreq::Int64)
	Dp,pFor = getData(vec(m),pFor)
	nfreq = length(omega);
	Dobs = fetch(Dp);
	Wd = 1.0./(Dobs .+ 1e-1*mean(abs.(Dobs)));

	return Dobs, Wd;
end

function wFourthOrderSmoothing(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
	dm = m.-mref;
	d2R = wdiffusionReg(m,mref,M,Iact = Iact,C = C)[3];
	clear!(M);
	d2R = d2R'*d2R;
	dR  = d2R*dm;
	Rc  = 0.5*dot(dm,dR);
   return Rc,dR,d2R
end

function getFreqContParams(pFor::Array{RemoteChannel}, Dobs::Array, Wd::Array,
	nfreq::Int64, nx::Int64, nz::Int64, mref::Array{Float64,2}, Mr::RegularMesh,
	boundsHigh::Float64, boundsLow::Float64, plotting::Bool=false,
	plottingFunc::Function=dummy)
	N = prod(Mr.n);
	Iact = sparse(I,N,N);
	mback   = zeros(Float64,N);

	########################################################################################################
	##### Set up remote workers ############################################################################
	########################################################################################################

	Ainv = getJuliaSolver();
	## Choose the workers for FWI (here, its a single worker)
	misfun = SSDFun;
	probsMax = ceil(Integer,nfreq/nworkers());
	nprobs   = zeros(maximum(workers()));
	pMis = getMisfitParam(pFor, Wd, Dobs, misfun, Iact, mback);
	boundsHigh = boundsHigh*ones(Float32,N);
	boundsLow = boundsLow*ones(Float32,N);

	maxStep				=0.05*maximum(boundsHigh);
	# regparams 			= [1.0,1.0,1.0,1e-5];
	cgit 				= 8;
	alpha 				= 1e+1;
	pcgTol 				= 1e-3;
	maxit 				= 5;
	HesPrec 			= getExactSolveRegularizationPreconditioner();
	# regfun(m,mref,M) 	= wdiffusionReg(m,mref,M,Iact=Iact,C=[]);
	regfun(m,mref,M) 	= wFourthOrderSmoothing(m,mref,M,Iact=Iact,C=[]);

	pInv = getInverseParam(Mr,identityMod,regfun,alpha,mref[:],boundsLow,boundsHigh,
	                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
							 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);

	function plotIntermediateResults(mc,Dc,iter,pInv,PMis, resultsFilename::String)
		# Models are usually shown in velocity.
		fullMc = reshape(Iact*pInv.modelfun(mc)[1] + mback,tuple((pInv.MInv.n)...));
		Temp = splitext(resultsFilename);
		if iter>0
			Temp = string(Temp[1],iter,Temp[2]);
		else
			Temp = resultsFilename;
		end
		if resultsFilename!=""
			DelimitedFiles.writedlm(Temp,convert(Array{Float16},fullMc));
		end
		if plotting
			plottingFunc(fullMc');
		end
	end
	return pInv, pMis, plotIntermediateResults, Iact, mback;
end

export solveInverseProblemExtendedSources
function solveInverseProblemExtendedSources(pFor::Array{RemoteChannel}, Dobs::Array, Wd::Array,
	nfreq::Int64, nx::Int64, nz::Int64, mref::Array{Float64,2}, Mr::RegularMesh,
	boundsHigh::Float64, boundsLow::Float64, resultsFilename::String, plotting::Bool=false,
	plottingFunc::Function=dummy)

	regfunFourthOrder(m,mref,M) = wFourthOrderSmoothing(m,mref,M,Iact=Iact,C=[]);
	pInv, pMis, plotIntermediateResults, Iact, mback = getFreqContParams(pFor, Dobs, Wd, nfreq, nx, nz, mref, Mr,
		boundsHigh, boundsLow, plotting, plottingFunc);

	regfunDiff(m,mref,M) = wdiffusionReg(m,mref,M,Iact=Iact,C=[]);

	# Run one sweep of a frequency continuation procedure.
	mc, Dc = freqContExtendedSources(copy(mref[:]), pInv, pMis[1:2], 2, 4,Iact,mback, plotIntermediateResults, resultsFilename, 1);

	pInv.regularizer = regfunDiff;
	# Run one sweep of a frequency continuation procedure.
	mc, Dc = freqContExtendedSources(copy(mc[:]), pInv, pMis, nfreq, 4,Iact,mback, plotIntermediateResults, resultsFilename, 1);

	return mc, Dc, pInv, Iact, mback, map(x->fetch(x), pMis);
end

export solveInverseProblem
function solveInverseProblem(pFor::Array{RemoteChannel}, Dobs::Array, Wd::Array,
	nfreq::Int64, nx::Int64, nz::Int64, mref::Array{Float64,2}, Mr::RegularMesh,
	boundsHigh::Float64, boundsLow::Float64, resultsFilename::String, plotting::Bool=false,
	plottingFunc::Function=dummy)

	regfunFourthOrder(m,mref,M) = wFourthOrderSmoothing(m,mref,M,Iact=Iact,C=[]);
	pInv, pMis, plotIntermediateResults, Iact, mback = getFreqContParams(pFor, Dobs, Wd, nfreq, nx, nz, mref, Mr,
		boundsHigh, boundsLow, plotting, plottingFunc);

	regfunDiff(m,mref,M) = wdiffusionReg(m,mref,M,Iact=Iact,C=[]);

	# Run one sweep of a frequency continuation procedure.
	mc, Dc = freqContBasic(copy(mref[:]), pInv, pMis[1:2], 2, 4, plotIntermediateResults, resultsFilename, 1);

	pInv.regularizer = regfunDiff;
	# Run one sweep of a frequency continuation procedure.
	mc, Dc = freqContBasic(copy(mc[:]), pInv, pMis, nfreq, 4, plotIntermediateResults, resultsFilename, 1);

	return mc, Dc, pInv, Iact, mback;
end

export solveInverseProblemTraceEstimation
function solveInverseProblemTraceEstimation(pFor::Array{RemoteChannel}, Dobs::Array, Wd::Array,
	nfreq::Int64, nx::Int64, nz::Int64, mref::Array{Float64,2}, Mr::RegularMesh,
	boundsHigh::Float64, boundsLow::Float64, resultsFilename::String, plotting::Bool=false,
	plottingFunc::Function=dummy)
	pInv, pMis, plotIntermediateResults, Iact, mback = getFreqContParams(pFor,
		Dobs, Wd, nfreq, nx, nz, mref, Mr,
		boundsHigh, boundsLow, plotting, plottingFunc);

	# Run one sweep of a frequency continuation procedure.
	mc, Dc = freqContTraceEstimation(copy(mref[:]), pInv, pMis, nfreq, 4,
			Iact, mback, plotIntermediateResults, resultsFilename, 1);

	return mc, Dc, pInv, Iact, mback;
end

export solveInverseProblemNoProcs
function solveInverseProblemNoProcs(pFor::BasicFWIparam, Dobs::Array, Wd::Array, nfreq::Int64,
	nx::Int64, nz::Int64, mref::Array{Float64,2}, Mr::RegularMesh,
	boundsHigh::Float64, boundsLow::Float64, plotting::Bool=false, plottingFunc::Function=dummy)
	N = prod(Mr.n);
	Iact = sparse(I,N,N);
	mback   = zeros(Float64,N);

	########################################################################################################
	##### Set up remote workers ############################################################################
	########################################################################################################

	Ainv = getJuliaSolver();
	## Choose the workers for FWI (here, its a single worker)
	misfun = SSDFun;
	probsMax = ceil(Integer,nfreq/nworkers());
	nprobs   = zeros(maximum(workers()));

	pMis = getMisfitParam(pFor, Wd, Dobs, misfun);

	boundsHigh = boundsHigh*ones(Float32,N);
	boundsLow = boundsLow*ones(Float32,N);

	maxStep				=0.05*maximum(boundsHigh);

	#regparams 			= [1.0,1.0,1.0,1e-5];
	cgit 				= 8;
	alpha 				= 1e+2;
	pcgTol 				= 1e-1;
	maxit 				= 30;
	HesPrec 			= getExactSolveRegularizationPreconditioner();
	regfun(m,mref,M) 	= wdiffusionReg(m,mref,M,Iact=Iact,C=[]);

	pInv = getInverseParam(Mr,identityMod,regfun,alpha,mref[:],boundsLow,boundsHigh,
	                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
							 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);

	 function plotIntermediateResults(mc,Dc,iter,pInv,PMis)
	 	# Models are usually shown in velocity.
	 	fullMc = reshape(Iact*pInv.modelfun(mc)[1] + mback,tuple((pInv.MInv.n)...));
	 	if plotting
			plottingFunc(fullMc');
	 	end
	 end

	mc, Dc, = projGNCG(copy(mref[:]),pInv,pMis);
	return mc, Dc, pInv, Iact, mback;
end


text:
#Visualization Functions

export visualize

using Plots

gr()

macro _call_func_(f,c, kwargs...)
    z = [esc(a) for a in kwargs]
    return :((eval($f))($c; $(z...)))
end

function _plot(c::confusion_matrix; func = nothing, type = nothing, title = "Visualization", labels = nothing)
    x = @_call_func_ Symbol(type)  c average = "binary"
    l = Array{typeof(labels[1])}(labels)
    x = func(l, x, title = title, labels = permutedims(l))
end

"""
```visualize(c::confusion_matrix, <keywords>)```

Visualize the given properties of the confusion matrix as specified

## Keywords
- `mode` : String, Array of String; represents properties given below of the confusion matrix, can be either an array containing different properties
or a single string.\n
    Supported Modes: \n
    - `matrix`\n
    - `condition_positive`\n
    - `condition_negative`\n
    - `predicted_positive\n
    - `predicted_negative`\n
    - `correctly_classified`\n
    - `incorrectly_classified`\n
    - `sensitivity_score`\n
    - `recall_score`\n
    - `specificity_score`\n
    - `precision_score`\n
    - `positive_predictive-value`\n
    - `accuracy_score`\n
    - `balanced_accuracy-score`\n
    - `negative_predictive-value`\n
    - `false_negative_rate`\n
    - `false_positive_rate`\n
    - `false_discovery_rate`\n
    - `false_omission_rate`\n
    - `f1_score`\n
    - `prevalence_threshold`\n
    - `threat_score`\n
    - `matthews_correlation_coeff`\n
    - `fowlkes_mallows_index`\n
    - `informedness`\n
    - `markedness`\n
    - `cohen_kappa_score`\n
    - `hamming_loss`\n
    - `jaccard_score`\n

- `seriestype::String = "heatmap"` :
    Supported visualization functions:
        - `heatmap`
        - `bar`
        - `histogram`
        - `scatter`
        - `line`

`title::String` : Denotes the title that will displayed above the drawn plot, default: nothing

`labels::Vector` : Denotes the labels that will be used for the plot. If equals to nothing, labels of the given confusion matrix will be used.
"""
function visualize(c::confusion_matrix; mode = "matrix", seriestype::String = "heatmap", title= nothing, labels = nothing)
    @assert seriestype in ["scatter", "heatmap", "line", "histogram", "bar"] "Unknown visualization format"
    labels = labels != nothing ? labels : convert(Array{typeof(c.Labels[1])}, c.Labels)
    if title == nothing; title = mode isa Array ? mode : String(Base.copymutable(mode)); end
    title = title isa Array ? title : [title]
    mode = mode isa Array ? mode : [mode]
    plt = []
    for i in 1:length(mode)
        @assert mode[i] in ["matrix", "condition_positive", "condition_negative", "predicted_positive","predicted_negative", "correctly_classified", "incorrectly_classified", "sensitivity_score", "recall_score", "specificity_score", "precision_score", "positive_predictive_value", "accuracy_score", "balanced_accuracy_score", "negative_predictive_value", "false_negative_rate", "false_positive_rate", "false_discovery_rate",
         "false_omission_rate", "f1_score", "prevalence_threshold", "threat_score", "matthews_correlation_coeff", "fowlkes_mallows_index",
         "informedness", "markedness", "cohen_kappa_score", "hamming_loss", "jaccard_score"] "Unknown visualization mode"
        if mode[i] != "matrix"; @assert seriestype in ["scatter", "line", "histogram", "bar"] "The given mode does not support this visualization format"; end
        x = nothing
        if mode[i] == "matrix"
            if seriestype == "histogram"; x = histogram(labels, c.matrix, labels = permutedims(labels), title = title[i])
            elseif seriestype == "scatter"; x = scatter(labels, c.matrix, labels = permutedims(labels), title = title[i])
            elseif seriestype == "line"; x = plot(labels, c.matrix, labels = permutedims(labels), title = title[i])
            elseif seriestype == "bar"; x = bar(labels, c.matrix, labels = permutedims(labels), title = title[i])
            elseif seriestype == "heatmap"; x = heatmap(labels, labels, c.matrix, labels = permutedims(labels), title = title[i])
            end
        else
            if seriestype == "histogram"; x = _plot(c; func = histogram, type = Symbol(mode[i]), title = title[i], labels = labels)
            elseif seriestype == "scatter"; x = _plot(c; func = scatter, type = Symbol(mode[i]), title = title[i], labels = labels)
            elseif seriestype == "line"; x =  _plot(c; func = plot, type = Symbol(mode[i]), title = title[i], labels = labels)
            elseif seriestype == "bar"; x =  _plot(c; func = bar, type = Symbol(mode[i]), title = title[i], labels = labels)
            elseif seriestype == "heatmap"; x =  _plot(c; func = heatmap, type = Symbol(mode[i]), title = title[i], labels = labels)
            end
        end
        push!(plt, x)
    end
    plot(plt..., layout = (length(plt), 1))
end

visualize(r::RocCurve, title = "Roc Curve", marker = (15, 0.2, :orange)) = scatter(r.ratios, annotations = [(r.ratios[i]..., r.thresholds[i]) for i in 1:length(r.ratios)], title = title, marker = marker)
visualize(r::PrecisionRecallCurve, title = "Precision Recall Curve", marker = (15, 0.2, :orange)) = scatter(r.ratios, annotations = [(r.ratios[i]..., r.thresholds[i]) for i in 1:length(r.ratios)], title = title, marker = marker)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletgeometryeditor

"""
    ''_leafletgeometryeditor(;kwargs...)
    ''_leafletgeometryeditor(children::Any;kwargs...)
    ''_leafletgeometryeditor(children_maker::Function;kwargs...)

A LeafletGeometryEditor component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `className` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `style` (Dict; optional)
"""
function ''_leafletgeometryeditor(; kwargs...)
        available_props = Symbol[:children, :id, :className, :loading_state, :style]
        wild_props = Symbol[]
        return Component("''_leafletgeometryeditor", "LeafletGeometryEditor", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

''_leafletgeometryeditor(children::Any; kwargs...) = ''_leafletgeometryeditor(;kwargs..., children = children)
''_leafletgeometryeditor(children_maker::Function; kwargs...) = ''_leafletgeometryeditor(children_maker(); kwargs...)

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
using MolecularIntegrals

open("src/HGPgen3.jl","w") do io
    for amax in 0:3
        for cmax in 0:(amax-1)
            f = MolecularIntegrals.vrr_autogen(amax,cmax)
            write(io,f)
            f = MolecularIntegrals.vrr_autogen(cmax,amax)
            write(io,f)
        end
        f = MolecularIntegrals.vrr_autogen(amax,amax)
        write(io,f)
    end
end


text:
using Statistics

#g::Global, p::Population,
function epoch(state::SimulationState, n::Int, report::Bool=true, save_best::Bool=false,
               checkpoint_interval::Int=15, checkpoint_generation=0)
    #= Runs NEAT's genetic algorithm for n epochs.
        Keyword arguments:
        report -- show stats at each epoch (default true)
        save_best -- save the best chromosome from each epoch (default False)
        checkpoint_interval -- time in minutes between saving checkpoints (default 15 minutes)
        checkpoint_generation -- time in generations between saving checkpoints
            (default 0 -- option disabled)
    =#
    Random.seed!(123456)
    t0 = time()
    tc = time() # for saving checkpoints
    #fitness = Evaluation(g.cf.dataset_id, g.cf.fitness_function_id)
    #p.evaluate()

    # Evaluate individuals
    #p.evaluate(p.setup)
    state.evaluate!(state.setup)
    g = state.g
    p = state.population

    initial_generation = p.generation == 0 ? 0 : 1

    for gen in initial_generation:n
        p.generation += 1

        print("\n ****** Running generation $(p.generation) ******")

        if report println("\n ****** Running generation $(p.generation) ******") end

        # Speciates the population
        speciate(g, p, report)

        # Current population's average fitness
        push!(p.avg_fitness, average_fitness(p))

        # Current generation's best chromosome
        bestfit, bestidx = findmax(map(ch-> ch.fitness, p.population))
        best = p.population[bestidx]
        previous_best_fitness = p.generation == 0 ? 0 : p.best_fitness[end].fitness
        push!(p.best_fitness, best)

        if previous_best_fitness == best.fitness
            p.best_no_improvement_age += 1
        else
            p.best_no_improvement_age = 0
        end

        # Which species has the best chromosome?
        for s in p.species
            s.hasBest = false
            if best.species_id == s.id
                s.hasBest = true
            end
        end

        limit_species!(g, p)

        #-----------------------------------------
        # Prints chromosome's parents id:  {dad_id, mon_id} -> child_id
#         map(ch-> @printf("{%3d; %3d} -> %3d   Nodes %3d   Connections %3d\n",
#                          ch.parent1_id, ch.parent2_id, ch.id, size(ch)[1], size(ch)[2]), p.population)
        #-----------------------------------------

        remove_stagnated_species!(p, g, report)

        # Compute spawn levels for each remaining species
        compute_spawn_levels(g, p)

        remove_unspawned_species!(p, report)

        # Logging speciation stats
        log_species(p)

        print_debugging_information(report, p, best)

        # Stops the simulation
        if best.fitness >= g.cf.max_fitness_threshold
            @printf("Best individual found in epoch %s - complexity: %s\n", p.generation, number_nodes(best))
            break
        end

        # -------------------------- Producing new offspring -------------------------- #
        new_population = Tuple{MLPChromosome, Bool}[] # next generation's population

        # Spawning new population
        for s in p.species new_population = vcat(new_population, reproduce(g, s, state.setup)) end

        update_population!(p, g, state.setup, new_population, report)

        if time() > tc + 60 * checkpoint_interval
            create_checkpoint(g, p, report)
            tc = time() # updates the counter
        elseif  checkpoint_generation != 0 && p.generation % checkpoint_generation == 0
            create_checkpoint(g, p, report)
        end
    end
    t1 = time()
    Δt = t1 - t0
    p.simulation_time += Δt
    t2 = time()
    create_checkpoint(g, p, report)
    total_time = Δt + (time() - t2)
    println("Time: $total_time seconds.")
end

function remove_stagnated_species!(p::Population, g::Global, report)
    # Remove stagnated species and its members (except if it has the best chromosome)
    number_species = length(p.species)
    speciesToKeep = trues(number_species)
    deletedSpeciesIds = Int[]

    if number_species > 2 && p.best_no_improvement_age > g.cf.best_max_stagnation
        sort!(p.species, by = s-> s.last_avg_fitness, rev = true)
        speciesToKeep = map((s) -> s[2].hasBest || s[1] == 1 || s[1] == 2, enumerate(p.species))
        #speciesToKeep[3:end] .= false
        p.best_no_improvement_age = 0
        append!(deletedSpeciesIds, map((s) -> s.id, p.species[.!speciesToKeep]))
    else
        number_kept_species = number_species
        for i = 1:length(p.species)
            if p.species[i].no_improvement_age > g.cf.max_stagnation
                if (!p.species[i].hasBest && p.species[i].no_improvement_age > 2 * g.cf.max_stagnation) && number_kept_species > 2
                    if report @printf("\n   Species %2d age %2s (with %2d individuals) is stagnated: removing it",
                                      p.species[i].id, p.species[i].age, length(p.species[i])) end
                    speciesToKeep[i] = false
                    number_kept_species -= 1
                    push!(deletedSpeciesIds, p.species[i].id)
                end
            end
        end
    end

    #if length(p.species) > 1 p.species = p.species[speciesToKeep] end
    p.species = p.species[speciesToKeep] # prune unwanted species

    # remove species' chromosomes from population
    # prune unwanted chromosomes
    p.population = filter((i) -> i.species_id ∉ deletedSpeciesIds, p.population)
end

function remove_unspawned_species!(p::Population, report)
    # Removing species with spawn amount = 0
    speciesToKeep = trues(length(p.species))
    deletedSpeciesIds = Int64[]
    for i = 1:length(p.species)

        # This rarely happens
        if p.species[i].spawn_amount == 0
            if report @printf("\n   Species %2d age %2s removed: produced no offspring",p.species[i].id, p.species[i].age) end
            speciesToKeep[i] = false
            push!(deletedSpeciesIds,p.species[i].id)
        end
    end
    p.species = p.species[speciesToKeep] # prune unwanted species

    # remove species' chromosomes from population
    p.population = filter((i) -> i.species_id ∉ deletedSpeciesIds, p.population)  # prune unwanted chromosomes
end

function print_debugging_information(report::Bool, p::Population, best::MLPChromosome)
    if report
        @printf("\nPopulation's average fitness: %3.5f stdev: %3.5f", p.avg_fitness[end], stdeviation(p))
        @printf("\nBest fitness: %2.12s - size: %s - species %s - id %s", best.fitness, number_nodes(best), best.species_id, best.id)

        # print some "debugging" information
        @printf("\nSpecies length: %d totalizing %d individuals", length(p.species), sum([length(s) for s in p.species]))
        @printf("\nSpecies ID       : %s",   [s.id for s in p.species])
        @printf("\nEach species size: %s",   [length(s) for s in p.species])
        @printf("\nAmount to spawn  : %s",   [s.spawn_amount for s in p.species])
        @printf("\nSpecies age      : %s",   [s.age for s in p.species])
        @printf("\nSpecies no improv: %s\n", [s.no_improvement_age for s in p.species]) # species no improvement age

        for s in p.species println(s) end
    end
end

function update_population!(p::Population, g::Global, fitness::Function, new_population, report)
    # ----------------------------------------------#
    # Controls target population under or overflow  #
    # ----------------------------------------------#
    fill = p.popsize - length(new_population)
    if fill < 0 # overflow
        if report println("\n   Removing $(abs(fill)) excess individual(s) from the new population") end
        # TODO: This is dangerous? I can't remove a species' representant!
        #new_population = new_population[1:end+fill] # Removing the last added members
        while fill < 0
            pos = rand((1:length(new_population)))
            # checks flag variable which indicates the individuals that must be kept
            if !new_population[pos][2]
                deleteat!(new_population, pos)
                fill += 1
            end
        end
    end

    new_population = getindex.(new_population, 1) # Delete flag variables

    if fill > 0 # underflow
        if report println("\n   Producing $fill more individual(s) to fill up the new population") end

        # TODO:
        # what about producing new individuals instead of reproducing?
        # increasing diversity from time to time might help
        while fill > 0
            # Selects a random chromosome from population
            parent1 = p.population[rand(1:length(p.population))]
            # Search for a mate within the same species
            found = false
            for parent2 in p.population
                if parent2.species_id == parent1.species_id && parent2.id != parent1.id
                    child = mutate!(crossover(g, parent1, parent2),g)
                    fitness(child)
                    push!(new_population, child)
                    found = true
                    break
                end
            end
            if !found
                child = mutate!(deepcopy(parent1),g)
                fitness(child)
                push!(new_population, child) # will irreversibly mutate parent. ok?
            end # If no mate was found, just mutate it

            fill -= 1
        end
    end

    @assert p.popsize == length(new_population) # Different population sizes!

    # Updates current population
    p.population = new_population
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
#semicircle.m
#Algorithm 1.2 of Random Eigenvalues by Alan Edelman

#Experiment:  Sample random symmetric Gaussian matrices
#Plot:        Histogram of the eigenvalues
#Theory:      Semicircle as n->infinity
#See also:    RandomMatrices.SemiCircle Distribution

## Parameters
n=1000       # matrix size
t=10         # trials
dx=.2        # binsize
function semicircle_experiment(n,t,dx)
    ## Experiment
    v=Float64[] # eigenvalue samples
    for i = 1:t
        a=randn(n,n)      # n by n matrix of random Gaussians
        s=(a+a')/sqrt(2*n)# symmetrize and normalize matrix
        append!(v,eig(s)[1])# eigenvalues
    end
    grid=-2.4:dx:2.4
    count=hist(v,grid)[2]/(n*t*dx)
    (grid, count)
end
grid,count = semicircle_experiment(n,t,dx)

## Theory
x=[-2:dx:2]
y=sqrt(4-x.^2)/(2*pi)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count, step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, y, "color", "blue", "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "semicircle.png")
end


text:
;; Performance test for JL

(define perf1 (lambda (n) (if n (+ n (perf1 (- n 1))) 0)))
(define perf (lambda (n) (if n (+ (perf1 n) (perf (- n 1))) 0)))

(perf 2048)

max_stars_repo_name: ambotaku/Pico-JL

max_issues_repo_name: ambotaku/Pico-JL

max_forks_repo_name: ambotaku/Pico-JL


text:
using OutOfSampleBootstraps
include("recursive_ols.jl")


text:
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Test

@testset "Basis" begin
    @variables u[1:3]
    @parameters w[1:2]
    h = [u[1]; u[2]; cos(w[1]*u[2]+w[2]*u[3]); u[3]+u[2]]
    h_not_unique = [1u[1]; u[1]; 1u[1]^1; h]
    basis = Basis(h_not_unique, u, parameters = w)

    @test isequal(variables(basis), u)
    @test isequal(parameters(basis), w)
    @test free_parameters(basis) == 6
    @test free_parameters(basis, operations = [+, cos]) == 7
    @test DataDrivenDiffEq.count_operation((ModelingToolkit.Constant(1) + cos(u[2])*sin(u[1]))^3, [+, cos, ^, *]) == 4

    basis_2 = unique(basis)
    @test isequal(basis, basis_2)
    @test size(basis) == size(h)
    @test basis([1.0; 2.0; π], p = [0. 1.]) ≈ [1.0; 2.0; -1.0; π+2.0]
    @test size(basis) == size(basis_2)
    push!(basis_2, sin(u[2]))
    @test size(basis_2)[1] == length(h)+1

    basis_3 = merge(basis, basis_2)
    @test size(basis_3) == size(basis_2)
    @test isequal(variables(basis_3), variables(basis_2))
    @test isequal(parameters(basis_3), parameters(basis_2))

    merge!(basis_3, basis)
    @test basis_3 == basis_2

    push!(basis, 1u[3]+1u[2])
    unique!(basis)
    @test size(basis) == size(h)

    @variables a
    g = [u[1]; u[3]; a]
    basis = Basis(g, [u; a])
    @test basis([1; 2; 3; 4]) == [1; 3; 4]
    g = [1.0*u[1]; 1.0*u[3]; 1.0*u[2]]
    basis = Basis(g, u, parameters = [])
    X = ones(Float64, 3, 10)
    X[1, :] .= 3*X[1, :]
    X[3, :] .= 5*X[3, :]
    # Check the array evaluation
    @test basis(X) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0] * X
    f = jacobian(basis)
    @test f([1;1;1], [0.0 0.0]) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0]
    @test_nowarn sys = ODESystem(basis)
    @test_nowarn [xi for xi in basis]
    @test_nowarn basis[2:end]; basis[2]; first(basis); last(basis); basis[:]

    @variables u[1:2] t

    g = [u[2]; -sin(u[1])*exp(-t)]
    basis = Basis(g, [u...; t])
    @test_nowarn ODESystem(basis, t)
end

@testset "DMD" begin
    # Create some linear data
    A = [0.9 -0.2; 0.0 0.2]
    y = [[10.; -10.]]
    for i in 1:10
        push!(y, A*y[end])
    end
    X = hcat(y...)
    @test_throws AssertionError ExactDMD(X[:, 1:end-2], dt = -1.0)
    estimator = ExactDMD(X[:,1:end-2])
    @test isstable(estimator)
    @test estimator.Ã ≈ A
    @test eigvals(estimator) ≈ eigvals(A)
    @test eigvecs(estimator) ≈ eigvecs(A)
    @test_nowarn dynamics(estimator)
    @test_throws AssertionError dynamics(estimator, discrete = false)
    @test_nowarn update!(estimator, X[:, end-1], X[:,end])
end

@testset "EDMD" begin
    # Test for linear system
    function linear_sys(u, p, t)
        x = -0.9*u[1]
        y = -0.3*u[2]
        return [x;y]
    end

    u0 = [π; 1.0]
    tspan = (0.0, 20.0)
    prob = DiscreteProblem(linear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())

    @variables u[1:2]
    h = [1u[1]; 1u[2]; sin(u[1]); cos(u[1]); u[1]*u[2]]
    basis = Basis(h, u)

    @test_throws AssertionError ExtendedDMD(sol[:,:], basis, dt = -1.0)
    estimator = ExtendedDMD(sol[:,:], basis)
    @test basis == estimator.basis
    basis_2 = reduce_basis(estimator, threshold = 1e-5)
    @test size(basis_2)[1] < size(basis)[1]
    estimator_2 = ExtendedDMD(sol[:,:], basis_2)
    p1 = DiscreteProblem(dynamics(estimator), u0, tspan, [])
    s1 = solve(p1,FunctionMap())
    p2 = DiscreteProblem(dynamics(estimator_2), u0, tspan, [])
    s2 = solve(p2,FunctionMap())
    p3 = DiscreteProblem(linear_dynamics(estimator_2), estimator_2(u0), tspan, [])
    s3 = solve(p3,FunctionMap())
    @test sol[:,:] ≈ s1[:,:]
    @test sol[:,:] ≈ s2[:,:]
    @test estimator_2.basis(sol[:,:])≈ s3[:,:]
    @test eigvals(estimator_2) ≈ [-0.9; -0.3]

    # Test for nonlinear system
    function nonlinear_sys(du, u, p, t)
        du[1] = sin(u[1])
        du[2] = -0.3*u[2] -0.9*u[1]
    end

    prob = DiscreteProblem(nonlinear_sys, u0, tspan)
    sol = solve(prob,FunctionMap())
    estimator = ExtendedDMD(sol[:,:], basis)
    p4 = DiscreteProblem(dynamics(estimator), u0, tspan, [])
    s4 = solve(p4,FunctionMap())
    @test sol[:,:] ≈ s4[:,:]
end

@testset "DMDc" begin
    # Define measurements from unstable system with known control input
    X = [4 2 1 0.5 0.25; 7 0.7 0.07 0.007 0.0007]
    U = [-4 -2 -1 -0.5]
    B = Float32[1; 0]

    # But with a little more knowledge
    sys = DMDc(X, U, B = B)
    @test isa(get_dynamics(sys), ExactDMD)
    @test sys.koopman.Ã ≈[1.5 0; 0 0.1]
    @test get_input_map(sys) ≈ [1.0; 0.0]
    @test !isstable(sys)
    @test_nowarn eigen(sys)

    # Check the solution of an unforced and forced system against each other
    dudt_ = dynamics(sys)
    prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
    sol_unforced = solve(prob,  FunctionMap())

    dudt_ = dynamics(sys, control = (u, p, t) -> -0.5u[1])
    prob = DiscreteProblem(dudt_, X[:, 1], (0., 10.))
    sol = solve(prob, FunctionMap())

    @test all(abs.(diff(sol[1,:])) .< 1e-5)
    @test sol[2,:] ≈ sol_unforced[2,:]
end

@testset "SInDy" begin

    # Create a nonlinear pendulum
    function pendulum(u, p, t)
        x = u[2]
        y = -9.81sin(u[1]) - 0.1u[2]^3 -0.2*cos(u[1])
        return [x;y]
    end

    u0 = [0.99π; -1.0]
    tspan = (0.0, 20.0)
    dt = 0.3
    prob = ODEProblem(pendulum, u0, tspan)
    sol = solve(prob, Tsit5(), saveat = dt)

    # Create the differential data
    DX = similar(sol[:,:])
    for (i, xi) in enumerate(eachcol(sol[:,:]))
        DX[:,i] = pendulum(xi, [], 0.0)
    end
    # Create a basis
    @variables u[1:2]

    # Lots of polynomials
    polys = Operation[1]
    for i ∈ 1:5
        push!(polys, u.^i...)
        for j ∈ 1:i-1
            push!(polys, u[1]^i*u[2]^j)
        end
    end

    # And some other stuff
    h = [cos(u[1]); sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys...]

    basis = Basis(h, u)

    opt = STRRidge(1e-2)
    basis = Basis(h, u, parameters = [])
    Ψ = SInDy(sol[:,:], DX, basis, opt = opt, maxiter = 2000)
    @test_nowarn set_threshold!(opt, 1e-2)
    @test size(Ψ)[1] == 2

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, [])
    sol_ = solve(estimator,Tsit5(), saveat = dt)
    @test sol[:,:] ≈ sol_[:,:]

    opt = ADMM(1e-2, 0.7)
    Ψ = SInDy(sol[:,:], DX, basis, maxiter = 5000, opt = opt)
    @test_nowarn set_threshold!(opt, 1e-2)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, [])
    sol_2 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_2[:,:], 2) < 2e-1

    opt = SR3(1e-2, 1.0)
    Ψ = SInDy(sol[:,:], DX, basis, maxiter = 5000, opt = opt)
    @test_nowarn set_threshold!(opt, 0.1)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, [])
    sol_3 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_3[:,:], 2) < 1e-1

    # Now use the threshold adaptation
    opt = SR3(1e-2, 20.0)
    λs = exp10.(-5:0.1:-1)
    Ψ = SInDy(sol[:,:], DX[:, :], basis, λs,  maxiter = 20, opt = opt)
    estimator = ODEProblem(dynamics(Ψ), u0, tspan, [])
    sol_4 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_4[:,:], 2) < 1e-1

    # Check for errors
    @test_nowarn SInDy(sol[:,:], DX[1,:], basis, λs, maxiter = 1, opt = opt)
    @test_nowarn SInDy(sol[:, :], DX[1, :], basis, λs, maxiter = 1, opt = opt, denoise = true, normalize = true)

    # Check with noise
    X = sol[:, :] + 1e-3*randn(size(sol[:,:])...)
    Ψ = SInDy(X, DX, basis, λs, maxiter = 10000, opt = opt, denoise = true, normalize = true)

    estimator = ODEProblem(dynamics(Ψ), u0, tspan, [])
    sol_4 = solve(estimator,Tsit5(), saveat = dt)
    @test norm(sol[:,:] - sol_4[:,:], 2) < 5e-1

    # Check sparse_regression
    X .= Array(sol)
    opt = SR3(1e-1, 1.0)
    maxiter = 5000
    θ = basis(X)
    Ξ1 = sparse_regression(X, DX, basis, parameters(basis), maxiter, opt, true, true)
    Ξ2 = similar(Ξ1)
    Ξ3 = similar(Ξ1)
    sparse_regression!(Ξ2, X, DX, basis, parameters(basis), maxiter, opt, true, true)
    sparse_regression!(Ξ3, θ, DX, maxiter, opt, true, true)

    @test Ξ1 ≈ Ξ2 ≈ Ξ3
    @test isapprox(norm(Ξ1'*θ - DX,2), 25.18; atol = 1e-1)

end

@testset "ISInDy" begin

    # Create a test problem
    function simple(u, p, t)
        return [(2.0u[2]^2 - 3.0)/(1.0 + u[1]^2); -u[1]^2/(2.0 + u[2]^2); (1-u[2])/(1+u[3]^2)]
    end

    u0 = [2.37; 1.58; -3.10]
    tspan = (0.0, 10.0)
    prob = ODEProblem(simple, u0, tspan)
    sol = solve(prob, Tsit5(), saveat = 0.1)

    # Create the differential data
    X = sol[:,:]
    DX = similar(X)
    for (i, xi) in enumerate(eachcol(X))
        DX[:, i] = simple(xi, [], 0.0)
    end

    # Create a basis
    @variables u[1:3]
    polys = ModelingToolkit.Operation[]
    # Lots of basis functions
    for i ∈ 0:5
        if i == 0
            push!(polys, u[1]^0)
        end
        for ui in u
            if i > 0
                push!(polys, ui^i)
            end
        end
    end

    basis= Basis(polys, u)

    opt = ADM(1e-3)
    Ψ = ISInDy(X, DX, basis, opt = opt, maxiter = 100, rtol = 0.1)

    # Simulate
    estimator = ODEProblem(dynamics(Ψ), u0, tspan)
    sol_ = solve(estimator, Tsit5(), saveat = 0.1)
    @test sol[:,:] ≈ sol_[:,:]
end

@testset "Utilities" begin
    t = collect(-2:0.01:2)
    U = [cos.(t).*exp.(-t.^2) sin.(2*t)]
    S = Diagonal([2.; 3.])
    V = [sin.(t).*exp.(-t) cos.(t)]
    A = U*S*V'
    σ = 0.5
    Â = A + σ*randn(401, 401)
    n_1 = norm(A-Â)
    B = optimal_shrinkage(Â)
    optimal_shrinkage!(Â)
    @test norm(A-Â) < n_1
    @test norm(A-B) == norm(A-Â)

    X = randn(3, 100)
    Y = randn(3, 100)
    k = 3

    @test AIC(k, X, Y) == 2*k-2*log(sum(abs2, X- Y))
    @test AICC(k, X, Y) == AIC(k, X, Y)+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)
    @test BIC(k, X, Y) == -2*log(sum(abs2, X -Y)) + k*log(size(X)[2])
    @test AICC(k, X, Y, likelyhood = (X,Y)->sum(abs, X-Y)) == AIC(k, X, Y, likelyhood = (X,Y)->sum(abs, X-Y))+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)

    # Numerical derivatives
    function lorenz(u,p,t)
        x, y, z = u
        ẋ = 10.0*(y - x)
        ẏ = x*(28.0-z) - y
        ż = x*y - (8/3)*z
        return [ẋ, ẏ, ż]
    end
    u0 = [1.0;0.0;0.0]
    tspan = (0.0,50.0)
    dt = 0.005
    prob = ODEProblem(lorenz,u0,tspan)
    sol = solve(prob, Tsit5(), saveat = dt)

    X = Array(sol)
    DX = similar(X)
    for (i, xi) in enumerate(eachcol(X))
        DX[:,i] = lorenz(xi, [], 0.0)
    end

    windowSize, polyOrder = 9, 4
    halfWindow = Int(ceil((windowSize+1)/2))
    DX_sg = similar(X)
    DX_sg_cropped = similar(X[:,halfWindow+1:end-halfWindow])
    X_cropped = similar(X[:,halfWindow+1:end-halfWindow])
    for i =1:size(X,1)
        DX_sg[i,:] = savitzky_golay(X[i,:], windowSize, polyOrder, deriv=1, dt=dt, crop=false)
        X_cropped[i,:], DX_sg_cropped[i,:] = savitzky_golay(X[i,:], windowSize, polyOrder, deriv=1, dt=dt)
    end
    DX_sg2 = savitzky_golay(X, windowSize, polyOrder, deriv=1, dt=dt, crop=false)
    X_cropped2, DX_sg_cropped2 =  savitzky_golay(X, windowSize, polyOrder, deriv=1, dt=dt)
    @test(DX_sg2 == DX_sg)
    DX_sg = DX_sg[:,halfWindow+1:end-halfWindow]
    @test X_cropped == X[:,halfWindow+1:end-halfWindow]
    @test DX_sg_cropped == DX_sg
    @test X_cropped2 == X[:,halfWindow+1:end-halfWindow]
    @test DX_sg_cropped2 == DX_sg

    DX = DX[:,halfWindow+1:end-halfWindow]
    @test isapprox(DX_sg, DX, rtol=1e-2)

    # Sampling
    X = randn(Float64, 2, 100)
    t = collect(0:0.1:9.99)
    Y = randn(size(X))
    xt = burst_sampling(X, 5, 10)
    @test 10 <= size(xt)[end] <= 60
    @test all([any(xi .≈ X) for xi in eachcol(xt)])
    xt, tt = burst_sampling(X, t, 5, 10)
    @test all(diff(tt) .> 0.0)
    @test size(xt)[end] == size(tt)[end]
    @test all([any(xi .≈ X) for xi in eachcol(xt)])
    @test !all([any(xi .≈ Y) for xi in eachcol(xt)])
    xs, ts = burst_sampling(X, t, 2.0, 1)
    @test all([any(xi .≈ X) for xi in eachcol(xs)])
    @test size(xs)[end] == size(ts)[end]
    @test ts[end]-ts[1] ≈ 2.0
    X2n = subsample(X, 2)
    t2n = subsample(t, 2)
    @test size(X2n)[end] == size(t2n)[end]
    @test size(X2n)[end] == Int(round(size(X)[end]/2))
    @test X2n[:, 1] == X[:, 1]
    @test X2n[:, end] == X[:, end-1]
    @test all([any(xi .≈ X) for xi in eachcol(X2n)])
    xs, ts = subsample(X, t, 0.5)
    @test size(xs)[end] == size(ts)[end]
    @test size(xs)[1] == size(X)[1]
    @test all(diff(ts) .≈ 0.5)
    # Loop this a few times to be sure its right
    @test_nowarn for i in 1:20
        xs, ts = burst_sampling(X, t, 2.0, 1)
        xs, ts = subsample(X, t, 0.5)
    end
  end


text:
module RoME

using
  IncrementalInference,
  Graphs,
  TransformUtils,
  CoordinateTransformations,
  Rotations,
  KernelDensityEstimate,
  Distributions,
  JLD,
  HDF5,
  ProgressMeter,
  DocStringExtensions,
  Compat

import Base: +, \, convert
import TransformUtils: ⊖, ⊕, convert, compare, ominus, veeQuaternion
import IncrementalInference: convert, getSample, reshapeVec2Mat  #, compare

export
  # pass throughs from TransformUtils
  SE2,
  se2vee,
  se2vee!,
  SE3,
  Euler,
  Quaternion,
  AngleAxis,
  SO3,
  so3,
  compare,
  convert,

  # pass throughs from IncrementalInference
  FunctorSingleton,
  FunctorPairwise,
  FunctorPairwiseNH,
  FunctorSingletonNH,
  ls,
  addFactor!,
  addNode!,
  getVert,
  getVertKDE,
  getVal,
  setVal!,
  getData,
  FNDencode,
  FNDdecode,
  localProduct,
  predictbelief,
  VNDencoder,
  VNDdecoder,
  GenericWrapParam,
  wipeBuildNewTree!,
  inferOverTree!,
  inferOverTreeR!,
  writeGraphPdf,
  savejld,
  loadjld,
  FactorGraph,
  initializeNode!,
  isInitialized,
  ensureAllInitialized!,
  # overloaded functions from IIF
  # decodefg,
  # convertfrompackedfunctionnode,

  # RoME specific functions
  initfg,
  measureMeanDist,
  predictBodyBR,
  getLastPose,
  getLastPose2D,
  odomKDE,
  initFactorGraph!,
  addOdoFG!,
  addposeFG!,
  newLandm!,
  addBRFG!,
  addMMBRFG!,
  addAutoLandmBR!,
  projNewLandm!,
  malahanobisBR,
  veePose3,
  veePose,
  \,
  RangeAzimuthElevation,

  # types
  BetweenPoses,

  # helper functions
  get2DSamples,
  getAll2D,
  get2DSampleMeans,
  getAll2DMeans,
  getAll2DPoses,
  get2DPoseSamples,
  get2DPoseMeans,
  getKDE,
  getVertKDE,
  get2DPoseMax,
  getAll2DLandmarks,
  get2DLandmSamples,
  get2DLandmMeans,
  get2DLandmMax,

  # helper functions
  getLastLandm2D,
  getLastPose2D,
  getNextLbl,

  # RobotUtils
  getRangeKDEMax2D,

  # some transform functions
  cart2pol,
  pol2cart,

  # Feature tracking code
  Feature,
  initTrackersFrom,
  propAllTrackers!,
  measUpdateTrackers!,
  assocMeasWFeats!,

  lsrBR,

  # solve with isam in pytslam
  doISAMSolve,
  drawCompPosesLandm,

  # Victoria Park data specific
  LaserFeatures,
  addLandmarksFactoGraph!,
  appendFactorGraph!,
  doBatchRun,
  rotateFeatsToWorld,

  # Didson model
  evalPotential,
  LinearRangeBearingElevation,
  project!,
  project,
  backprojectRandomized!,
  residual!,
  residualLRBE!,
  reuseLBRA,
  ominus,
  ominus!,
  +,
  evalPotential,
  getSample!,
  getSample,
  # obsolete
  WrapParam,
  WrapParamArray,

  # Didson convenience function
  addLinearArrayConstraint,

  # camera model -- TODO --separate out
  CameraIntrinsic,
  CameraExtrinsic,
  CameraModelFull,
  project!,
  project,
  backprojectRandomized!,
  # keep
  cameraResidual!,

  # Point2D
  Point2DPoint2D,
  PackedPoint2DPoint2D,
  Pose2DPoint2DBearingRange,
  PackedPose2DPoint2DBearingRange,
  Pose2DPoint2DBearing,
  Pose2DPoint2DRange,
  Point2DPoint2DRange,
  PackedPoint2DPoint2DRange,
  PriorPoint2D,
  PackedPriorPoint2D,
  solveLandm,
  solvePose2,
  solveSetSeps,
  addPose2Pose2!,

  # Point2D with null hypotheses
  PriorPoint2DensityNH,
  PackedPriorPoint2DensityNH,

  # acoustics
  Pose2DPoint2DBearingRangeDensity,
  PackedPose2DPoint2DBearingRangeDensity,
  Pose2DPoint2DRangeDensity,
  PackedPose2DPoint2DRangeDensity,

  Pose2,
  Point2,
  Pose3,
  Point3,
  Prior,
  Pose2Pose2_NEW,

  # Pose2D
  PriorPose2,
  PackedPriorPose2,
  Pose2Pose2,
  PackedPose2Pose2,
  addPose2Pose2,

  # packedTypes
  passTypeThrough,

  # Pose3D
  PriorPose3,
  PackedPriorPose3,
  Pose3Pose3,
  PackedPose3Pose3,
  projectParticles,
  ⊕,
  Pose3Pose3NH,
  PackedPose3Pose3NH,

  # partial Pose3
  PartialPriorRollPitchZ,
  PackedPartialPriorRollPitchZ,
  PartialPose3XYYaw,
  PackedPartialPose3XYYaw,
  PartialPose3XYYawNH,
  PackedPartialPose3XYYawNH,

  # MultipleFeatures2D constraint functions
  MultipleFeatures2D,
  getUvecScaleFeature2D,
  getUvecScaleBaseline2D,

  # SLAM specific functions
  SLAMWrapper,

  # FG Analysis tools
  rangeErrMaxPoint2,
  rangeCompAllPoses,
  rangeCompAllPoses,

  # new robot navigation functionality
  triggerPose,
  GenericInSituSystem,
  InSituSystem,
  makeInSituSys,
  makeGenericInSituSys,
  advOdoByRules,
  poseTrigAndAdd!,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  addSoftEqualityPoint2D,
  vectoarr2

  ## Visualization tools have been moved to RoMEPlotting.jl
  # draw pose beliefs etc
  # Some vizualization tools
  # togglePrtStbLines,
  # plotLsrScanFeats,
  # drawFeatTrackers,
  # saveImgSeq,
  # drawPoses,
  # drawLandms,
  # drawPosesLandms,
  # drawSubmaps,
  # investigatePoseKDE,
  # plotPose3Pairs,
  # drawMarginalContour,
  # accumulateMarginalContours,
  # progressExamplePlot,
  # plotTrckStep,

@compat abstract type BetweenPoses <: IncrementalInference.FunctorPairwise end

@compat const VoidUnion{T} = Union{Void, T}

@compat const CTs = CoordinateTransformations
@compat const TUs = TransformUtils

vectoarr2(v) = reshape(v, length(v),1)

include("BayesTracker.jl")

include("SensorModels.jl")
include("CameraModel.jl")
include("Point2D.jl")
include("Pose2D.jl")
include("Pose3D.jl")
include("BearingRange2D.jl")

include("BearingRangeDensity2D.jl")

include("Pose3Pose3.jl")
include("PartialPose3.jl")
include("MultipleFeaturesConstraint.jl")

include("InertialPose3.jl")

include("Slam.jl")

# include("RobotViz.jl")
include("RobotUtils.jl")

include("SimulationUtils.jl")

include("FactorGraphAnalysisTools.jl")

include("WheeledRobotUtils.jl")
include("NavigationSystem.jl")

# include("dev/ISAMRemoteSolve.jl")

end


text:
using Base.Test
using GARCH

filename = Pkg.dir("GARCH", "test","data","SPY.csv")
close = readcsv(filename,Float64)[:,2]
ret = diff(log(close))
ret = ret - mean(ret)
fit = garchFit(ret)
param = [2.469347e-06, 1.142268e-01, 8.691734e-01] #R fGarch garch(1,1) estimated params
@test_approx_eq_eps(fit.params,param,1e-3)


text:
import PyCall, SymPy

using MLStyle

stats = PyCall.pyimport_conda("sympy.stats", "sympy")
SymPy.import_from(stats)

sym(x) = SymPy.symbols(x)

macro ℓ(expr)
    args = @match expr begin
        Expr(head, args...) => args
    end
    d = args[1]
    ps = args[2:end]
    quote
        L = SymPy.density(stats.$d(:foo, sym.($ps)...)).pdf(sym(:x))
        SymPy.sympy.expand_log(log(L), force=true)
    end
end

@ℓ Cauchy(μ, σ)


text:
using ITensors,
      Test

@testset "TagSet" begin
  ts = TagSet("t3,t2,t1")
  ts2 = copy(ts)
  @test ts == ts2
  @test hastags(ts,"t1")
  @test hastags(ts,"t2")
  @test hastags(ts,"t3")
  @test hastags(ts,"t3,t1")
  @test !hastags(ts,"t4")
  @test TagSet(ts) === ts

  t1 = TagSet("t1")
  t2 = TagSet("t2")
  t3 = TagSet("t3")
  @test ts[1] == t1[1]
  @test ts[2] == t2[1]
  @test ts[3] == t3[1]

  @testset "Empty TagSet" begin
    ts1 = TagSet()
    @test length(ts1) == 0

    ts2 = TagSet("")
    @test ts2 == ts1
    @test length(ts2) == 0
  end

  @testset "Ignore Whitespace" begin
    ts = TagSet(" aaa , bb bb  , ccc    ")
    @test hastags(ts," aaa ")
    @test hastags(ts,"aaa")
    @test hastags(ts," aa a ")
    @test hastags(ts,"bbbb")
  end

  @testset "Remove tags" begin
    ts1 = TagSet("x,y,z")
    ts2 = TagSet("x,z")
    @test removetags(ts1,"y") == ts2
  end

  @testset "Unicode tags" begin
    ts = TagSet("α")
    @test length(ts) == 1
    @test hastags(ts, "α")
    @test ts[1] == ITensors.SmallString("α")

    ts = TagSet("α,β")
    @test length(ts) == 2
    @test hastags(ts, "β")
    @test hastags(ts, "α")
    @test ts[1] == ITensors.SmallString("α")
    @test ts[2] == ITensors.SmallString("β")

    ts = TagSet("αβγδϵζηθ,ijklmnop,qrstuvwx,ΑΒΓΔΕΖΗΘ")
    @test length(ts) == 4
    @test hastags(ts, "αβγδϵζηθ")
    @test hastags(ts, "ijklmnop")
    @test hastags(ts, "qrstuvwx")
    @test hastags(ts, "ΑΒΓΔΕΖΗΘ")
    @test ts[1] == ITensors.SmallString("ijklmnop")
    @test ts[2] == ITensors.SmallString("qrstuvwx")
    @test ts[3] == ITensors.SmallString("ΑΒΓΔΕΖΗΘ")
    @test ts[4] == ITensors.SmallString("αβγδϵζηθ")
  end

  @testset "Tag too long" begin
    @test_throws ErrorException TagSet("ijklmnopq")
    @test_throws ErrorException TagSet("abcd,ijklmnopq")
    @test_throws ErrorException TagSet("ijklmnopqr,abcd")
  end

  @testset "Check for Integer Tags" begin
    @test_throws ErrorException TagSet("123")
  end

  @testset "Show TagSet" begin
    ts = TagSet("Site,n=2")
    @test length(sprint(show,ts)) > 1
  end

  @testset "addtags" begin
    ts = TagSet("Blue")
    @test hastags(ts,"Blue")

    ts = addtags(ts,"Red")
    @test hastags(ts,"Blue")
    @test hastags(ts,"Red")

    ts = addtags(ts,"Green")
    @test hastags(ts,"Blue")
    @test hastags(ts,"Red")
    @test hastags(ts,"Green")

    ts = addtags(ts,"Yellow")
    @test hastags(ts,"Blue")
    @test hastags(ts,"Red")
    @test hastags(ts,"Green")
    @test hastags(ts,"Yellow")

    @test_throws ErrorException addtags(ts,"Orange")
  end
end

nothing


text:
module EscapeSimulator
##
using Random
using Distributions
using SpecialFunctions
using DataFrames
using StatsBase
using BitBasis
using HypergeometricFunctions

include("count_likelihood.jl") # analyzes count data using mle and bayesian methods
include("sensitivity_analysis.jl") # Minimum disparity functions
include("WrightSampler.jl")
include("virus_gillespie.jl")
include("virus_trace.jl")
include("equilibrium.jl")

end # module

max_stars_repo_name: StatPhysBio/EscapeSimulator

max_issues_repo_name: StatPhysBio/EscapeSimulator

max_forks_repo_name: StatPhysBio/EscapeSimulator


text:
a = 1
b = 5
println(a  b)


text:
function histogram_name(prefix::String, ppi::Int64, N::Int64, K2::Int64, K3::Int64, rmax::Float64, max_pictures::Int64, suffix::String, gamma::Float64=0.0, sigma::Float64=1.0, lambda::Float64=0.0)
    return "$(prefix)$(ppi)p_N$(N)_K2_$(K2)_K3_$(K3)_R$(rmax)_P$(max_pictures)$(gamma > 0.0 ? "_G$(gamma)_S$(sigma)" : "")_lambda_$(lambda)$(suffix)"
end

"""Starts a cluster job for synthetic correlation generation"""
function generate_histograms(; max_pictures::Integer=Integer(0), N::Integer=32, number_incident_photons::Integer=500, incident_photon_variance::Integer = 0, gamma::Float64=0.0, sigma::Float64=1.0, noise_photons::Int64=0, Ncores::Integer=8, batchsize::Integer = Integer(1e4), successive_jobs::Integer=1, prefix::String="correlations_", suffix::String="", use_cube::Bool=true, qcut_ratio::Float64=1.0, K2::Int64=38, K3::Int64=26, rmax::Float64=35.0, histogram_method="histogramCorrelationsInPicture_alltoall", structure_pdb_path::String="", number_particles::Int64=1, lambda::Float64=0.0, beamstop_width::Float64=0.0, beamstop_only::Bool=false)
    name = histogram_name(prefix, calculate_ppi(number_incident_photons), N, K2, K3, rmax, max_pictures, suffix, gamma, sigma, lambda)

    julia_script = """
    using ThreePhotons

    volume = 0.0

    if $(use_cube)
        # _,_,volume = createCubicStructure("$(structure_pdb_path)", 4*$K2+1, 2.0*$rmax)
        volume = loadCube("$(ENV["DETERMINATION_DATA"])/histograms/synthetic_crambin/intensityCube_high.mrc")
    else
        _,_,intensity = createSphericalHarmonicsStructure("$(structure_pdb_path)", 35, $K2, $rmax)
        volume = getSurfaceVolume(intensity)
        volume.radial_interp = false
    end

    generateHistogram(volume; qcut=$(qcut_ratio)*volume.rmax, K2=$K2, K3=$K3, N=$N, max_pictures=$max_pictures, number_incident_photons=$number_incident_photons, incident_photon_variance=$incident_photon_variance, numprocesses=$(Ncores), file="histo.dat", noise=GaussianNoise($gamma, $sigma, $(noise_photons)), batchsize = $batchsize, histogramMethod=$histogram_method, number_particles=$(number_particles), lambda=$(lambda), beamstop_width=$(beamstop_width), beamstop_only=$(beamstop_only))
    """

    launch_job("data_generation/$name", Ncores, false, julia_script, successive_jobs, architecture="")
end

function generate_histogram_image(img::Int64, ppi::Int64, K2::Int64, K3::Int64, N::Int64, lambda::Float64=0.0; setsize::Int64=Integer(2*2.048e7), name::String="", beamstop_width::Float64=0.0)
    numbersets = ceil(Int64, img / setsize)
    for i = 1:numbersets
        generate_histograms(; max_pictures = setsize, Ncores=8, N=N, number_incident_photons=calculate_incident_photons(ppi), batchsize = Integer(setsize/8), successive_jobs=1, prefix="parts/$(name)SH_", suffix="_$(i)", use_cube=false, qcut_ratio=1.0, K2=K2, K3=K3, rmax=float(K2), histogram_method="histogramCorrelationsInPicture_alltoall", structure_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", lambda=lambda, beamstop_width=beamstop_width)
    end
end

function generate_histogram_set_ppi(ppi::Int64=10; K2::Int64=38, K3::Int64=26, N::Int64=32, lambda::Float64=0.0)
    for img in calculate_images_ppi(ppi)
        generate_histogram_image(img, ppi, K2, K3, N, lambda)
    end
end

# generate_histogram_set_ppi(10, 38, 26, 32)
# generate_histogram_set_ppi(25)
# generate_histogram_set_ppi(50)

function generate_noisy_histograms()
    noise_photons = Dict(0.5=>14, 0.75=>14, 1.125=>14, 2.5=>25)
    let K2 = 38, K3 = 26, N = 32
        # for sigma in [0.5, 0.75, 1.125]
        let sigma = 2.5
            for gamma in collect(0.1:0.1:0.5)
                setsize = Integer(2*2.048e7)
                numbersets = Integer(3.2768e9 / setsize)
                for i = 1:numbersets
                    generate_histograms(; max_pictures = setsize, Ncores=8, N=N, number_incident_photons=calculate_incident_photons(10), batchsize = Integer(2.56e5), successive_jobs=3, prefix="parts/SH_", suffix="_$(i)", use_cube=false, qcut_ratio=1.0, K2=K2, K3=K3, rmax=float(K2), histogram_method="histogramCorrelationsInPicture_alltoall", gamma=gamma, sigma=sigma, noise_photons=round(Int64, noise_photons[sigma] * gamma), structure_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb")
                end
            end
        end
    end
end

# generate_noisy_histograms()

function generate_single_multiparticle_histogram(number_images::Int64, setsize::Int64, number_particles::Int64=2, ppi::Int64=10; K2::Int64=38, K3::Int64=26, N::Int64=32, lambda::Float64=0.0)
    numbersets = ceil(Int64, number_images / setsize)
    for i = 1:numbersets
        generate_histograms(; max_pictures = setsize, Ncores=8, N=N, number_incident_photons=calculate_incident_photons(ppi), batchsize = Integer(setsize/8), successive_jobs=1, prefix="parts/multi_$(number_particles)_SH_", suffix="_$(i)", use_cube=false, qcut_ratio=1.0, K2=K2, K3=K3, rmax=float(K2), histogram_method="histogramCorrelationsInPicture_alltoall", structure_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", number_particles=number_particles, lambda=lambda)
    end
end

using HDF5
"""Distributes the calculation of correlations among many jobs"""
function run_calculate_correlation_from_images(particle_name::String, images_path::String, images_per_job::Int64, K2::Int64, K3::Int64, N::Int64; Ncores::Int64=8, symmetrize::Bool=false, sample_photons::Int64=Integer(1e6))

    file = h5open(images_path, "r")
    photonConverter = read(file["photonConverter"])
    _,_,images_in_file = Base.size(photonConverter["pnccdBack"]["photonCount"])
    number_runs = floor(Int64, images_in_file / images_per_job)

    for n in 1:number_runs
        julia_script = """
        using HDF5
        using ThreePhotons
        using Images

        K2 = $K2
        K3 = $K3
        N = $N

        photonConverter = h5open("$(images_path)", "r") do file
            read(file["photonConverter"])
        end

        #With resize
        image_list = [ Images.imresize(convert(Array{Float64},photonConverter["pnccdBack"]["photonCount"][:,:,i]), (2*K2+1, 2*K2+1)) for i=$((n-1)*images_per_job+1):$(n*images_per_job)]

        # image_list = [ convert(Array{Float64},photonConverter["pnccdBack"]["photonCount"][:,:,i]) for i=$((n-1)*images_per_job+1):$(n*images_per_job)]

        #Calculate maximum intensity integral for scaling
        overall_maximum = maximum(Float64[sum(abs, image_list[i]) for i = 1:length(image_list)])

        calculate_correlations_in_image_using_single_photons(image_list, K2, K3, N, "histo.dat", overall_maximum, Integer($(sample_photons)), $(symmetrize))
        """
        launch_job("exp_data/parts/$(particle_name)_$(n)", Ncores, false, julia_script, 1)#, memory="$(Ncores*1.5)G")
    end
end

"""Distributes the calculation of correlations among many jobs"""
function run_calculate_beamstop_correlation(jobname::String, images_path::String, K2::Int64, K3::Int64, N::Int64; Ncores::Int64=8)

    julia_script = """
    using HDF5
    using ThreePhotons
    using Images

    K2 = $K2
    K3 = $K3
    N = $N

    file = h5open("$(images_path)", "r")
    photonConverter = read(file["photonConverter"])

    sum = 0
    for i = 1:500
        sum += convert(Array{Float64},photonConverter["pnccdBack"]["photonCount"][:,:,i])
    end

    # beamstop_map = (convert(Array{Float64},(abs(sum) .< eps())))

    #With resize
    beamstop_map = Images.imresize(convert(Images.Image,beamstop_map), (2*K2+1, 2*K2+1)).data

    calculate_correlations_in_image_using_single_photons([1.0 - beamstop_map], K2, K3, N, "histo.dat", 4.0, Integer(2e5), false)
    """
    launch_job("exp_data/$(jobname)_K2_$(K2)_K3_$(K3)_N$(N)", Ncores, false, julia_script, 1)
end

# """Distributes the calculation of correlations among many jobs"""
# function run_calculate_generic_beamstop_correlation(jobname::String, K2::Int64, K3::Int64, N::Int64; Ncores::Int64=8)
#
#     julia_script = """
#     using ThreePhotons
#     using Images
#
#     K2 = $K2
#     K3 = $K3
#     N = $N
#
#     function beamstop(K::Vector{Float64}, beamstop_width::Float64)
#         (abs(K[1]) > beamstop_width && abs(K[2]) > beamstop_width) ? 1.0 : 0.0
#     end
#
#     qm = 3.141592653589793
#     delq = 0.08267349088394192
#     range = -qm:delq:qm
#     beamstop_width = qmax(K2,float(K2))/20.0
#     crambin_beamstop = [beamstop([k1,k2], beamstop_width) for k1 in range, k2 in range]
#
#     calculate_correlations_in_image([crambin_beamstop], K2, K3, N)
#     """
#     launch_job("exp_data/$(jobname)", Ncores, false, julia_script, 1)
# end

#Split datasets into approx. 7500 images each
exp_filelist1 = String[
"amo86615_186_PR772_single.h5",
"amo86615_188_PR772_single.h5",
"amo86615_190_PR772_single.h5",
"amo86615_191_PR772_single.h5",
"amo86615_192_PR772_single.h5"
]

exp_filelist2 = String[
"amo86615_193_PR772_single.h5",
"amo86615_194_PR772_single.h5",
"amo86615_196_PR772_single.h5",
"amo86615_197_PR772_single.h5"
]

exp_filelist = String[
"amo86615_186_PR772_single.h5",
"amo86615_188_PR772_single.h5",
"amo86615_190_PR772_single.h5",
"amo86615_191_PR772_single.h5",
"amo86615_192_PR772_single.h5",
"amo86615_193_PR772_single.h5",
"amo86615_194_PR772_single.h5",
"amo86615_196_PR772_single.h5",
"amo86615_197_PR772_single.h5"
]

# for file in exp_filelist
#     run_calculate_correlation_from_images("coliphage_single_photons_K2_40_K3_30_N32/$file", environment_path("exp_data/Coliphage_PR772/$file"), 100, 40, 30, 32, symmetrize=false)
# end

#Calculate beamstop of Coliphage_PR772
# run_calculate_beamstop_correlation("coliphage_beamstop", environment_path("exp_data/Coliphage_PR772/amo86615_186_PR772_single.h5"), 40, 30, 32)

#Calculate beamstop of crambin
# run_calculate_generic_beamstop_correlation("crambin_beamstop", 38, 26, 32)

function combine_histograms(dir::String, num::Int64)
    println("Processing $dir")

    println("\tLoading file #1")
    params,c2_full,c3_full,c1_full = deserializeFromFile("$(dir)_1/histo.dat")
    part_images = params["num_pictures"]

    for i = 2:num
        println("\tLoading file #$i")
        histofile = "$(dir)_$(i)/histo.dat"
        if isfile(histofile)
            params_part, c2_part, c3_part, c1_part =  deserializeFromFile(histofile)
            part_triplets = sum(c3_part)
            if part_triplets < float(1.0e20)
                params["num_pictures"] += params_part["num_pictures"]
                c1_full += c1_part
                c2_full += c2_part
                c3_full += c3_part
                println("\tAdded $(part_triplets) triplets and $(params_part["num_pictures"]) images.")
            else
                println("\t------ERROR: $i file.")
            end
        end
    end

    finalpath = replace(dir, "parts/", "")
    finalpath = replace(finalpath, "P$(part_images)", "P$(params["num_pictures"])")

    println("Writing to $(finalpath)/histo.dat")
    try mkdir(finalpath) end
    serializeToFile("$(finalpath)/histo.dat", (params, c2_full, c3_full, c1_full))
end

function combine_histograms(filelist::Array{String}, output_dir::String)

    println("Processing $(filelist[1]).")
    params,c2_full,c3_full,c1_full = deserializeFromFile("$(filelist[1])")
    part_images = params["num_pictures"]

    for i = 2:length(filelist)
        histofile = filelist[i]
        println("Processing $histofile.")
        if isfile("$histofile")
            try
                params_part, c2_part, c3_part, c1_part =  deserializeFromFile("$histofile")
                part_triplets = sum(c3_part)
                if part_triplets < float(1.0e20)
                    params["num_pictures"] += params_part["num_pictures"]
                    c1_full += c1_part
                    c2_full += c2_part
                    c3_full += c3_part
                    println("\tAdded $(part_triplets) triplets and $(params_part["num_pictures"]) images.")
                else
                    println("\t------ERROR: $histofile.")
                end
            catch err
                println("\t------ERROR: $histofile.\n$err")
            end
        end
    end

    println("Writing to $(output_dir)/histo.dat")
    try mkdir(output_dir) end
    serializeToFile("$(output_dir)/histo.dat", (params, c2_full, c3_full, c1_full))
end

function create_exp_filelist(name::String="coliphage")
    root = environment_path("exp_data/parts/$(name)")
    list = readdir(root)
    map((p)->"$root/$p/histo.dat", list)
end

function process_exp_data(name::String="coliphage", beamstop::String="coliphage_beamstop_K2_38_K3_30_N_32", Gauss_filter::Float64=0.0, symmetrize::Bool=false)
    combine_histograms(create_exp_filelist(name), environment_path("exp_data/$(name)"))
    p,c2,c3,c1 = deserializeFromFile(environment_path("exp_data/$(name)/histo.dat"))
    if beamstop != ""
        _,c2_beamstop,c3_beamstop,_ = deserializeFromFile(environment_path("exp_data/$(beamstop)/histo.dat"))
        c2,c3 = correct_beamstop(c2, c3, c2_beamstop, c3_beamstop)
    end
    if symmetrize
        c2 = symmetrize_correlation(c2)
        c3 = symmetrize_correlation(c3)
    end
    if Gauss_filter > 0.0
        c2 = add_Gaussian_filter(c2, Gauss_filter)
        c3 = add_Gaussian_filter(c3, Gauss_filter)
    end
    newname = "$(name)_$(beamstop == "" ? "nobsc" : "bsc")_processed$(Gauss_filter > 0.0 ? "_smoothed" : "_notsmoothed")$(symmetrize ? "_symmetrized" : "")"
    try mkdir(environment_path("exp_data/$(newname)")) end
    serializeToFile(environment_path("exp_data/$(newname)/histo.dat"), (p, c2, c3, c1))
end

function combine_set_noise(img::Int64, setsize::Int64, sigmavals::Vector{Float64}=Float64[0.5, 0.75, 1.125], gammavals::Vector{Float64}=[0.1, 0.2, 0.3, 0.4, 0.5], ppi::Int64=10, K::Int64=38, N::Int64=32, lambda::Float64=0.0)
    for sigma in sigmavals
        for gamma in gammavals
            name = histogram_name("parallel/data_generation/parts/SH_", ppi, N, K, float(K), setsize, "", gamma, sigma, lambda)
            combine_histograms(name, img, K, N, setsize)
        end
    end
end

#Combine noisy histograms
# combine_set_noise(Integer(3.2768e9), Integer(2*2.048e7), [2.5], [0.4, 0.5], 10, 38, 32)

function combine_set(images::Array{Int64}, setsize::Int64, ppi::Int64=10, K::Int64=38, N::Int64=32, lambda::Float64=0.0, prefix::String="$(ENV["DETERMINATION_DATA"])/data_generation/parts/SH_")
    for pic in images
        name = histogram_name(prefix, ppi, N, K, float(K), setsize, "", 0.0, 0.0, lambda)
        combine_histograms(name, pic, K, N, setsize)
    end
end

# generate_single_multiparticle_histogram(Integer(3.2768e9), Integer(2.048e7), 2)
# combine_histograms(environment_path("data_generation/parts/multi_2_SH_10p_N32_K2_38_K3_26_R38.0_P40960000"), 80)

# generate_histogram_image(Integer(3.2768e9), 10, 38, 26, 32, 2.5; setsize=Integer(2*2.048e7), name="Ewald_large_")
# generate_histogram_image(Integer(3.2768e8), 10, 38, 26, 32, 2.5; setsize=Integer(2*2.048e7), name="Ewald_medium_")
# combine_histograms(environment_path("data_generation/parts/Ewald_lambda_2.0_SH_10p_N32_K2_38_K3_26_R38.0_P40960000"), 80)

#With beamstop
# generate_histogram_image(Integer(3.2768e9), 10, 38, 26, 32; setsize=Integer(2*2.048e7), name="Ewald_lambda_0.0_beamstop_", lambda=0.0, beamstop_width=qmax(38,38.0)/20.0)
# combine_histograms(environment_path("data_generation/parts/Ewald_lambda_0.0_beamstop_SH_10p_N32_K2_38_K3_26_R38.0_P4096000"), 80)

#beamstop only
# generate_histograms(; max_pictures = Integer(2*2.048e7), Ncores=8, N=32, number_incident_photons=100, batchsize = Integer(2.048e6), successive_jobs=1, prefix="beamstop_only_3_", suffix="", use_cube=false, qcut_ratio=1.0, K2=38, K3=26, rmax=float(38), histogram_method="histogramCorrelationsInPicture_alltoall", structure_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", lambda=0.0, beamstop_width=qmax(38,38.0)/20.0, beamstop_only=true)

#Coliphage
# combine_histograms(environment_path("exp_data/parts/coliphage_symmetric"), 42)