text:
##################################################################
"""
$(TYPEDEF)

Extendable sparse matrix. A nonzero  entry of this matrix is contained
either in cscmatrix, or in extmatrix, never in both.

$(TYPEDFIELDS)
"""
mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    """
    Final matrix data
    """
    cscmatrix::SparseMatrixCSC{Tv,Ti}

    """
    Intermediate structure holding data of extension
    """
    extmatrix::SparseMatrixExtension{Tv,Ti}
end

"""
$(TYPEDSIGNATURES)

Create empty ExtendablSparseMatrix.

This is a pendant to spzeros.
"""
ExtendableSparseMatrix{Tv,Ti}(m::Integer, n::Integer) where {Tv,Ti<:Integer}=ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),SparseMatrixExtension{Tv, Ti}(m,n))

"""
$(TYPEDSIGNATURES)

Return index corresponding to entry (i,j) in the array of nonzeros,
if the entry exists, otherwise, return 0.
"""
function findindex(S::SparseMatrixCSC{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= S.m && 1 <= j <= S.n); throw(BoundsError()); end
    r1 = Int(S.colptr[j])
    r2 = Int(S.colptr[j+1]-1)
    if r1>r2
        return zero(T)
    end

    # See sparsematrix.jl
    r1 = searchsortedfirst(S.rowval, i, r1, r2, Base.Forward)
    if (r1>r2 ||S.rowval[r1] != i)
        return zero(T)
    end
    return r1
end

"""
$(SIGNATURES)

Find index in CSC matrix and set value if it exists. Otherwise,
set index in extension.
"""
function Base.setindex!(M::ExtendableSparseMatrix, v, i::Integer, j::Integer)
    k=findindex(M.cscmatrix,i,j)
    if k>0
        M.cscmatrix.nzval[k]=v
    else
        M.extmatrix[i,j]=v
    end
end

"""
$(SIGNATURES)

Find index in CSC matrix and return value, if it exists.
Otherwise, return value from extension.
"""
function Base.getindex(M::ExtendableSparseMatrix,i::Integer, j::Integer)
    k=findindex(M.cscmatrix,i,j)
    if k>0
        return M.cscmatrix.nzval[k]
    else
        return M.extmatrix[i,j]
    end
end

"""
$(SIGNATURES)

Size of ExtendableSparseMatrix.
"""
Base.size(E::ExtendableSparseMatrix) = (E.cscmatrix.m, E.cscmatrix.n)

"""
$(SIGNATURES)

Number of nonzeros of ExtendableSparseMatrix.
"""
SparseArrays.nnz(E::ExtendableSparseMatrix)=(nnz(E.cscmatrix)+nnz(E.extmatrix))

# Struct holding pair of value and row
# number, for sorting
mutable struct ColEntry{Tv,Ti<:Integer}
    i::Ti
    v::Tv
end

# Comparison method for sorting
Base.isless(x::ColEntry{Tv, Ti},y::ColEntry{Tv, Ti}) where {Tv,Ti<:Integer} = (x.i<y.i)

function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    # Create new CSC matrix with sorted entries from  CSC matrix S and matrix extension E.
    #
    # This method assumes that there are no entries with the same
    # indices in E and S, therefore  it appears too dangerous for general use and
    # so we don't export it. Generalizations appear to be possible, though.

    @assert(S.m==E.m)
    @assert(S.n==E.n)

    xnnz=nnz(S)+nnz(E)
    colptr=Vector{Ti}(undef,S.n+1)
    rowval=Vector{Ti}(undef,xnnz)
    nzval=Vector{Tv}(undef,xnnz)

    # Detect the maximum column length of E
    E_maxcol=0
    for j=1:S.n
        lcol=0
        k=j
        while k>0
            lcol+=1
            k=E.colptr[k]
        end
        E_maxcol=max(lcol,E_maxcol)
    end

    # pre-allocate column
    col=[ColEntry{Tv,Ti}(0,0) for i=1:E_maxcol]

    inz=1 # counts the nonzero entries in the new matrix

    # loop over all columns
    for j=1:S.n
        # put extension entries into col and sort them
        k=j
        lxcol=0
        while k>0
            if E.rowval[k]>0
                lxcol+=1
                col[lxcol].i=E.rowval[k]
                col[lxcol].v=E.nzval[k]
            end
            k=E.colptr[k]
        end
        sort!(col,1,lxcol, Base.QuickSort, Base.Forward)

        # jointly sort E and S entries  into new matrix data
        colptr[j]=inz
        jcol=1 # counts the entries in col
        k=S.colptr[j]  # counts entries in S
        while true

            # Check if there are entries in S preceding col[jcol]
            if ( nnz(S)>0 &&  (k<S.colptr[j+1]) ) &&
                 (
                     (jcol<=lxcol && S.rowval[k]<col[jcol].i) ||
                     (jcol>lxcol)
                 )
                rowval[inz]=S.rowval[k]
                nzval[inz]=S.nzval[k]
                k+=1
                inz+=1

            # Otherwise, insert next entry from col
            elseif jcol<=lxcol
                rowval[inz]=col[jcol].i
                nzval[inz]=col[jcol].v
                jcol+=1
                inz+=1
            else
                break
            end
        end
    end
    colptr[S.n+1]=inz
    S1=SparseMatrixCSC{Tv,Ti}(S.m,S.n,colptr,rowval,nzval)
    return S1
end

"""
$(TYPEDSIGNATURES)

If there are new entries in extension, create new CSC matrix
and reset extension.
"""
function flush!(M::ExtendableSparseMatrix{Tv,Ti}) where {Tv, Ti<:Integer}
    if nnz(M.extmatrix)>0
        M.cscmatrix=_splice(M.extmatrix,M.cscmatrix)
        M.extmatrix=SparseMatrixExtension{Tv,Ti}(M.cscmatrix.m, M.cscmatrix.n)
    end
    return M
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function SparseArrays.nonzeros(E::ExtendableSparseMatrix)
    @inbounds flush!(E)
    return nonzeros(E.cscmatrix)
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function SparseArrays.rowvals(E::ExtendableSparseMatrix)
    @inbounds flush!(E)
    return rowvals(E.cscmatrix)
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function xcolptrs(E::ExtendableSparseMatrix)
    @inbounds flush!(E)
    return E.cscmatrix.colptr
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function colptrs(E::ExtendableSparseMatrix)
    flush!(E)
    return E.cscmatrix.colptr
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function SparseArrays.findnz(E::ExtendableSparseMatrix)
    flush!(E)
    return findnz(E.cscmatrix)
end

"""
$(TYPEDSIGNATURES)

Drop in replacement for LU factorization.

"""
function LinearAlgebra.lu(E::ExtendableSparseMatrix)
    @inbounds flush!(E)
    return LinearAlgebra.lu(E.cscmatrix)
end

"""
$(SIGNATURES)

Flush and delegate to cscmatrix.
"""
function LinearAlgebra.mul!(r::AbstractArray{Tv,1},
                            E::ExtendableSparse.ExtendableSparseMatrix{Tv,Ti},
                            x::AbstractArray{Tv,1}) where{Tv,Ti<:Integer}
    @inbounds flush!(E)
    return LinearAlgebra.mul!(r,E.cscmatrix,x)
end


text:
using Random
using Plots
import Distributions
include("functions.jl")

mutable struct Net
    weights::Vector{Matrix}
    biases::Vector{Vector}
    f::Vector{Function}

    function Net(sizes::Int...)
        dist = Distributions.Uniform(-0.1, 0.1)

        weights = Vector{Matrix}()
        biases = Vector{Vector}()

        for i in 1:length(sizes)-1
            w = rand(dist, sizes[i+1], sizes[i])
            b = rand(dist, sizes[i+1])

            push!(weights, w)
            push!(biases, b)
        end

        new(weights, biases, [fill(vtanh, length(sizes) - 2); softmax])
    end
end

mean(x) = sum(x)/length(x)
relu(x::Vector) = max.(x, 0)

function softmax(x::Vector)
    m = maximum(x)
    s = exp.(x .- m)
    s/sum(s)
end

function chop(v::Vector, s::Int)
    res = Vector{Vector}()
    l = length(v)

    for i in 1:s:l
        stop = i+s-1
        stop = stop <= l ? stop : l

        push!(res, v[i:stop])
    end

    res
end

function forward(net::Net, x::AbstractVector)
    outputs = Vector{AbstractVector}()
    push!(outputs, x)

    for (W, b, f) in zip(net.weights, net.biases, net.f)
        raw = W*outputs[end] + b
        push!(outputs, f(raw))
    end

    outputs
end

predict(net::Net, x::AbstractVector) = argmax(forward(net, x)[end])

function accuracy(net, x::Vector, t::Vector)
    predictions = [predict(net, i) for i in x]
    cmp = isequal.(predictions, t)
    sum(cmp)/length(x)
end

function update_weights(net::Net, batch_x::Vector, batch_t::Vector, α::Number, regularization::Number)
    out_grad(o, t) = o - [Int(==(i, t)) for i in 1:length(o)]
    hid_grad(W, g, h) = dtanh(h) * (W' * g)
    dW(g, h) = g*h'

    outputs = map(x -> forward(net, x), batch_x)
    layers = collect(zip(outputs...))
    output_layer_gradient = out_grad.(layers[end], batch_t)

    net.biases[end] -= α*mean(output_layer_gradient)
    net.weights[end] -= α*(mean(dW.(output_layer_gradient, layers[end-1])) + regularization * net.weights[end])

    grad = output_layer_gradient
    for i in (length(net.weights)-1):-1:1
        grad = [hid_grad(net.weights[i+1], g, l) for (g, l) in zip(grad, layers[i+1])]

        net.biases[i] -= α*mean(grad)
        net.weights[i] -= α*(mean(dW.(grad, layers[i])) + regularization*net.weights[i])
    end
end

function cost(net, x, t)
    o = forward(net, x)[end]
    -log(o[t])
end

function train(net::Net, x::Vector, t::Vector, epochs::Int, batch_size::Int, rate::Number, decay::Number = 1; regularization::Number = 0)
    plotly()

    indices = 1:length(x)
    test_indices = shuffle(indices)[1:2*batch_size]
    train_indices = filter(i -> !(i in test_indices), indices)

    mc = Vector{Float64}()
    test_accs = Vector{Float64}()
    train_accs = Vector{Float64}()

    for epoch in 1:epochs
        batch_indices = chop(shuffle(train_indices), batch_size)

        costs = nothing
        for bi in batch_indices
            update_weights(net, x[bi], t[bi], rate, regularization)
            costs = [cost(net, x[bi][j], t[bi][j]) for j in 1:length(bi)]
        end

        mean_cost = mean(costs)
        test_acc = accuracy(net, x[test_indices], t[test_indices])
        train_acc = accuracy(net, x[train_indices], t[train_indices])

        println("Epoch $epoch/$epochs")
        println("\tTrain accuracy: $(train_acc * 100) %")
        println("\tTest accuracy: $(test_acc * 100) %")
        println("\tMean cost: $mean_cost")
        println("\tLearning rate: $rate")

        push!(mc, mean_cost)
        push!(test_accs, test_acc)
        push!(train_accs, train_acc)

        rate /= decay
    end

    plt = plot(1:epochs, mc, xlabel = "Epoch", ylabel = "Value", label = "Mean cost", line = [:scatter], grid = false) # dpi = 600, size = (1440, 810)
    plot!(plt, 1:epochs, test_accs, line = [:scatter], label = "Test accuracy")
    plot!(plt, 1:epochs, train_accs, line = [:scatter], label = "Train accuracy")
    display(plt)
end

using Serialization
using Mmap
data = deserialize("trainset/data2")
targets = deserialize("trainset/targets2")

#data = [Float64.(Gray.(d)) for d in data]

N = 310
l = 900
#data = mmap("trainset/trainst3", BitMatrix, (l, N))
#targets = mmap("trainset/targets3", BitMatrix, (3, N))

#data = [data[:, i] for i in 1:N]
#targets = [argmax(targets[:, i]) for i in 1:N]

net = Net(l, 10, 3)
train(net, data, targets, 750, 50, 2e-2, regularization = 1e-2)

ext: jl

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition


text:
@testset "recv_transmit.jl" begin
    @testset begin receiver = TCPReceiver(4000)
        segment_arrives(receiver; with_syn=true, seqno=0, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(1)
        segment_arrives(receiver; data="abcd", seqno=1, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(5)
        @test read!(stream_out(receiver)) == "abcd"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 4
    end

    @testset begin receiver = TCPReceiver(4000)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        segment_arrives(receiver; data="abcd", seqno=isn+1, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn+5)
        @test read!(stream_out(receiver)) == "abcd"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 4

        segment_arrives(receiver; data="efgh", seqno=isn+5, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn+9)
        @test read!(stream_out(receiver)) == "efgh"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 8
    end

    @testset begin receiver = TCPReceiver(4000)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        segment_arrives(receiver; data="abcd", seqno=isn+1, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn+5)

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 4

        segment_arrives(receiver; data="efgh", seqno=isn+5, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn+9)
        @test read!(stream_out(receiver)) == "abcdefgh"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 8
    end

    @testset begin receiver = TCPReceiver(10000)
        max_block_size = 10
        n_rounds = 100;
        bytes_sent = UInt32(0);
        isn = rand(UInt32(0):typemax(UInt32))

        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)
        all_data = ""
        for i in 1:n_rounds
            data = ""
            block_size = rand(1:max_block_size)
            for j in 0:block_size-1
                c = 'a' + ((i + j) % 26)
                data *= c
            end
            all_data *= data
            @test ackno(receiver) == WrappingInt32(isn+bytes_sent+1)
            @test assembled_bytes(receiver) == bytes_sent

            segment_arrives(receiver; seqno=isn+bytes_sent+1, data=data, result=SegmentArrives_Result_OK)
            bytes_sent += block_size
        end
        @test read!(stream_out(receiver)) == all_data

    end
end

main()

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using DataFrames, DataFramesMeta, CSV
using PGFPlots
using Statistics

""" Functions for Fig.3 """

# ---- Average crossvalidated LL across multiple runs
# For results of Fig. 3
function getCrossValResults_multipleruns_incolumnsandrows(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/crossvalidation/"),
        filestr = "",
        # savepath = "./"
        )
    savepath = makehomesavepath("getCrossValResults")
    mkdir(savepath)

    prefix = split(filepath, "/")[end-1]
    runs = getsubfolders(filepath)
    dfresults = DataFrame()
    dfresults[!, :runs] = 1:length(runs)
    for irun in 1:length(runs)
        algorithms = getsubfolders(joinpath(filepath, runs[irun]))
        for i in algorithms
            if irun == 1 # Initialize
                 dfresults[!, Symbol(i)] .= 0.0
            end
            dfalgo = CSV.read(joinpath(filepath, runs[irun], i, "cross_validation/out_of_sample_LL_SUM.csv"), header=false)
            dfresults[irun, Symbol(i)] = dfalgo[1, :Column1]
        end
    end
    CSV.write(joinpath(savepath, prefix*"_resultscrossval_multipleruns"*filestr*".csv"), dfresults,
            delim = " ");

    dfresultsmean = DataFrame()
    algorithms = getsubfolders(joinpath(filepath, runs[1]))
    for i in algorithms
        dfresultsmean[!, Symbol(i*"_mean")] = [mean(dfresults[!, Symbol(i)])]
        dfresultsmean[!, Symbol(i*"_std")] = [std(dfresults[!, Symbol(i)])]
        dfresultsmean[!, Symbol(i*"_stderror")] = [std(dfresults[!, Symbol(i)]) / sqrt(length(dfresults[!, Symbol(i)]))]
    end
    CSV.write(joinpath(savepath, prefix*"_resultscrossval_multipleruns_mean"*filestr*".csv"), dfresultsmean,
            delim = " ");

    dfrows = describe(dfresults, :mean, :std)
    CSV.write(joinpath(savepath, prefix*"_resultscrossval_multipleruns_mean"*filestr*"_rows.csv"), dfrows,
                   delim = " ");

    dfresults, dfresultsmean
end

# Transform the Klaas results that are saved from Matlab in the form we need for
# plotting them as barplots
function getKlassresults_selection(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/KlaasAnalysis/"),
        )

    savepath = makehomesavepath("getKlassresults")
    mkdir(savepath)

    dfalgolist = CSV.read(joinpath(filepath, "algolist.csv"), header=false)

    df_selection_PG = CSV.read(joinpath(filepath, "selection_PG.csv"), header=false)
    df_pxp_PG = CSV.read(joinpath(filepath, "pxp_PG.csv"), header=false)
    df_exp_r_PG = CSV.read(joinpath(filepath, "exp_r_PG.csv"), header=false)
    df_alpha_PG = CSV.read(joinpath(filepath, "alpha_PG.csv"), header=false)

    dfresults_pxp = DataFrame()
    dfresults_exp_r = DataFrame()
    dfresults_alpha = DataFrame()
    nAlgo = size(df_selection_PG,2)
    for i in 1:nAlgo
        dfresults_pxp[!, Symbol(dfalgolist[df_selection_PG[1,i],1])] = [df_pxp_PG[1,i]]
        dfresults_exp_r[!, Symbol(dfalgolist[df_selection_PG[1,i],1])] = [df_exp_r_PG[1,i]]
        dfresults_alpha[!, Symbol(dfalgolist[df_selection_PG[1,i],1])] = [df_alpha_PG[1,i]]
    end

    CSV.write(joinpath(savepath, "resultsKlaas_pxp_PG.csv"), dfresults_pxp, delim = " ");
    CSV.write(joinpath(savepath, "resultsKlaas_exp_r_PG.csv"), dfresults_exp_r, delim = " ");
    CSV.write(joinpath(savepath, "resultsKlaas_alpha_PG.csv"), dfresults_alpha, delim = " ");
end

# ---- Average LLpop across multiple runs and calculate BICpop
# For results of Fig. A.4
function getLLpop_multiplefits_incolumnsandrows(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/fits_multipleruns/"),
        filestr = "",
        # savepath = "./"
        )

    savepath = makehomesavepath("getLLruns")
    mkdir(savepath)

    nr_actions = 3666
    nr_subjects = 21

    prefix = split(filepath, "/")[end-1]
    runs = getsubfolders(filepath)
    dfresultsLL = DataFrame()
    dfresultsLL[!, :runs] = 1:length(runs)
    dfresultsBICactions = DataFrame()
    dfresultsBICactions[!, :runs] = 1:length(runs)
    dfresultsBICsubj = DataFrame()
    dfresultsBICsubj[!, :runs] = 1:length(runs)
    for irun in 1:length(runs)
        algorithms = getsubfolders(joinpath(filepath, runs[irun]))
        for i in algorithms
            if irun == 1 # Initialize
                 dfresultsLL[!, Symbol(i)] .= 0.0
                 dfresultsBICactions[!, Symbol(i)] .= 0.0
                 dfresultsBICsubj[!, Symbol(i)] .= 0.0
            end
            dfalgo = CSV.read(joinpath(filepath, runs[irun], i, "maxLL_fit/LogLikelihood_MAP.csv"), header=false)
            dfalgoParams = CSV.read(joinpath(filepath, runs[irun], i, "maxLL_fit/MAP_params.csv"), header=false)
            nr_params = size(dfalgoParams,1)

            dfresultsLL[irun, Symbol(i)] = dfalgo[1, :Column1]
            dfresultsBICactions[irun, Symbol(i)] = RL_Fit.get_BIClogevidence(dfalgo[1, :Column1], nr_params, nr_actions)
            dfresultsBICsubj[irun, Symbol(i)] = RL_Fit.get_BIClogevidence(dfalgo[1, :Column1], nr_params, nr_subjects)
        end
    end
    CSV.write(joinpath(savepath, prefix*"_fitLL_multipleruns"*filestr*".csv"), dfresultsLL, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICactions_multipleruns"*filestr*".csv"), dfresultsBICactions, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICsubj_multipleruns"*filestr*".csv"), dfresultsBICsubj, delim = " ");

    dfresultsmeanLL = DataFrame()
    dfresultsmeanBICactions = DataFrame()
    dfresultsmeanBICsubjects = DataFrame()
    algorithms = getsubfolders(joinpath(filepath, runs[1]))
    for i in algorithms
        dfresultsmeanLL[!, Symbol(i*"_mean")] = [mean(dfresultsLL[!, Symbol(i)])]
        dfresultsmeanLL[!, Symbol(i*"_std")] = [std(dfresultsLL[!, Symbol(i)])]
        dfresultsmeanLL[!, Symbol(i*"_stderror")] = [std(dfresultsLL[!, Symbol(i)]) / sqrt(length(dfresultsLL[!, Symbol(i)]))]
        #
        dfresultsmeanBICactions[!, Symbol(i*"_mean")] = [mean(dfresultsBICactions[!, Symbol(i)])]
        dfresultsmeanBICactions[!, Symbol(i*"_std")] = [std(dfresultsBICactions[!, Symbol(i)])]
        dfresultsmeanBICactions[!, Symbol(i*"_stderror")] = [std(dfresultsBICactions[!, Symbol(i)]) / sqrt(length(dfresultsBICactions[!, Symbol(i)]))]
        #
        dfresultsmeanBICsubjects[!, Symbol(i*"_mean")] = [mean(dfresultsBICsubj[!, Symbol(i)])]
        dfresultsmeanBICsubjects[!, Symbol(i*"_std")] = [std(dfresultsBICsubj[!, Symbol(i)])]
        dfresultsmeanBICsubjects[!, Symbol(i*"_stderror")] = [std(dfresultsBICsubj[!, Symbol(i)]) / sqrt(length(dfresultsBICsubj[!, Symbol(i)]))]
    end

    CSV.write(joinpath(savepath, prefix*"_fitLL_multipleruns_mean"*filestr*".csv"), dfresultsmeanLL, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICactions_multipleruns_mean"*filestr*".csv"), dfresultsmeanBICactions, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICsubjects_multipleruns_mean"*filestr*".csv"), dfresultsmeanBICsubjects, delim = " ");

    dfrowsLL = describe(dfresultsLL, :mean, :std)
    dfrowsBICactions = describe(dfresultsBICactions, :mean, :std)
    dfrowsBICsubj = describe(dfresultsBICsubj, :mean, :std)
    CSV.write(joinpath(savepath, prefix*"_fitLL_multipleruns_mean"*filestr*"_rows.csv"), dfrowsLL, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICactions_multipleruns_mean"*filestr*"_rows.csv"), dfrowsBICactions, delim = " ");
    CSV.write(joinpath(savepath, prefix*"_fitBICsubjects_multipleruns_mean"*filestr*"_rows.csv"), dfrowsBICsubj, delim = " ");

    dfresultsLL, dfresultsBICactions, dfresultsBICsubj,
    dfresultsmeanLL, dfresultsmeanBICactions, dfresultsmeanBICsubjects
end


text:
module IRTest

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) IRTest

using Test

function run_ir_test_script(script::AbstractString, include, dir::AbstractString)
    if (
        Base.JLOptions().can_inline == 0 ||  # --inline=no
        Base.JLOptions().check_bounds == 1 ||  # --check-bounds=no
        Base.JLOptions().code_coverage != 0  # not --code-coverage=none
    )
        script = isabspath(script) ? script : joinpath(dir, script)
        code = """
        $(Base.load_path_setup_code())
        include($(repr(script)))
        """
        cmd = Base.julia_cmd()
        cmd = `$cmd --check-bounds=no --code-coverage=none --inline=yes`
        if Base.JLOptions().color == 1
            cmd = `$cmd --color=yes`
        end
        @info "Running IR test in a subprocess..." cmd script
        @test success(pipeline(`$cmd -e $code`; stdout=stdout, stderr=stderr))
        @info "Running IR test in a subprocess...DONE"
    else
        include(script)
    end
end

"""
    IRTest.@include(script)

Include a test `script` or run it in an external process if one of the
following flags is specified for the current process:

```sh
--inline=no
--check-bounds=yes
--code-coverage=user
--code-coverage=all
```

Test `script` should contain `@test` to appropriately fail when there
is a failing test.
"""
macro include(script)
    dir = dirname(string(__source__.file))
    esc(:($run_ir_test_script($script, include, $dir)))
end

end # module