text:
#DECORRELATION SEARCH PATTERNS
function permute_source(m::ICA_PWM_Model, models::Vector{Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat, source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},Bool}}; iterates::Integer=length(m.sources)*2, weight_shift_freq::AbstractFloat=PWM_SHIFT_FREQ, weight_shift_dist::Distributions.ContinuousUnivariateDistribution=PWM_SHIFT_DIST, length_change_freq::AbstractFloat=PWM_LENGTHPERM_FREQ, length_perm_range::UnitRange{<:Integer}=LENGTHPERM_RANGE, remote=false)
#weight_shift_dist is given in decimal probability values- converted to log space in permute_source_lengths!
    new_log_Li=-Inf;  iterate = 1
    O = size(obs_array,2);S = length(m.sources)
    new_sources=deepcopy(m.sources);

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources);
        s = rand(1:S)
        clean=Vector{Bool}(trues(O))
        clean[m.mix_matrix[:,s]].=false #all obs with source are dirty

        weight_shift_freq > 0 && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits, length_perm_range))

        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate", "PS from $(m.name)", new_sources, m.source_length_limits, m.mix_matrix, new_log_Li)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate", "PS from $(m.name)", new_sources, m.source_length_limits, m.mix_matrix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function permute_mix(m::ICA_PWM_Model, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; iterates::Integer=10, mix_move_range::UnitRange=1:length(m.mix_matrix), remote=false)
    new_log_Li=-Inf;  iterate = 1; O = size(obs_array,2);
    new_mix=falses(size(m.mix_matrix))

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)
    dirty=false

    while (new_log_Li <= contour || !dirty) && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        mix_moves=rand(mix_move_range)
        mix_moves > length(m.mix_matrix) && (mix_moves = length(m.mix_matrix))

        new_mix, clean = mix_matrix_decorrelate(m.mix_matrix, mix_moves) #generate a decorrelated candidate mix
        c_log_li, c_cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #calculate the model with the candidate mix
        positive_indices=c_cache.>(cache) #obtain any obs indices that have greater probability than we started with
        clean=Vector{Bool}(trues(O)-positive_indices)
        if any(positive_indices) #if there are any such indices
            new_log_Li, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #calculate the new model
            dirty=true
        end
        iterate += 1
    end

    return ICA_PWM_Model("candidate","PM from $(m.name)", m.sources, m.source_length_limits, new_mix, new_log_Li, Vector{Function}()) #no consolidate check is necessary as sources havent changed
end

function perm_src_fit_mix(m::ICA_PWM_Model,  models::Vector{Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat},contour::AbstractFloat,  source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},Bool}}; iterates::Integer=length(m.sources)*2, weight_shift_freq::AbstractFloat=PWM_SHIFT_FREQ, length_change_freq::AbstractFloat=PWM_LENGTHPERM_FREQ, length_perm_range::UnitRange{<:Integer}=LENGTHPERM_RANGE,weight_shift_dist::Distributions.ContinuousUnivariateDistribution=PWM_SHIFT_DIST, remote=false)
    new_log_Li=-Inf;  iterate = 1
    O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    fit_mix in m.permute_blacklist ? (new_bl=m.permute_blacklist) : (new_bl=Vector{Function}())

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix); tm_one=deepcopy(m.mix_matrix);
        clean=Vector{Bool}(trues(O))
        s = rand(1:S);
        clean[new_mix[:,s]].=false #all obs starting with source are dirty

        new_mix[:,s].=false;tm_one[:,s].=true

        weight_shift_freq > 0 && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits, length_perm_range))

        l,zero_cache=IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true)
        l,one_cache=IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, tm_one, true, true)
        fit_mix=one_cache.>=zero_cache

        if REVCOMP #if we're looking at reverse strands, we want to see if sources fit better on the reverse strand, and if so switch over to the revcomp source and use its fitted mix
            revsrc=(revcomp_pwm(new_sources[s][1]),new_sources[s][2])
            revsrcs=deepcopy(new_sources)
            revsrcs[s]=revsrc
            l, revsrc_cache=IPM_likelihood(revsrcs, obs_array, obs_lengths, bg_scores, tm_one, true, true)
            rfit_mix=revsrc_cache.>=zero_cache
            lps(revsrc_cache[rfit_mix])>lps(one_cache[fit_mix]) && (new_sources=revsrcs;fit_mix=rfit_mix) #if the total likelihood contribution for the revsource fit is greater than the source fit, take the revsource
        end

        new_mix[:,s]=fit_mix

        clean[fit_mix].=false #all obs ending with source are dirty

        new_log_Li = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, false, zero_cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","PSFM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li, new_bl)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","PSFM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li, new_bl), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function fit_mix(m::ICA_PWM_Model, models::Vector{Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}; remote=false)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix); test_mix=falses(size(m.mix_matrix))
    new_bl=[fit_mix]

    l,zero_cache=IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, test_mix, true, true)
    for (s,source) in enumerate(m.sources)
        test_mix=falses(size(m.mix_matrix))
        test_mix[:,s].=true
        l,src_cache=IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, test_mix, true, true)
        fit_mix=src_cache.>=zero_cache

        if REVCOMP #if we're looking at reverse strands, we want to see if sources fit better on the reverse strand, and if so switch over to the revcomp source and use its fitted mix
            revsrc=(revcomp_pwm(source[1]),source[2])
            revsrcs=deepcopy(new_sources)
            revsrcs[s]=revsrc
            l, revsrc_cache=IPM_likelihood(revsrcs, obs_array, obs_lengths, bg_scores, test_mix, true, true)
            rfit_mix=revsrc_cache.>=zero_cache
            lps(revsrc_cache[rfit_mix])>lps(src_cache[fit_mix]) && (new_sources=revsrcs;fit_mix=rfit_mix) #if the total likelihood contribution for the revsource fit is greater than the source fit, take the revsource
        end

        new_mix[:,s]=fit_mix
    end
    new_mix==m.mix_matrix ? (new_log_Li=-Inf) : (new_log_Li = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, new_mix, true)) #bail if the mix matrix doesnt change
    new_log_Li in [model.log_Li for model in models] && (new_log_Li=-Inf) #an existing identical log_Li is almost certainly a fitted model from another origin with the same sources, in this case abort to prevent dupes

    return ICA_PWM_Model("candidate","FM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li, new_bl) #no consolidate check necessary as no change to sources
end

function random_decorrelate(m::ICA_PWM_Model, models::Vector{Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat, source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},Bool}}; iterates::Integer=length(m.sources)*2, source_permute_freq::AbstractFloat=SRC_PERM_FREQ, weight_shift_freq::AbstractFloat=PWM_SHIFT_FREQ, length_change_freq::AbstractFloat=PWM_LENGTHPERM_FREQ, weight_shift_dist::Distributions.ContinuousUnivariateDistribution=PWM_SHIFT_DIST, mix_move_range::UnitRange=1:size(m.mix_matrix,1), remote=false)
    new_log_Li=-Inf;  iterate = 1
    O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)

    a, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))
        s = rand(1:S)
        clean[new_mix[:,s]].=false #all obs starting with source are dirty
        rand() < source_permute_freq && (new_sources[s]=permute_source_weights(new_sources[s], weight_shift_freq, weight_shift_dist))
        rand() < length_change_freq && (new_sources[s]=permute_source_length(new_sources[s], source_priors[s], m.source_length_limits))
        new_mix[:,s]=mixvec_decorrelate(new_mix[:,s],rand(mix_move_range))
        clean[new_mix[:,s]].=false #all obs ending with source are dirty
        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true, cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate", "RD from $(m.name)", new_sources, m.source_length_limits, new_mix, new_log_Li)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","RD from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function shuffle_sources(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; remote=false)
    new_log_Li=-Inf; O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    remote ? (merger_m = deserialize(rand(models).path)) : (merger_m = remotecall_fetch(deserialize, 1, rand(models).path))#randomly select a model to merge

    svec=[1:S...]

    while new_log_Li <= contour && length(svec)>0 #until we produce a model more likely than the lh contour or no more sources to attempt merger
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        s = popat!(svec,rand(1:length(svec))) #randomly select a source to merge

        new_sources[s]=merger_m.sources[s]; new_mix[:,s].=merger_m.mix_matrix[:,s] #shuffle in the merger model source from this index

        clean[m.mix_matrix[:,s]].=false #mark dirty any obs that have the source
        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","SS from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li,Vector{Function}())) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","SS from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function accumulate_mix(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; remote=false)
    new_log_Li=-Inf; O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    new_bl=Vector{Function}()

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    remote ? (merger_m = deserialize(rand(models).path)) : (merger_m = remotecall_fetch(deserialize, 1, rand(models).path))#randomly select a model to merge

    svec=[1:S...]

    while new_log_Li <= contour && length(svec)>0 #until we produce a model more likely than the lh contour or no more sources to attempt merger
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        s = popat!(svec,rand(1:length(svec))) #randomly select a source to merge
        distvec=[pwm_distance(src[1],m_src[1]) for src in m.sources, m_src in merger_m.sources]
        S > 1 ? merge_s=findmin(distvec)[2][2] : merge_s=findmin(distvec)[2]

        new_mix[:,s].=[new_mix[o,s] || merger_m.mix_matrix[o,merge_s] for o in 1:size(new_mix,1)] #accumulate the mix vector for this source

        clean[m.mix_matrix[:,s]].=false #mark dirty any obs that have the source
        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood

        if new_log_Li <= contour #if accumulating on the host source doesnt work, try copying over the merger source
            new_sources[s]=merger_m.sources[merge_s]
            new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
        end
    end

    (new_log_Li == m.log_Li || new_log_Li in [model.log_Li for model in models]) && (new_log_Li=-Inf) #accumulate can sometimes duplicate models if source mix is identical

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","AM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","AM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function distance_merge(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; remote=false)
    new_log_Li=-Inf; O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)

    a, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true)

    remote ? (merger_m = deserialize(rand(models).path)) : (merger_m = remotecall_fetch(deserialize, 1, rand(models).path)) #randomly select a model to merge

    svec=[1:S...]

    while new_log_Li <= contour && length(svec)>0 #until we produce a model more likely than the lh contour or no more sources to attempt merger
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        s = popat!(svec,rand(1:length(svec))) #randomly select a source to merge
        merge_s=most_dissimilar(new_mix[:,s],merger_m.mix_matrix) #find the source in the merger model whose mixvec is most dissimilar to the one selected

        clean[new_mix[:,s]].=false #mark dirty any obs that start with the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source
        new_mix[:,s] = merger_m.mix_matrix[:,merge_s] #copy the mixvector (without which the source will likely be highly improbable)
        clean[new_mix[:,s]].=false #mark dirty any obs that end with the source

        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #assess likelihood
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","DM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li,Vector{Function}())) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","DM from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function similarity_merge(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; remote=false)
    new_log_Li=-Inf; O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources)

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    remote ? (merger_m = deserialize(rand(models).path)) : (merger_m = remotecall_fetch(deserialize, 1, rand(models).path))#randomly select a model to merge

    svec=[1:S...]

    while new_log_Li <= contour && length(svec)>0 #until we produce a model more likely than the lh contour or no more sources to attempt merger
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        s = popat!(svec,rand(1:length(svec))) #randomly select a source to merge
        merge_s=most_similar(m.mix_matrix[:,s],merger_m.mix_matrix) #obtain the source in the merger model whose mixvec is most similar to the one in the original

        clean[m.mix_matrix[:,s]].=false #mark dirty any obs that have the source
        new_sources[s] = merger_m.sources[merge_s] #copy the source, but dont copy the mixvector
        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","SM from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li,Vector{Function}())) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","SM from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function reinit_src(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat},contour::AbstractFloat, source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},Bool}}; iterates::Integer=length(m.sources)*2, remote=false)
    new_log_Li=-Inf;  iterate = 1
    O = size(obs_array,2); S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix); tm_one=deepcopy(m.mix_matrix);
        clean=Vector{Bool}(trues(O))
        s = rand(1:S);
        clean[new_mix[:,s]].=false #all obs starting with source are dirty

        new_mix[:,s].=false;tm_one[:,s].=true

        new_sources[s] = init_logPWM_sources([source_priors[s]], m.source_length_limits)[1] #reinitialize the source

        l,zero_cache=IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true)
        l,one_cache=IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, tm_one, true, true)

        fit_mix=one_cache.>=zero_cache

        new_mix[:,s]=fit_mix

        clean[fit_mix].=false #all obs ending with source are dirty

        new_log_Li = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, false, zero_cache, clean)
        iterate += 1
    end

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","RS from $(m.name)", new_sources, m.source_length_limits, new_mix, new_log_Li,Vector{Function}())) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","RS from $(m.name)",new_sources, m.source_length_limits, new_mix, new_log_Li), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function erode_model(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; info_thresh::AbstractFloat=EROSION_INFO_THRESH, remote=false)
    new_log_Li=-Inf; O = size(obs_array,2);
    new_sources=deepcopy(m.sources)

    erosion_sources=Vector{Integer}()
    for (s,src) in enumerate(m.sources)
        pwm,pi=src
        if size(pwm,1)>m.source_length_limits[1]  #do not consider eroding srcs at min length limit
            infovec=get_pwm_info(pwm)
            any(info-><(info, info_thresh),infovec) && push!(erosion_sources,s)
        end
    end

    length(erosion_sources)==0 && return ICA_PWM_Model("candidate","EM from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li,[erode_model])#if we got a model we cant erode bail out with a model marked -Inf lh and with EM blacklisted

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    while new_log_Li <= contour && length(erosion_sources) > 0 #until we produce a model more likely than the lh contour or there are no more sources to erode
        clean=Vector{Bool}(trues(O))
        s=popat!(erosion_sources,rand(1:length(erosion_sources)))

        new_sources[s]=erode_source(new_sources[s], m.source_length_limits, info_thresh)
        clean[m.mix_matrix[:,s]].=false

        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
    end

    new_log_Li <= contour ? (blacklist=[erode_model]) : (blacklist=Vector{Function}())

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","EM from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li, blacklist)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","EM from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li, blacklist), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

function info_fill(m::ICA_PWM_Model, models::AbstractVector{<:Model_Record}, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat; remote=false)
    new_log_Li=-Inf; O = size(obs_array,2); S = length(m.sources); iterate=1
    new_sources=deepcopy(m.sources)

    a, cache = IPM_likelihood(m.sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true)

    svec=[1:S...]

    while new_log_Li <= contour && length(svec)>0 #until we produce a model more likely than the lh contour or no more sources to attempt infofill
        new_sources=deepcopy(m.sources)
        clean=Vector{Bool}(trues(O))
        s = popat!(svec,rand(1:length(svec)))
        pwm=new_sources[s][1]
        fill_idx=findmin(get_pwm_info(pwm))[2]
        fill_bases=[1,2,3,4]
        fill_candidate=findmax(pwm[fill_idx,:])[2]
        deleteat!(fill_bases, findfirst(b->b==fill_candidate,fill_bases))
        new_sources[s][1][fill_idx,:].=-Inf; new_sources[s][1][fill_idx,fill_candidate]=0.
        clean[m.mix_matrix[:,s]].=false
        new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood

        while new_log_Li <=contour && length(fill_bases)>0
            fill_candidate=popat!(fill_bases,rand(1:length(fill_bases)))
            new_sources[s][1][fill_idx,:].=-Inf; new_sources[s][1][fill_idx,fill_candidate]=0.
            new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, m.mix_matrix, true, true, cache, clean) #assess likelihood
        end
    end

    new_log_Li <= contour ? (blacklist=[info_fill]) : (blacklist=Vector{Function}())

    cons_check, cons_idxs = consolidate_check(new_sources)
    cons_check ? (return ICA_PWM_Model("candidate","IF from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li, blacklist)) : (return consolidate_srcs(cons_idxs, ICA_PWM_Model("candidate","IF from $(m.name)",new_sources, m.source_length_limits, m.mix_matrix, new_log_Li, blacklist), obs_array, obs_lengths, bg_scores, contour, models; remote=remote))
end

full_perm_funcvec=[permute_source, permute_mix, perm_src_fit_mix, fit_mix, random_decorrelate, shuffle_sources, accumulate_mix, distance_merge, similarity_merge, reinit_src, erode_model, info_fill]


text:
include("./src/MiniPL.jl")
using .MiniPL
const m = MiniPL

if length(ARGS) != 1
  println("USAGE: julia run.jl sourcefile")
  exit(1)
end

filename = ARGS[1]

if !isfile(filename)
  println("File $(filename) not found.")
  exit(1)
end

source = open(ARGS[1]) do file
  read(file, String)
end

source = try
  ascii(source)
catch e
  println("The source code contains non-ascii characters.")
  exit(1)
end

try
  m.run(source)
catch e
  showerror(stdout, e)
end

println()
exit(0)

ext: jl

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
"""
Given an Array, `a`, and a face, `f`, return a tuple of numbers
interpreting the values and offset of the face as indices into `A`.

Note: This is not bounds checked. It is recommended that you use
`checkbounds` to confirm the indices are safe for loops.
Also be aware when writing generic code that faces may be of more than
3 vertices.
"""
@generated function Base.getindex{T,N,FD,FT,Offset}(a::Array{T,N},
                                               f::Face{FD, FT, Offset})
    v = Expr(:tuple)
    for i = 1:FD
        push!(v.args, Expr(:call, Base.unsafe_getindex, :a, :(Int(f[$i])-Int(Offset))))
    end
    :($(v)::NTuple{FD,T})
end

function setindex!{T,N,FD,FT,Offset}(a::Array{T,N}, b::Array{T,N}, f::Face{FD, FT, Offset})
	for i=1:FD
    	a[onebased(f,i)] = b[i]
    end
end

convert{T1<:Face}(::Type{T1}, f::T1) = f
convert{T1<:Face, T2<:Face}(::Type{T1}, f::T2) = T1(f)

# Silly duplication, but call(::FixedVector, ::Any) = convert is overloaded in FixedSizeArrays
@compat (::Type{F}){F<:Face}(f::F) = f

@compat function (::Type{Face{N, T, O}}){T, T2, O, N}(f::Face{N, T2, O})
    Face{N, T, O}(convert(NTuple{N, T}, getfield(f, 1)))
end
immutable IndexConvertFunc{T1, T2}
	f::T2
end
@compat function (ifunc::IndexConvertFunc{Face{N,T1,O1}, Face{N,T2,O2}}){N,T1,T2,O1,O2}(i)
	Int(ifunc.f[i])+Int(O1)-Int(O2)
end
@compat function (T::Type{Face{N, T1, O1}}){N, T1, O1, F<:Face}(f::F)
	map(IndexConvertFunc{T,F}(f), T)
end

function Face{T<:Number}(vals::T...)
    Face{length(vals), T, 0}(vals...)
end
onebased{N,T,Offset}(face::Face{N,T,Offset}, i) = T(Int(face[i]) - Int(Offset))
zerobased{N,T,Offset}(face::Face{N,T,Offset}, i) = T(Int(face[i]) - Int(Offset) - 1)
offsetbased{N,T,Offset}(face::Face{N,T,Offset}, i, offset) = T(Int(face[i]) - Int(Offset) + Int(offset))

export onebased, zerobased


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

function test_linear_cut_cell_assembly()
    polyorder = 1
    numqp = 2

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)
    nf = CutCell.number_of_basis_functions(basis)
    mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [1, 1], nf)

    normal = [1.0, 0.0]
    x0 = [0.5, 0.0]
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> plane_distance_function(x, normal, x0), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    dx = 0.1
    e11 = dx / 2.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22
    s11 = (lambda + 2mu) * e11 + lambda * e22

    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    CutCell.apply_dirichlet_bc!(matrix, rhs, 9, 0.0)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 10, 0.0)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 11, 0.0)

    CutCell.apply_dirichlet_bc!(matrix, rhs, 5, dx)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 7, dx)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    testdisp = [
        0.0 0.0 dx dx 0.0 0.0 dx dx
        0.0 dy 0.0 dy 0.0 dy 0.0 dy
    ]
    @test allapprox(disp, testdisp, 1e2eps())
end

function test_curved_interface_assembly()
    polyorder = 2
    numqp = 4

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)
    nf = CutCell.number_of_basis_functions(basis)
    mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [1, 1], nf)

    xc = [2.0, 0.5]
    rad = 1.0
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> circle_distance_function(x, xc, rad), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    dx = 0.1
    e11 = dx / 2.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22
    s11 = (lambda + 2mu) * e11 + lambda * e22

    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    CutCell.apply_dirichlet_bc!(matrix, rhs, [10, 11, 12], 1, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [10], 2, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [7, 8, 9], 1, dx, 2)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    testdisp = [
        0.0 0.0 0.0 dx/2 dx/2 dx/2 dx dx dx
        0.0 dy/2 dy 0.0 dy/2 dy 0.0 dy/2 dy
    ]
    testdisp = repeat(testdisp, outer = (1, 2))

    @test allapprox(disp, testdisp, 1e2eps())
end

test_linear_cut_cell_assembly()
test_curved_interface_assembly()


text:
"""
    LangString <: AbstractString

String with a language tag.

    LangString{locale}(content)

A language tag is the locale name (without format), e.g
`:zh_CN`, `:fr_FR`. A full list of locales can be found
here:

https://docs.moodle.org/dev/Table_of_locales
"""
struct LangString{L, T <: AbstractString} <: AbstractString
    content::T

    LangString{L}(content::T) where {L, T} = new{L, T}(content)
    LangString(content::T) where T = new{:en_US, T}(content) # default language is English
end

"""
    language(s) -> String

Returns the name of language in "<locale/local language>" format
"""
function language end

"""
    locale(s) -> String

Return the locale, e.g
`:zh_CN.UTF-8`, `:fr_FR.UTF-8`. A full list of locales can be found
here:

https://docs.moodle.org/dev/Table_of_locales
"""
function locale end

language(s::LangString{:zh_CN}) = "zh_CN/中文"
language(s::LangString{:zh_TW}) = "zh_TW/中文"
language(s::LangString{:fr_FR}) = "fr_FR/français"

locale(s::LangString{L, String}) where L = "$L.UTF-8"

# forward String interface
Base.write(io::IO, s::LangString) = Base.write(io::IO, s.content)
Base.show(io::IO, s::LangString) = Base.show(io, s.content)

Compat.firstindex(s::LangString) = Compat.firstindex(s.content)
Compat.lastindex(s::LangString) = Compat.lastindex(s.content)

Base.start(s::LangString) = start(s.content)
Base.next(s::LangString, i) = next(s.content, i)
Base.done(s::LangString, i) = done(s.content, i)
if isdefined(Base, :iterate)
    Base.iterate(s::LangString, i::Int) = iterate(s.content, i)
end

Base.nextind(s::LangString, i::Int) = nextind(s.content, i)
Base.prevind(s::LangString, i::Int) = prevind(s.content, i)
Base.eachindex(s::LangString) = eachindex(s.content)
Base.length(s::LangString) = length(s.content)

Base.getindex(s::LangString, i::Integer) = getindex(s.content, i)
Base.getindex(s::LangString, i::Int) = getindex(s.content, i) # for method ambig in Julia 0.6
Base.getindex(s::LangString, i::UnitRange{Int}) = getindex(s.content, i)
Base.getindex(s::LangString, i::UnitRange{<:Integer}) = getindex(s.content, i)
Base.getindex(s::LangString, i::AbstractVector{<:Integer}) = getindex(s.content, i)
Compat.codeunit(s::LangString, i::Integer) = codeunit(s.content, i)
Compat.codeunit(s::LangString) = codeunit(s.content)
Compat.ncodeunits(s::LangString) = ncodeunits(s.content)
Compat.codeunits(s::LangString) = codeunits(s.content)
Base.sizeof(s::LangString) = sizeof(s.content)
Base.isvalid(s::LangString, i::Integer) = isvalid(s.content, i)
Base.pointer(s::LangString) = pointer(s.content)
Base.IOBuffer(s::LangString) = IOBuffer(s.content)
Base.unsafe_convert(T::Union{Type{Ptr{UInt8}},Type{Ptr{Int8}},Cstring}, s::LangString) = Base.unsafe_convert(T, s.content)

struct LangMismatchError{T <: Tuple} <: Exception
    got::T
end

LangMismatchError(langs::Symbol...) = LangMismatchError(langs)

Base.show(io::IO, e::LangMismatchError{L}) where L = print(io, "Language mismatch, got: ", e.got)

macro zh_str(s)
    LangString{:zh_CN}(s)
end

macro fr_str(s)
    LangString{:fr_FR}(s)
end

macro lang(locale, s)
    LangString{locale}(s)
end


text:
include("IntSets.jl")
include("iterate_jkp.jl")
include("ldlt.jl")
include("ldlt_api.jl")
include("selinv.jl")
include("selinv_api.jl")

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
#####
##### NetParams Options
#####

# baseline_network = ResNetHP(
#   num_filters=128,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.1)

# more_filters = ResNetHP(
#   num_filters=256,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.1)

# less_filters = ResNetHP(
#   num_filters=64,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.1)

# more_blocks = ResNetHP(
#   num_filters=128,
#   num_blocks=30,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.1)

# less_blocks = ResNetHP(
#   num_filters=128,
#   num_blocks=10,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.1)

# more_head_filters = ResNetHP(
#   num_filters=128,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=128,
#   num_value_head_filters=128,
#   batch_norm_momentum=0.1)

# less_head_filters = ResNetHP(
#   num_filters=128,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=32,
#   num_value_head_filters=32,
#   batch_norm_momentum=0.1)

# more_batch_momentum = ResNetHP(
#   num_filters=128,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.25)

# less_batch_momentum = ResNetHP(
#   num_filters=128,
#   num_blocks=20,
#   conv_kernel_size=(3, 3),
#   num_policy_head_filters=64,
#   num_value_head_filters=64,
#   batch_norm_momentum=0.05)

#####
##### LearningParams Options
#####

# baseline_learning = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.0075,
#     lr_high=0.05,
#     lr_low=0.0025,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# mom_99 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.0075,
#     lr_high=0.05,
#     lr_low=0.0025,
#     momentum_high=0.99,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# mom_97 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.0075,
#     lr_high=0.05,
#     lr_low=0.0025,
#     momentum_high=0.97,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

mom_95 = LearningParams(
  use_gpu=true,
  use_position_averaging=true,
  samples_weighing_policy=LOG_WEIGHT,
  batch_size=128,
  loss_computation_batch_size=256,
  optimiser=CyclicNesterov(
    lr_base=0.0075,
    lr_high=0.05,
    lr_low=0.0025,
    momentum_high=0.95,
    momentum_low=0.8),
  l2_regularization=1e-5,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=1,
  max_batches_per_checkpoint=25000,
  num_checkpoints=1)

# mom_9 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.0075,
#     lr_high=0.05,
#     lr_low=0.0025,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# lr_1_0 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.5,
#     lr_high=1.0,
#     lr_low=0.1,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# lr_0_5 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.25,
#     lr_high=0.5,
#     lr_low=0.05,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# lr_0_25 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.125,
#     lr_high=0.25,
#     lr_low=0.025,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

# lr_0_1 = LearningParams(
#   use_gpu=true,
#   use_position_averaging=true,
#   samples_weighing_policy=LOG_WEIGHT,
#   batch_size=128,
#   loss_computation_batch_size=256,
#   optimiser=CyclicNesterov(
#     lr_base=0.05,
#     lr_high=0.1,
#     lr_low=0.01,
#     momentum_high=0.9,
#     momentum_low=0.8),
#   l2_regularization=1e-5,
#   nonvalidity_penalty=1.,
#   min_checkpoints_per_epoch=1,
#   max_batches_per_checkpoint=25000,
#   num_checkpoints=1)

#####
##### Other Params
#####

self_play = SelfPlayParams(
  num_games=1000,
  reset_mcts_every=1,
  mcts=MctsParams(
    use_gpu=true,
    num_workers=32,
    num_iters_per_turn=1000,
    cpuct=4.0,
    prior_temperature=0.7,
    temperature=PLSchedule([0, 24], [1.0, 0.5]),
    dirichlet_noise_ϵ=0.25,
    dirichlet_noise_α=1.25))

arena = ArenaParams(
  num_games=2,
  reset_mcts_every=1,
  flip_probability=0.1,
  update_threshold=0.5,
  play_swiss_style=true,
  mcts=MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.2),
    dirichlet_noise_ϵ=0.05))

#####
##### Full Params Options
#####

# baseline_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=baseline_learning,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# mom_99_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=mom_97,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# mom_97_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=mom_99,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

mom_95_params = Params(
  arena=arena,
  self_play=self_play,
  learning=mom_95,
  num_iters=100,
  ternary_rewards=false,
  use_symmetries=true,
  memory_analysis=MemAnalysisParams(
    num_game_stages=4),
  mem_buffer_size=PLSchedule(
  [      5,        57],
  [208_000, 1_560_000]))

# mom_9_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=mom_9,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# lr_1_0_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=lr_1_0,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# lr_0_5_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=lr_0_5,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# lr_0_25_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=lr_0_25,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

# lr_0_1_params = Params(
#   arena=arena,
#   self_play=self_play,
#   learning=lr_0_1,
#   num_iters=100,
#   ternary_rewards=false,
#   use_symmetries=true,
#   memory_analysis=MemAnalysisParams(
#     num_game_stages=4),
#   mem_buffer_size=PLSchedule(
#   [      5,        57],
#   [208_000, 1_560_000]))

#####
##### Trials
#####

optimal_network = ResNetHP(
  num_filters=256,
  num_blocks=25,
  conv_kernel_size=(3, 3),
  num_policy_head_filters=128,
  num_value_head_filters=128,
  batch_norm_momentum=0.1)

trials = []
push!(trials, (label="opt_params", params=mom_95_params, netparams=optimal_network))

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
using SphericalGeometry
using Test

@testset "SphericalGeometry.jl" begin
    tests = ["utility", "types", "distances", "angles", "points", "operations", "area"]
    for t in tests
        include("$(t).jl")
    end
end


text:
mutable struct MaximinObjective <: AbstractObjective
    name::String
    sense::String
    cores::Vector{MaximinObjectiveCore}

    MaximinObjective(cores::Vector{MaximinObjectiveCore}) = new("maximin", "max", cores)
    MaximinObjective() = new("maximin", "max", MaximinObjectiveCore[])
end


text:
"""
Structure returned by correlation functions (`l2`, `s2` and `c2`).

This structure holds correlation data computated along specified
directions. To access those data use the `mean` function.
"""
struct CorrelationData
    directions :: Vector{Symbol}
    success    :: Dict{Symbol, Vector{Float64}}
    total      :: Dict{Symbol, Vector{Float64}}
end

function CorrelationData(len        :: Integer,
                         directions :: Vector{Symbol})
    success = Dict(map(x -> x => zeros(len), directions))
    total   = Dict(map(x -> x => zeros(len), directions))
    return CorrelationData(directions, success, total)
end

# Copier
function CorrelationData(data :: CorrelationData)
    directions = data.directions
    success    = data.success
    total      = data.total

    return CorrelationData(directions,
                           Dict(x => success[x] for x in directions),
                           Dict(x => total[x]   for x in directions))
end

function Base.getindex(x :: CorrelationData, i)
    directions = x.directions
    if i ∉ directions
        error("Requested direction $i is not among calculated directions $directions")
    end

    success = x.success[i]
    # Define correlation function f(a, x) = 0 for x > dimension of a
    total = (x -> max(x, 1)).(x.total[i])
    return success ./ total
end

function Base.show(io :: IO, x :: CorrelationData)
    directions = x.directions
    corr = reduce(hcat, x[direction] for direction in directions)
    pretty_table(io, corr; header = directions)
end

function Base.iterate(x :: CorrelationData, state = nothing)
    next = state != nothing ? iterate(x.directions, state) : iterate(x.directions)
    if next != nothing
        direction, next = next
        return (direction, x[direction]), next
    else
        return nothing
    end
end

Base.length(x :: CorrelationData) = let dir = first(x.directions); length(x[dir]) end

import StatsBase: mean
"""
    mean(data :: CorrelationData, directions::Vector{Symbol})
    mean(data :: CorrelationData)

Return mean of correlation function stored in `data` calculated along
directions stored in array `directions`. Elements of `directions` must
be members of `known_directions`. If `directions` contain only one
direction, values of correlation function computed along that
direction are returned. If `mean` called without `directions`
argument, mean of all computed directions is returned.
"""
function mean end

function mean(data :: CorrelationData, directions::Vector{Symbol})
    if !all(x -> x ∈ data.directions, directions)
        error("Correlation function is not computed for directions $directions")
    end

    success = mapreduce(x -> data.success[x],                  +, directions)
    # Once again, define correlation function f(a, x) = 0
    # for x > dimension of a.
    total   = mapreduce(x -> (y -> max(y, 1)).(data.total[x]), +, directions)

    return success ./ total
end

function mean(data :: CorrelationData)
    return mean(data, data.directions)
end

"""
    directions(data :: CorrelationData)

Return directions along which a correlation function is computed.

# Examples
```jldoctest
julia> directions(l2(rand(0:1, (50, 10)), 1))
2-element Vector{Symbol}:
 :x
 :y
```
"""
directions(data :: CorrelationData) = data.directions

"""
    direction ∈ corrdata

Return `true` in correlation data `corrdata` is computed for a
direction `direction`.
"""
Base.in(direction :: Symbol, data :: CorrelationData) = direction ∈ directions(data)


text:
@testset "aspect ratio" begin
    @test spy(rand(2, 2)).graphics.grid |> size == (5, 2)  # 2 + (2x1) = 5 (columns padding)
    @test spy(rand(4, 4)).graphics.grid |> size == (5, 2)
    @test spy(rand(8, 4)).graphics.grid |> size == (5, 2)
    @test spy(rand(12, 4)).graphics.grid |> size == (5, 3)
    @test spy(rand(8, 8)).graphics.grid |> size == (5, 2)
    @test spy(rand(8, 12)).graphics.grid |> size == (6, 2)
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 80, 80, 0.15), fix_ar = true)
    test_ref(
        "references/spy/fix_aspect_ratio_80x80_.txt",
        @show_col(p, :displaysize => T_SZ)
    )
end

@testset "sizing" begin
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 0, 0, 0.1))
    test_ref("references/spy/default_0x0.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 10, 10, 0.15))
    @test p isa Plot
    test_ref("references/spy/default_10x10.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 10, 15, 0.15))
    test_ref("references/spy/default_10x15.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 15, 10, 0.15))
    test_ref("references/spy/default_15x10.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 2000, 200, 0.0001))
    test_ref("references/spy/default_2000x200.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 2000, 0.0001))
    test_ref("references/spy/default_200x2000.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 2000, 2000, 0.1))
    test_ref("references/spy/default_overdrawn.txt", @show_col(p, :displaysize => T_SZ))
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 200, 0.001))
    test_ref(
        "references/spy/default_200x200_normal.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    test_ref(
        "references/spy/default_200x200_normal_nocolor.txt",
        @show_nocol(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = @inferred spy(Matrix(_stable_sprand(RNG, 200, 200, 0.001)))
    test_ref(
        "references/spy/default_200x200_normal.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 200, 0.001), width = 10)
    test_ref(
        "references/spy/default_200x200_normal_small.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 200, 0.001), height = 5)
    test_ref(
        "references/spy/default_200x200_normal_small.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 200, 0.001), height = 5, width = 20)
    test_ref(
        "references/spy/default_200x200_normal_misshaped.txt",
        @show_col(p, :displaysize => T_SZ)
    )
end

@testset "parameters" begin
    seed!(RNG, 1337)
    p = @inferred spy(_stable_sprand(RNG, 200, 200, 0.001), color = :green)
    test_ref(
        "references/spy/parameters_200x200_green.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    test_ref(
        "references/spy/parameters_200x200_green_nocolor.txt",
        @show_nocol(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = spy(
        _stable_sprand(RNG, 200, 200, 0.001),
        title = "Custom Title",
        canvas = DotCanvas,
        border = :ascii,
    )
    test_ref(
        "references/spy/parameters_200x200_dotcanvas.txt",
        @show_col(p, :displaysize => T_SZ)
    )
    seed!(RNG, 1337)
    p = spy(_stable_sprand(RNG, 200, 200, 0.99), show_zeros = true)
    test_ref(
        "references/spy/parameters_200x200_zeros.txt",
        @show_col(p, :displaysize => T_SZ)
    )
end


text:
using DAQnidaqmx
using Test

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


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ a93d84c9-400f-4624-a67d-cec990f3d822
begin
	using PlutoUI, LinearAlgebra
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ afe8da5c-4e4e-4ce7-87de-e3b201a4c133
using Arpack

# ╔═╡ 8db1c1b2-8b5f-402c-8441-768c1a943c38
using SparseArrays

# ╔═╡ f04cb86f-94aa-418d-82f9-d9ef12ca90fd
md"""
# Singular Value Decomposition - Jacobi and Lanczos Methods

Since computing the SVD of $A$ can be seen as computing the EVD of the symmetric matrices $A^*A$, $AA^*$, or $\begin{bmatrix}0 & A \\ A^* & 0 \end{bmatrix}$, simple modifications of the corresponding EVD algorithms yield version for computing the SVD.

For more details on one-sided Jacobi method, see [Z. Drmač, Computing Eigenvalues and Singular Values to High Relative Accuracy, pp. 59.1-59.21](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

__Prerequisites__

The reader should be familiar with concepts of singular values and vectors, related perturbation theory, and algorithms, and Jacobi and Lanczos methods for the symmetric eigenvalue decomposition.

__Competences__

The reader should be able to recognise matrices which warrant high relative accuracy and to apply Jacobi method to them. The reader should be able to recognise matrices to which Lanczos method can be efficiently applied and do so.
"""

# ╔═╡ 2eff5f1f-0989-4292-8f04-80983c1bd133
md"""
# One-sided Jacobi method

Let $A\in\mathbb{R}^{m\times n}$ with $\mathop{\mathrm{rank}}(A)=n$ (therefore, $m\geq n$) and $A=U\Sigma V^T$ its thin SVD.

## Definition

Let $A=BD$, where $D=\mathop{\mathrm{diag}} (\| A_{:,1}\|_2, \ldots, \|A_{:,n}\|_2)$ is a __diagonal scaling__ , and $B$ is the __scaled matrix__ of $A$ from the right.  Then $[B^T B]_{i,i}=1$.

## Facts

1. Let $\tilde U$, $\tilde V$ and $\tilde \Sigma$ be the approximations of $U$, $V$ and $\Sigma$, respectively, computed by a backward stable algorithm as $A+\Delta A=\tilde U\tilde \Sigma \tilde V^T$. Since the orthogonality of $\tilde U$ and $\tilde V$ cannot be guaranteed, this product in general does not represent and SVD. There exist nearby orthogonal matrices $\hat U$ and $\hat V$ such that $(I+E_1)(A+\Delta A)(I+E_2)=\hat U \tilde \Sigma \hat V^T$, where departures from orthogonalithy, $E_1$ and $E_2$, are small in norm.

2. Standard algorithms compute the singular values with backward error $\| \Delta A\|\leq \phi\varepsilon \|A\|_2$, where $\varepsilon$ is machine precision and $\phi$ is a slowly growing function of $n$. The best error bound for the singular values is $|\sigma_j-\tilde \sigma_j|\leq \| \Delta A\|_2$, and the best relative error bound is

$$
\max_j \frac{|\sigma_j-\tilde\sigma_j|}{\sigma_j}\leq \frac{\| \Delta A \|_2}{\sigma_j} \leq \phi \varepsilon \kappa_2(A).$$

3. Let $\|[\Delta A]_{:,j}\|_2\leq \varepsilon \| A_{:,j}\|_2$ for all $j$. Then $A+\Delta A=(B+\Delta B)D$ and $\|\Delta B\|_F\leq \sqrt{n}\varepsilon$, and

$$
\max_j \frac{|\sigma_j-\tilde\sigma_j|}{\sigma_j}\leq
\| (\Delta B) B^{\dagger} \|_2\leq
\sqrt{n}\varepsilon \| B^{\dagger}\|_2.$$

This is Fact 3 from the [Relative perturbation theory](https://ivanslapnicar.github.io/GIAN-Applied-NLA-Course/L5b%20Singular%20Value%20Decomposition%20-%20Perturbation%20Theory%20.jl.html).

4. It holds

$$
\| B^{\dagger} \| \leq \kappa_2(B) \leq \sqrt{n} \min_{S=\mathop{\mathrm{diag}}}
\kappa_2(A\cdot S)\leq \sqrt{n}\,\kappa_2(A).$$

Therefore, numerical algorithm with column-wise small backward error computes singular values more accurately than an algorithm with small norm-wise backward error.

5. In each step, one-sided Jacobi method computes the Jacobi rotation matrix from the pivot submatrix of the current Gram matrix $A^TA$. Afterwards, $A$ is multiplied with the computed rotation matrix from the right (only two columns are affected). Convergence of the Jacobi method for the symmetric matrix $A^TA$ to a diagonal matrix, implies that the matrix $A$ converges to the matrix $AV$ with orthogonal columns and $V^TV=I$. Then $AV=U\Sigma$, $\Sigma=\mathop{\mathrm{diag}}(\| A_{:,1}\|_2, \ldots, \| A_{:,n}\|_2)$, $U=AV\Sigma^{-1}$, and  $A=U\Sigma V^T$ is the SVD of $A$.

6. One-sided Jacobi method computes the SVD with error bound from Facts 2 and 3, provided that the condition of the intermittent scaled matrices does not grow much. There is overwhelming numerical evidence for this. Alternatively, if $A$ is square, the one-sided Jacobi method can be applied to the transposed matrix $A^T=DB^T$ and the same error bounds apply, but the condition of the scaled matrix  (_this time from the left_) does not change. This approach is slower.

7. One-sided Jacobi method can be preconditioned by applying one QR factorization with full pivoting and one QR factorization withour pivoting to $A$, to obtain faster convergence, without sacrifying accuracy. This method is implemented in the LAPACK routine [DGESVJ](http://www.netlib.org/lapack/explore-html-3.3.1/d1/d5e/dgesvj_8f_source.html). _Writing the wrapper for `DGESVJ` is a tutorial assignment._

"""

# ╔═╡ 21541d5a-d787-438f-ada3-56e07e62ae4e
md"
## Examples

### Random matrix
"

# ╔═╡ 9fd7b6a0-ed2f-486b-99ec-806c1b19f033
function JacobiR(A₁::AbstractMatrix)
    A=deepcopy(A₁)
    m,n=size(A)
    T=typeof(A[1,1])
    V=Matrix{T}(I,n,n)
    # Tolerance for rotation
    tol=√n*eps(T)
    # Counters
    p=n*(n-1)/2
    sweep=0
    pcurrent=0
    # First criterion is for standard accuracy, second one is for relative accuracy
    # while sweep<30 && vecnorm(A-diagm(diag(A)))>tol
    while sweep<30 && pcurrent<p
        sweep+=1
        # Row-cyclic strategy
        for i = 1 : n-1
            for j = i+1 : n
                # Compute the 2 x 2 sumbatrix of A'*A
                # F=A[:,[i,j]]'*A[:,[i,j]]
				F=view(A,:,[i,j])'*view(A,:,[i,j])
                # Check the tolerance - the first criterion is standard,
                # the second one is for relative accuracy
                # if A[i,j]!=zero(T)
                #
                if abs(F[1,2])>tol*sqrt(F[1,1]*F[2,2])
                    # Compute c and s
                    τ=(F[2,2]-F[1,1])/(2*F[1,2])
                    t=sign(τ)/(abs(τ)+√(1+τ^2))
                    c=1/sqrt(1+t^2)
                    s=c*t
                    G=LinearAlgebra.Givens(i,j,c,s)
                    # A*=G'
                    rmul!(A,G)
                    # V*=G'
                    rmul!(V,G)
                    pcurrent=0
                else
                    pcurrent+=1
                end
            end
        end
    end
    σ=[norm(A[:,k]) for k=1:n]
    for k=1:n
        A[:,k]./=σ[k]
    end
    # A, σ, V
    SVD(A,σ,adjoint(V))
end

# ╔═╡ 04b83647-fde5-40ad-a0e1-aebfd99c40b5
begin
	m=8
	n=5
	import Random
	Random.seed!(432)
	A=rand(-9.0:9,m,n)
end

# ╔═╡ 751c0242-28cf-4e85-a28c-adc2b096511d
U,σ,V=JacobiR(A)

# ╔═╡ 3b813d5b-b27f-41bd-8c28-b7a4f3591ca8
# Residual and orthogonality
norm(A*V-U*Diagonal(σ)), norm(U'*U-I),norm(V'*V-I)

# ╔═╡ 25b11022-ddd4-48c7-91a8-8306c8f04d08
md"""
### Strongly scaled matrix
"""

# ╔═╡ 90059aec-9467-4548-9ef3-05a36d8e730e
begin
	m₁=20
	n₁=15
	B₁=rand(m₁,n₁)
	D₁=exp.(50*(rand(n₁).-0.5))
	A₁=B₁*Diagonal(D₁)
end

# ╔═╡ d71f46c3-0265-4e9b-9e11-bad11ae617ca
cond(B₁), cond(A₁)

# ╔═╡ 9b015091-c8b6-416f-bfa0-cd8742feae53
U₁,σ₁,V₁=JacobiR(A₁);

# ╔═╡ fb528048-4886-476d-acf6-bb0d1105f3c3
[sort(σ₁,rev=true) svdvals(A₁)]

# ╔═╡ 2c628908-f9e2-4b04-a7b8-ebe6ff97d207
# Relative errors
(sort(σ₁,rev=true)-svdvals(A₁))./sort(σ₁,rev=true)

# ╔═╡ f1734a00-d452-4adf-8125-dfbcee74d338
# Residual
norm(A₁*V₁-U₁*Diagonal(σ₁))/norm(A₁)

# ╔═╡ caa171b1-e8cd-4fbe-a064-964c1a896f7a
# Diagonalization
U₁'*A₁*V₁

# ╔═╡ 98fa4641-2d3f-48f5-ae1d-a9f3b96f798e
md"""
In the alternative approach, we first apply QR factorization with column pivoting to obtain the square matrix.
"""

# ╔═╡ 19f3f447-ca93-4a6a-8b3f-403afd273fdd
# ?qr

# ╔═╡ 9dc37b20-3b94-4cbe-9c96-0bf12bfbe005
Q=qr(A₁,Val(true))

# ╔═╡ 38954e60-3f69-407e-8000-8a6e4477c245
# Residual
norm(Q.Q*Q.R-A₁[:,Q.p])

# ╔═╡ bda5ca94-e567-4615-a3f4-bce3c88cfb39
S=JacobiR(Q.R')

# ╔═╡ aeb0105f-53e6-437b-a1f0-376ffcfdb5d7
(sort(σ₁)-sort(S.S))./sort(σ₁)

# ╔═╡ 7c3ffb94-6403-43da-81fa-6fd8ef8a24b8
md"
Now $QRP^T=A$ and $R^T=U_R\Sigma_R V^T_R$, so

$$
A=(Q V_R) \Sigma_R (U_R^T P^T)$$

is an SVD of $A$.
"

# ╔═╡ 11d8a586-361f-42f4-bd9a-77e483fc7a74
# Residual
norm(A₁*S.U[invperm(Q.p),:]-Q.Q*S.V*Diagonal(S.S))

# ╔═╡ e1dffb13-ed1b-44b4-9cb1-bd5227962cc9
md"
# Lanczos method

The function `svds()` is based on the Lanczos method for symmetric matrices. Input can be matrix, but also an operator which defines the product of the given matrix with a vector.

## Examples

### Random matrix
"

# ╔═╡ e6f4eb1b-2a3c-4785-8645-df5c72aa1ad7
# ?svds

# ╔═╡ 8c66ea99-3a59-4923-b73c-4ae3c9d7062e
begin
	m₂=20
	n₂=15
	A₂=rand(m₂,n₂)
end

# ╔═╡ 8af3203b-9e58-460a-a7e0-81c60d6dcad4
σ₂=svdvals(A₂)

# ╔═╡ 5f2c6284-21c4-4bdc-8545-50bf4944f1ef
begin
	# Some largest singular values
	k₂=6
	S₂,rest=svds(A₂,nsv=k₂)
	(σ₂[1:k₂]-S₂.S)./σ₂[1:k₂]
end

# ╔═╡ 559e2a09-70d4-4d1e-8941-45ed9ccbae6b
md"
### Large matrix
"

# ╔═╡ 91952aea-a239-4e2d-80f1-df6aefc58c4c
begin
	m₃=2000
	n₃=1500
	A₃=randn(m₃,n₃)
end

# ╔═╡ a7068313-39ae-4e83-a5cd-68aec37712c9
@time U₃,σ₃,V₃=svd(A₃);

# ╔═╡ 2e1d70eb-9751-4f4e-826f-e2aa05612cca
begin
	# This is rather slow
	k₃=10
	@time S₃,rest₃=svds(A₃,nsv=k₃);
end

# ╔═╡ 11a45bbe-a9f9-44e5-a7c2-559bddebefae
(σ₃[1:k₃]-S₃.S)./σ₃[1:k₃]

# ╔═╡ 5a75f57f-1293-4c76-9095-4d32e00cc093
md"""
### Very large sparse matrix
"""

# ╔═╡ 37f296ac-8853-480d-b424-b9037a7e50bc
#$ ?sprand

# ╔═╡ 151e13fb-3c48-40fb-8699-17ebb1569f2d
A₄=sprand(10000,3000,0.05)

# ╔═╡ 017a66c9-bbc2-464a-905b-b5316c49292a
begin
	# No vectors, this takes about 5 sec.
	k₄=10
	@time S₄,rest₄=svds(A₄,nsv=k₄,ritzvec=false)
end

# ╔═╡ 435d9693-996e-4f35-ac8d-e1f70077ce94
# Full matrix, no vectors, about 19 sec
@time σ₄=svdvals(Matrix(A₄))

# ╔═╡ 1b45c928-f1a4-437d-9980-a70a8e1896db
maximum(abs,(S₄.S-σ₄[1:k₄])./σ₄[1:k₄])

# ╔═╡ Cell order:
# ╟─a93d84c9-400f-4624-a67d-cec990f3d822
# ╟─f04cb86f-94aa-418d-82f9-d9ef12ca90fd
# ╟─2eff5f1f-0989-4292-8f04-80983c1bd133
# ╟─21541d5a-d787-438f-ada3-56e07e62ae4e
# ╠═9fd7b6a0-ed2f-486b-99ec-806c1b19f033
# ╠═04b83647-fde5-40ad-a0e1-aebfd99c40b5
# ╠═751c0242-28cf-4e85-a28c-adc2b096511d
# ╠═3b813d5b-b27f-41bd-8c28-b7a4f3591ca8
# ╟─25b11022-ddd4-48c7-91a8-8306c8f04d08
# ╠═90059aec-9467-4548-9ef3-05a36d8e730e
# ╠═d71f46c3-0265-4e9b-9e11-bad11ae617ca
# ╠═9b015091-c8b6-416f-bfa0-cd8742feae53
# ╠═fb528048-4886-476d-acf6-bb0d1105f3c3
# ╠═2c628908-f9e2-4b04-a7b8-ebe6ff97d207
# ╠═f1734a00-d452-4adf-8125-dfbcee74d338
# ╠═caa171b1-e8cd-4fbe-a064-964c1a896f7a
# ╟─98fa4641-2d3f-48f5-ae1d-a9f3b96f798e
# ╠═19f3f447-ca93-4a6a-8b3f-403afd273fdd
# ╠═9dc37b20-3b94-4cbe-9c96-0bf12bfbe005
# ╠═38954e60-3f69-407e-8000-8a6e4477c245
# ╠═bda5ca94-e567-4615-a3f4-bce3c88cfb39
# ╠═aeb0105f-53e6-437b-a1f0-376ffcfdb5d7
# ╟─7c3ffb94-6403-43da-81fa-6fd8ef8a24b8
# ╠═11d8a586-361f-42f4-bd9a-77e483fc7a74
# ╟─e1dffb13-ed1b-44b4-9cb1-bd5227962cc9
# ╠═afe8da5c-4e4e-4ce7-87de-e3b201a4c133
# ╠═e6f4eb1b-2a3c-4785-8645-df5c72aa1ad7
# ╠═8c66ea99-3a59-4923-b73c-4ae3c9d7062e
# ╠═8af3203b-9e58-460a-a7e0-81c60d6dcad4
# ╠═5f2c6284-21c4-4bdc-8545-50bf4944f1ef
# ╟─559e2a09-70d4-4d1e-8941-45ed9ccbae6b
# ╠═91952aea-a239-4e2d-80f1-df6aefc58c4c
# ╠═a7068313-39ae-4e83-a5cd-68aec37712c9
# ╠═2e1d70eb-9751-4f4e-826f-e2aa05612cca
# ╠═11a45bbe-a9f9-44e5-a7c2-559bddebefae
# ╟─5a75f57f-1293-4c76-9095-4d32e00cc093
# ╠═8db1c1b2-8b5f-402c-8441-768c1a943c38
# ╠═37f296ac-8853-480d-b424-b9037a7e50bc
# ╠═151e13fb-3c48-40fb-8699-17ebb1569f2d
# ╠═017a66c9-bbc2-464a-905b-b5316c49292a
# ╠═435d9693-996e-4f35-ac8d-e1f70077ce94
# ╠═1b45c928-f1a4-437d-9980-a70a8e1896db

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L6b Singular Value Decomposition - Jacobi and Lanczos Methods.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
#############################################################################
# gradient.jl
# calculates the gradients of functions
#############################################################################

export grad
export sign, curvature, monotonicity, evaluate

mutable struct GradientAtom <: AbstractExpr
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr}
  size::Tuple{Int, Int}
  variable::Variable
  func::Function

  function GradientAtom(x::AbstractExpr)
    var = get_cv(x)
    if x.size != (1, 1) || sign(x) == ComplexSign()
      error("A function has to be real-valued to calculate its gradient.")
    elseif length(var) > 1
      error("Only able to calculate gradients of functions with single variable.")
    elseif var[1].size[2] != 1
      error("Only able to calculate gradients of functions with variable's size being (n, 1).")
    end
    children = (x,)
    f(z) = evaluate(x, z)[1]
    return new(:grad, hash(children), children, var[1].size, var[1], f)
  end
end

grad(x::AbstractExpr) = GradientAtom(x)

function sign(x::GradientAtom)
  if monotonicity(x.children[1]) == (Nonincreasing(), )
    return Positive()
  elseif monotonicity(x.children[1]) == (Nondecreasing(), )
    return Negative()
  elseif monotonicity(x.children[1]) == (ConstMonotonicity(), )
    return Positive()
  else
    return NoSign()
  end
end

function curvature(x::GradientAtom)
  return NotDcp()
end

function monotonicity(x::GradientAtom)
  if curvature(x.children[1]) == ConvexVexity()
    return (Nondecreasing(), )
  elseif curvature(x.children[1]) == ConcaveVexity()
    return (Nonincreasing(), )
  elseif curvature(x.children[1]) == AffineVexity
    return (ConstMonotonicity(), )
  else
    return (NoMonotonicity(), )
  end
end

function evaluate(x::GradientAtom)
  var = copy(x.variable.value)
  val = gradient(x.func, var)
  x.variable.value = var
  return val
end


text:
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using Plots

@parallel function compute_τ!(Vx, Vy, Vz, ∇V, τxx, τyy, τzz, τxy, τyz, τzx, dt1_3, dtμ_dx, dtμ_dy, dtμ_dz, dt2μ_dx, dt2μ_dy, dt2μ_dz)
    @all(τxx) = @all(τxx) + @d_xi(Vx)*dt2μ_dx - dt1_3*@inn_yz(∇V)
    @all(τyy) = @all(τyy) + @d_yi(Vy)*dt2μ_dy - dt1_3*@inn_xz(∇V)
    @all(τzz) = @all(τzz) + @d_zi(Vz)*dt2μ_dz - dt1_3*@inn_xy(∇V)
    @all(τxy) = @all(τxy) + @d_yi(Vx)*dtμ_dy + @d_xi(Vy)*dtμ_dx
    @all(τyz) = @all(τyz) + @d_zi(Vy)*dtμ_dz + @d_yi(Vz)*dtμ_dy
    @all(τzx) = @all(τzx) + @d_xi(Vz)*dtμ_dx + @d_zi(Vx)*dtμ_dz
    return
end

@parallel function compute_V!(P, Vx, Vy, Vz, ∇V, τxx, τyy, τzz, τxy, τyz, τzx, dt_ρ_dx, dt_ρ_dy, dt_ρ_dz)
	# dVxdt = -1/ρ*(@d_xi(P)*_dx - @d_xa(τxx)*_dx - @d_ya(τxy)*_dy - @d_za(τzx)*_dz)
	# Vx    = Vx + dt*dVxdt
    @inn(Vx) = @inn(Vx) - (@d_xi(P) - @d_xa(τxx))*dt_ρ_dx + @d_ya(τxy)*dt_ρ_dy + @d_za(τzx)*dt_ρ_dz
    @inn(Vy) = @inn(Vy) - (@d_yi(P) - @d_ya(τyy))*dt_ρ_dy + @d_za(τyz)*dt_ρ_dz + @d_xa(τxy)*dt_ρ_dx
    @inn(Vz) = @inn(Vz) - (@d_zi(P) - @d_za(τzz))*dt_ρ_dz + @d_xa(τzx)*dt_ρ_dx + @d_ya(τyz)*dt_ρ_dy
    return
end

@parallel function compute_∇V!(Vx, Vy, Vz, ∇V, _dx, _dy, _dz)
    @all(∇V) = @d_xa(Vx)*_dx + @d_ya(Vy)*_dy + @d_za(Vz)*_dz
    return
end

@parallel function compute_P!(P, ∇V, K, dt)
	# dPdt = -K*∇V
	# P    = P + dt*dPdt
    @all(P) = @all(P) - dt*K*@all(∇V)
    return
end

"""
Runs the 3D elastic wave simulation\\
`elastic_wave_3D(res, ttot) -> xc, P`\\
    `res`    (Int)     : xyz-resolutions of the simulation\\
    `ttot`   (Float64) : total simulation time\\
    `xc` and `P` are returned for testing purposes
"""
@views function elastic_wave_3D(res::Int, ttot::Float64)
    # Physics
    Lx, Ly, Lz = 10.0, 10.0, 10.0
    ρ          = 1.0
    μ          = 1.0
    K          = 1.0

    # Numerics
    nx, ny, nz = res, res, res

    # Derived numerics
    dx, dy, dz = Lx/nx, Ly/ny, Lz/nz
    dt         = min(dx,dy,dz)/sqrt((K + 4/3*μ)/ρ)/2.1
    nt         = cld(ttot, dt)
    xc, yc, zc = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny), LinRange(dz/2, Lz-dz/2, nz)

	# Optimized numerics
	_dx, _dy, _dz = 1.0/dx, 1.0/dy, 1.0/dz
    dtμ_dx, dtμ_dy, dtμ_dz = dt*μ*_dx, dt*μ*_dy, dt*μ*_dz
    dt2μ_dx, dt2μ_dy, dt2μ_dz = 2.0*dtμ_dx, 2.0*dtμ_dy, 2.0*dtμ_dz
	dt_ρ = dt/ρ
    dt_ρ_dx, dt_ρ_dy, dt_ρ_dz = dt_ρ*_dx, dt_ρ*_dy, dt_ρ*_dz
    dt1_3 = dt/3.0

    # Array initialisation
    P   = Data.Array(exp.([-1.0*((x-Lx/2)^2 + (y-Ly/2)^2 + (z-Lz/2)^2) for x=xc, y=yc, z=zc]))
    Vx  = @zeros(nx+1,ny  ,nz  )
    Vy  = @zeros(nx  ,ny+1,nz  )
    Vz  = @zeros(nx  ,ny  ,nz+1)
    ∇V  = @zeros(nx  ,ny  ,nz  )
    τxx = @zeros(nx  ,ny-2,nz-2)
    τyy = @zeros(nx-2,ny  ,nz-2)
    τzz = @zeros(nx-2,ny-2,nz  )
    τxy = @zeros(nx-1,ny-1,nz-2)
    τyz = @zeros(nx-2,ny-1,nz-1)
    τzx = @zeros(nx-1,ny-2,nz-1)

    # Benchmark variables
    @static if BENCHMARK
        warmup = 10     # Number of warmup-iterations
        t_tic  = 0.0    # Holds start time
    end

    # Animation object for visualization
    @static if VISUALIZE
        anim = Animation()
    end

    # Time loop
    for it = 1:nt
        @parallel compute_τ!(Vx, Vy, Vz, ∇V, τxx, τyy, τzz, τxy, τyz, τzx, dt1_3, dtμ_dx, dtμ_dy, dtμ_dz, dt2μ_dx, dt2μ_dy, dt2μ_dz)
        @parallel compute_V!(P, Vx, Vy, Vz, ∇V, τxx, τyy, τzz, τxy, τyz, τzx, dt_ρ_dx, dt_ρ_dy, dt_ρ_dz)
        @parallel compute_∇V!(Vx, Vy, Vz, ∇V, _dx, _dy, _dz)
        @parallel compute_P!(P, ∇V, K, dt)

        # Start the clock after warmup-iterations
        @static if BENCHMARK
            if (it == warmup+1)
                t_tic = Base.time();
            end
        end

        # Render a slice of the 3D pressure-map and save as an animation frame
        @static if VISUALIZE
            opts = (aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), clims=(-0.2, 0.6), c=:davos, xlabel="Lx", ylabel="Ly", title="time = $(round(it*dt, sigdigits=3))")
            heatmap(xc, yc, P[:,:,cld(nz,2)]'; opts...)
            frame(anim)
        end
    end

    # Calculate performance
    @static if BENCHMARK
        t_toc = Base.time() - t_tic                 # Stop the clock
        A_eff = 1e-9 * sizeof(Float64) * (          # Effective main memory access per iteration [GB]
              2 * (length(τxx) + length(τyy) + length(τzz) + length(τxy) + length(τyz) + length(τzx))
            + 2 * (length(Vx)  + length(Vy)  + length(Vz))
            + 2 * (length(P))
        )
        t_it  = t_toc/(nt-warmup)                   # Execution time per iteration [s]
        T_eff = A_eff/t_it                          # Effective memory throughput [GB/s]
        println("size=$res, t=$t_toc s, T_eff=$T_eff GB/s")
    end

    # Save the animation as a GIF
    @static if VISUALIZE
		fps = min(60.0, (float(nt)/ttot))
        gif(anim, "./docs/img/elastic_wave_3D.gif", fps=fps)
    end

	@static if BENCHMARK
		return T_eff
	else
		return xc, P
	end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""

    RadioMenu

A menu that allows a user to select a single option from a list.

# Sample Output

```julia
julia> request(RadioMenu(options, pagesize=4))
Choose your favorite fruit:
^  grape
   strawberry
 > blueberry
v  peach
Your favorite fruit is blueberry!
```

"""
mutable struct RadioMenu{C} <: _ConfiguredMenu{C}
    options::Array{String,1}
    keybindings::Vector{Char}
    pagesize::Int
    pageoffset::Int
    selected::Int
    config::C
end

"""

    RadioMenu(options::Array{String,1}; pagesize::Int=10, kwargs...)

Create a RadioMenu object. Use `request(menu::RadioMenu)` to get user input.
`request()` returns an `Int` which is the index of the option selected by the
user.

# Arguments

  - `options::Array{String, 1}`: Options to be displayed
  - `pagesize::Int=10`: The number of options to be displayed at one time, the menu will scroll if length(options) > pagesize

Any additional keyword arguments will be passed to [`TerminalMenus.Config`](@ref).
"""
function RadioMenu(options::Array{String,1}; pagesize::Int=10, warn::Bool=true, keybindings::Vector{Char}=Char[], kwargs...)
    length(options) < 1 && error("RadioMenu must have at least one option")
    length(keybindings) in [0, length(options)] || error("RadioMenu must have either no keybindings, or one per option")

    # if pagesize is -1, use automatic paging
    pagesize = pagesize == -1 ? length(options) : pagesize
    # pagesize shouldn't be bigger than options
    pagesize = min(length(options), pagesize)
    # after other checks, pagesize must be greater than 1
    pagesize < 1 && error("pagesize must be >= 1")

    pageoffset = 0
    selected = -1 # none

    if !isempty(kwargs)
        RadioMenu(options, keybindings, pagesize, pageoffset, selected, Config(; kwargs...))
    else
        warn && Base.depwarn("Legacy `RadioMenu` interface is deprecated, set a configuration option such as `RadioMenu(options; charset=:ascii)` to trigger the new interface.", :RadioMenu)
        RadioMenu(options, keybindings, pagesize, pageoffset, selected, CONFIG)
    end
end

# AbstractMenu implementation functions
# See AbstractMenu.jl
#######################################

options(m::RadioMenu) = m.options

cancel(m::RadioMenu) = m.selected = -1

function pick(menu::RadioMenu, cursor::Int)
    menu.selected = cursor
    return true #break out of the menu
end

function writeline(buf::IOBuffer, menu::RadioMenu{Config}, idx::Int, iscursor::Bool)
    print(buf, replace(menu.options[idx], "\n" => "\\n"))
end

function keypress(m::RadioMenu, i::UInt32)
    isempty(m.keybindings) && return false
    i = findfirst(isequal(i), Int.(m.keybindings))
    isnothing(i) && return false
    m.selected = i
    return true
end

# Legacy interface
function writeLine(buf::IOBuffer, menu::RadioMenu{<:Dict}, idx::Int, cursor::Bool)
    # print a ">" on the selected entry
    cursor ? print(buf, menu.config[:cursor] ," ") : print(buf, "  ")

    print(buf, replace(menu.options[idx], "\n" => "\\n"))
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    abstract type FlipSignBridge{T,S1,S2} <: SetMapBridge{T,S2,S1} end

An abstract type that simplifies the creation of other bridges.
"""
abstract type FlipSignBridge{T,S1<:MOI.AbstractSet,S2<:MOI.AbstractSet} <:
              SetMapBridge{T,S2,S1} end

function MOI.Bridges.map_function(::Type{<:FlipSignBridge{T}}, func) where {T}
    return MOI.Utilities.operate(-, T, func)
end

# The map is an involution
function MOI.Bridges.inverse_map_function(BT::Type{<:FlipSignBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is symmetric
function MOI.Bridges.adjoint_map_function(BT::Type{<:FlipSignBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is a symmetric involution
function MOI.Bridges.inverse_adjoint_map_function(
    BT::Type{<:FlipSignBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end

function MOI.delete(
    model::MOI.ModelLike,
    bridge::FlipSignBridge,
    i::MOI.Bridges.IndexInVector,
)
    MOI.delete(model, bridge.variables[i.value])
    deleteat!(bridge.variables, i.value)
    return
end

"""
    NonposToNonnegBridge{T} <: Bridges.Variable.AbstractBridge

`NonposToNonnegBridge` implements the following reformulation:

* ``x \\in \\mathbb{R}_-`` into ``y \\in \\mathbb{R}_+`` with the substitution
  rule ``x = -y``,

where `T` is the coefficient type of `-y`.

## Source node

`NonposToNonnegBridge` supports:

* [`MOI.VectorOfVariables`](@ref) in [`MOI.Nonpositives`](@ref)

## Target nodes

`NonposToNonnegBridge` creates:

* One variable node: [`MOI.VectorOfVariables`](@ref) in [`MOI.Nonnegatives`](@ref),
"""
struct NonposToNonnegBridge{T} <:
       FlipSignBridge{T,MOI.Nonpositives,MOI.Nonnegatives}
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}
end

const NonposToNonneg{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonposToNonnegBridge{T},OT}

function MOI.Bridges.map_set(
    ::Type{<:NonposToNonnegBridge},
    set::MOI.Nonnegatives,
)
    return MOI.Nonpositives(set.dimension)
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:NonposToNonnegBridge},
    set::MOI.Nonpositives,
)
    return MOI.Nonnegatives(set.dimension)
end


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

emit(identity, Bool)

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Base.identity), Bool}"(%arg0: !jlir<"typeof(Base.identity)">, %arg1: i1) -> i1 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     return %arg1 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Base.identity), Bool}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i1) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     llvm.return %arg1 : i1
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Base.identity), Bool}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i1) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Base.identity), Bool}"(%arg0, %arg1) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i1) -> i1
# CHECK-NEXT:     llvm.return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
module Examples

using Compat

function small()
  P = cat(3, [0.5 0.8; 0.5 0.2], [0 0.1; 1 0.9])
  R = [5 10; -1 2]
  (P, R)
end

include("random.jl")
include("sparse_random.jl")

# function bellmanmachine(k::Int, M::Int, p::Array{Float64,2}, q::Array{Float64,2})
#     #
# end

end # module


text:
ENV["MPLBACKEND"]="agg"
using Documenter, ExtendableGrids, Literate, GridVisualize
import PyPlot

example_md_dir  = joinpath(@__DIR__,"src","examples")

examples1d=joinpath(@__DIR__,"..","examples","examples1d.jl")
include(examples1d)
examples2d=joinpath(@__DIR__,"..","examples","examples2d.jl")
include(examples2d)
examples3d=joinpath(@__DIR__,"..","examples","examples3d.jl")
include(examples3d)

include("makeplots.jl")

function mkdocs()

    Literate.markdown(examples1d, example_md_dir, documenter=false,info=false)
    Literate.markdown(examples2d, example_md_dir, documenter=false,info=false)
    Literate.markdown(examples3d, example_md_dir, documenter=false,info=false)

    generated_examples=joinpath.("examples",filter(x->endswith(x, ".md"),readdir(example_md_dir)))

    makeplots(example_md_dir)

    makedocs(sitename="ExtendableGrids.jl",
    modules = [ExtendableGrids],
    doctest = false,
    clean = true,
             authors = "J. Fuhrmann, Ch. Merdon",
             repo="https://github.com/j-fu/ExtendableGrids.jl",
             pages=[
                 "Home"=>"index.md",
                 "adjacency.md",
                 "vectorofconstants.md",
                 "typehierarchy.md",
                 "elementgeometry.md",
                 "coordinatesystem.md",
                 "extendablegrid.md",
                 "subgrid.md",
                 "regionedit.md",
                 "simplexgrid.md",
                 "more.md",
                 "tokenstream.md",
                 "allindex.md",
                 "Examples" => generated_examples
             ])
end

mkdocs()

deploydocs(repo = "github.com/j-fu/ExtendableGrids.jl.git")


text:
using JetPack, Jets, Test

@testset "JopProjection, 2D, first component" begin
    u = zeros(256,2)
    u[:,1] .= 1.0
    u[:,2] .= 0.0

    P = JopProjection(u)
    m = zeros(256,2)
    m[:,1] .= 1.0
    m[:,2] .= 2.0

    d = P*m
    a = P'*d

    d_expected = 1.0*ones(256)
    a_expected = zeros(256,2)
    a_expected[:,1] .= 1.0

    @test d ≈ 1.0*ones(256)
    @test a ≈ a_expected

    lhs,rhs = dot_product_test(P, rand(domain(P)), rand(range(P)))
    @test lhs ≈ rhs
end

@testset "JopProjection, 2D, second component" begin
    u = zeros(256,2)
    u[:,1] .= 0.0
    u[:,2] .= 1.0

    P = JopProjection(u)
    m = zeros(256,2)
    m[:,1] .= 1.0
    m[:,2] .= 2.0

    d = P*m
    a = P'*d

    d_expected = 2.0*ones(256)
    a_expected = zeros(256,2)
    a_expected[:,2] .= 2.0

    @test d ≈ 2.0*ones(256)
    @test a ≈ a_expected

    lhs,rhs = dot_product_test(P, rand(domain(P)), rand(range(P)))
    @test lhs ≈ rhs
end


text:
# # Numerical Stability

using Distributions
using HMMBase
using PyPlot
using Seaborn

rc("axes", xmargin = 0) # hide
set_style("whitegrid")  # hide

# Let's consider the case of a Normal distribution with null variance,
# such a case can appear during maximum likelihood estimation if only
# one observation is associated to a state:

A = [0.99 0.01; 0.1 0.9]
B = [Normal(0, 1), Normal(10.5, 0)]
hmm = HMM(A, B)
#-

y = rand(hmm, 500)
figure(figsize = (9,2)) # hide
plot(y)
gcf() # hide
#-

# The likelihood of a Normal distribution with null variance goes to infinity for `y = μ`,
# as there is a division by zero in the density function:
println(extrema(likelihoods(hmm, y)))
println(extrema(likelihoods(hmm, y, logl = true)))

# To avoid propagating these non-finite quantities (for example in the forward-backward algorithm),
# you can use the `robust` option:
println(extrema(likelihoods(hmm, y, robust = true)))
println(extrema(likelihoods(hmm, y, logl = true, robust = true)))

# This truncates `+Inf` to the largest Float64, and `-Inf` to the smallest Float64:
prevfloat(Inf), nextfloat(-Inf)

# In the log. case we use `log(prevfloat(Inf))` to avoid overflows when taking
# the exp. of the log-likelihood.

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
module PlutoEducation

function hello(name="world")
    println("Hello, ", name, "!")
end

end # module PlutoEducation


text:
# This file contains functions used in the high level service definitions
_merge(a::AbstractDict, b::AbstractDict) = merge(a, b)
_merge(a::AbstractString, b::AbstractString) = b

function _flatten_query(service::String, query::AbstractDict{String, <: Any})
    return _flatten_query!(Pair{String, String}[], service, query)
end

function _flatten_query!(result::Vector{Pair{String, String}}, service::String, query::AbstractDict{String, <:Any}, prefix::String="")
    for (k, v) in query
        if v isa AbstractDict
            _flatten_query!(result, service, v, string(prefix, k, "."))
        elseif v isa AbstractArray
            for (i, j) in enumerate(v)
                suffix = service in ("ec2", "sqs") ? "" : ".member"
                prefix_key = string(prefix, k, suffix, ".", i)

                if j isa AbstractDict
                    _flatten_query!(result, service, j, string(prefix_key, "."))
                else
                    push!(result, Pair(prefix_key, string(j)))
                end
            end
        else
            push!(result, Pair(string(prefix, k), string(v)))
        end
    end

    return result
end

"""
Escape special AWS S3 characters in the path of the uri properly.

AWS S3 allows for various special characters in file names, these characters are not being
properly escaped before we make the requests.

We cannot call `HTTP.escapeuri(request.uri)` because this will escape `/` characters which
are used in the filepathing for sub-directories.
"""
function _clean_s3_uri(uri::AbstractString)
    chars_to_clean = (
        ' ' => "%20",
        '!' => "%21",
        ''' => "%27",
        '(' => "%28",
        ')' => "%29",
        '*' => "%2A",
        '+' => "%2B",
        ',' => "%2C",
        ':' => "%3A",
        '=' => "%3D",
        '@' => "%40",
    )
    parsed_uri = URIs.URI(uri)
    cleaned_path = reduce(replace, chars_to_clean, init=parsed_uri.path)
    return string(URIs.URI(parsed_uri; path=cleaned_path))
end

function _extract_common_kw_args(service, args)
    return (
        service=service.name,
        api_version=service.api_version,
        return_stream=_pop!(args, "return_stream", false),
        return_raw=_pop!(args, "return_raw", false),
        response_stream=_pop!(args, "response_stream", nothing),
        headers=LittleDict{String, String}(_pop!(args, "headers", [])),
        http_options=_pop!(args, "http_options", LittleDict{Symbol, String}()),
        response_dict_type=_pop!(args, "response_dict_type", LittleDict),
        backend=_pop!(args, "backend", DEFAULT_BACKEND[]),
    )
end

# Use this until the three arg pop! is available for LittleDict
# https://github.com/JuliaCollections/OrderedCollections.jl/pull/59
function _pop!(dict::AbstractDict{String, <: Any}, kw, default)
    if haskey(dict,kw)
        val = dict[kw]
        delete!(dict, kw)

        return val
    else
        return default
    end
end

function _generate_rest_resource(request_uri::String, args::AbstractDict{String, <:Any})
    # There maybe a time where both $k and $k+ are in the request_uri, in which case this needs to be updated
    # From looking around, I have not seen an example yet, however that doesn't mean it doesn't exist

    for (k,v) in args
        if occursin("{$k}", request_uri)
            request_uri = replace(request_uri, "{$k}" => v)
        elseif occursin("{$k+}", request_uri)
            request_uri = replace(request_uri, "{$k+}" => HTTP.escapepath(v))
        end
    end

    return request_uri
end


text:
using Test

"
Function shootingmethod() for 2nd order BVP ODEs.

Arguments to shootingmethod(f, bvs, initg, steps, tol, limit):

f::Function           : Function to be estimated (see below).\n
bvs::Vector{Float64}  : Vector of boundary values [xa, ya, xb, yb].\n
initg::Matrix{Float64}: Vector of nitial gradients [a0, a1].\n
steps::Int64          : Number of calculations steps.\n
tol::Float64          : Tolerance for convergence.\n
limit::Int64          : Max number of iterations.

Example of how to define the function f(x, y) for y'' = 3x^3+4y:

Decompose in two first order equations:

dy[1] = y[2]\n
dy[2] = 3x^3 + 4y[1]

```
function f(x::Float64, y::Vector{Float64})
  [y[2], 3x^3 +4y[1]]
end
```
"
function shootingmethod(f::Function, bvs::Vector{Float64}, initg::Vector{Float64},
  steps::Int64, tol::Float64=0.00001, limit::Int64=25)
  #println("\n--- Shooting method for Second Order ODEs ---\n")
  @test size(bvs, 1) == 4
  @test size(initg, 1) == 2
  nsteps = steps + 1
  y0 = zeros(Float64, nsteps, 2)
  ystar = zeros(Float64, nsteps)
  y = zeros(Float64, 2)
  (xa, ya, xb, yb) = bvs[:, 1]
  a0 = initg[:, 1]
  h = (xb - xa)/steps
  res = zeros(Float64, nsteps, 2)
  for j in 1:2
    x = xa
    y[1] = ya
    y[2] = a0[j]
    for i in 1:nsteps
      y0[i, j] = y[1]
      k0 = h * f(x, y)
      k1 = h * f(x+h/2, y+k0/2)
      k2 = h * f(x+h/2, y+k1/2)
      k3 = h * f(x+h, y+k2)
      y += (k0+2k1+2k2+k3)/6
      x = x + h
    end
  end
  if (y0[nsteps,1]-yb)*(y0[nsteps,2]-yb) > 0.0
    println("Try new gradients ... ?")
    converged = abs((ystar[nsteps]-yb)/yb) < tol
    return((0, nsteps, converged, res))
  end
  iter = 0
  while iter <= limit
    iter += 1
    if iter > limit
      println("Max number of iterations exceeded ...")
      println("Iterations performed: $iter")
      for i in 1:nsteps
        res[i,:] = [xa+(i-1)*h ystar[i]]
      end
      println("Try new gradients ... ?")
      converged = abs((ystar[nsteps]-yb)/yb) < tol
      return((iter, nsteps, converged, res))
    end
    astar = a0[1] + (yb-y0[nsteps, 1])*(a0[2]-a0[1])/(y0[nsteps,2])
    x = xa
    y[1] = ya
    y[2] = astar
    for i in 1:nsteps
      ystar[i] = y[1]
      k0 = h * f(x, y)
      k1 = h * f(x+h/2, y+k0/2)
      k2 = h * f(x+h/2, y+k1/2)
      k3 = h * f(x+h, y+k2)
      y += (k0+2k1+2k2+k3)/6
      x = x + h
    end
    if abs((ystar[nsteps]-yb)/yb) < tol
      #println("Iterations to convergence: $iter")
      for i in 1:nsteps
        res[i,:] = [xa+(i-1)*h ystar[i]]
      end
      break
    end
    if (ystar[nsteps]-yb)*(y0[nsteps, 1]-yb) > 0.0
      y0[:, 1] = ystar
      a0[1] = astar
    else
      y0[:, 2] = ystar
      a0[2] = astar
    end
  end
  converged = abs((ystar[nsteps]-yb)/yb) < tol
  (iter, nsteps, converged, res)
end