text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences
import BioSequences: iscompatible

################################################
#OligoRL Toolset
#################################################

# Overload the concatenation operator to combine a sequence and
# a single base, i.e. dna"AGCGTGC" * DNA_T
function *(seq::LongDNASeq, base::DNA)
    seq * LongDNASeq([base])
end

# Overload sort to order DNA codes by degeneracy:
#   A,G,C,T,M,R,W,S,Y,K,V,H,D,B,N
function sort(seq::LongSequence{DNAAlphabet{4}}; rev::Bool=false)
    sort(convert(Vector{DNA}, seq), by=degeneracy, rev=rev)
end

#Convert an array of DNA sequences to an array of strings
function dnaarray2stringarray(array)
    array2=[]
    for i=1:length(array)
        push!(array2,convert(String,array[i]))
    end
    return array2
end

#Convert an array of strings to an array of DNA sequences
function stringarray2dnaarray(array)
    array2=[]
    for i=1:length(array)
        push!(array2,LongDNASeq(array[i]))
    end
    return array2
end
# Determine if a given sequence is palindromic. Works both for unambiguous and ambigous bases. Returns TRUE if the sequence is a palindrome.
function is_palindrome(randomer)
    return occursin(randomer,reverse_complement(randomer))
end

"""
    degeneracy(base::DNA[, uselog2=true])

Calculate the degeneracy of a single base. The degeneracy
is the number of bases in the code, so
    degeneracy(DNA_N) == 4
and
    degeneracy(DNA_A) == 1

If `uselog2`, the degeneracy is log2 transformed. By convention,
the degneracy of `DNA_Gap` is zero in both the log2-transformed
and the untransformed case.
"""
function degeneracy(base::DNA; uselog2=true)
    deg = 0.0
    if base == DNA_N
        deg = 4.0
    elseif base == DNA_A || base == DNA_C || base == DNA_G || base == DNA_T
        deg = 1.0
    elseif base == DNA_B || base == DNA_D || base == DNA_H || base == DNA_V
        deg = 3.0
    elseif (base == DNA_K || base == DNA_Y || base == DNA_S || base == DNA_W
            || base == DNA_R || base == DNA_M)
        deg = 2.0
    else # DNA_Gap
        deg = 0.0
    end

    if uselog2
        if base == DNA_Gap
            deg = 2^-100  #log2(0) is negative infinity, so we use something close to zero
        end

        return log2(deg)

    else
        return deg
    end
end

"""
    degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)

Calculate the degeneracy of a DNA sequence.
"""
function degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)
    if uselog2
        return sum(degeneracy.(seq, uselog2=true))
    else
        return prod(degeneracy.(seq, uselog2=false))
    end
end

#Calculate the degeneracy of an entire oligo pool
function pooled_degeneracy(pool)
    deg=0;
    for randomer in pool
        deg +=degeneracy(randomer;uselog2=false);
    end
    return deg
end

#Determine if two oligos contain at least one common sequence. Used for assessing ambigous oligos.
#Extends the iscompatible function in the BioSequences package from bases to sequences.
#The two sequences must be the same length, otherwise the function returns false
function iscompatible(seq_A::LongSequence{DNAAlphabet{4}},seq_B::LongSequence{DNAAlphabet{4}})
    if length(seq_A) == length(seq_B)
        return prod(iscompatible.(seq_A,seq_B))
    else # Sequences are different lenghts
        return false
    end
end

#Decompress an ambiguous oligo into all of its possible sequences. Input is a single oligo.
#
function decompress_oligo(compressed_oligo::LongSequence{DNAAlphabet{4}})
    basedict=Dict(DNA_A => dna"A", DNA_T => dna"T", DNA_C => dna"C", DNA_G => dna"G",DNA_R =>dna"AG",DNA_Y=>dna"CT",DNA_S=>dna"GC",
    DNA_W=>dna"AT",DNA_K=>dna"GT",DNA_M=>dna"AC",DNA_B=>dna"CGT",DNA_D=>dna"AGT",DNA_H=>dna"ACT",DNA_V=>dna"ACG",DNA_N=>dna"ACGT",DNA_Gap=>dna"-")
    nseqs=trunc(Int128,degeneracy(compressed_oligo;uselog2=false))
    randomerlen=length(compressed_oligo)
    blockarray=zeros(Int128,randomerlen+1)
    blockarray[1]=nseqs
    degsarray=zeros(Int128,randomerlen)
    repsarray=zeros(Int128,randomerlen+1)
    repsarray[1]=1
    uncompressed_oligos=[dna"-"^randomerlen for _ in 1:nseqs]
    for i = 1:randomerlen
        nucleotides=basedict[compressed_oligo[i]]
        degsarray[i]=length(nucleotides)
        blockarray[i+1]=trunc(Int128,blockarray[i]/degsarray[i]);
        bases=dna"-"^blockarray[i];
        repsarray[i+1]=trunc(Int128,prod(degsarray[1:i]));
            for  k in 1:degsarray[i]
                for l in 1:blockarray[i+1]
                    bases[l+blockarray[i+1]*(k-1)]=nucleotides[k]
                end
            end

        bases=bases^repsarray[i];
        for m in 1:nseqs
            uncompressed_oligos[m][i]=bases[m];
        end
    end
    return uncompressed_oligos
end

function decompress_pool(compressed_pool)
    decompressed_pool=decompress_oligo(compressed_pool[1])
    for i = 2:length(compressed_pool)
        subset=decompress_oligo(compressed_pool[i])
        for i in subset
            push!(decompressed_pool,i)
        end
    end
    return decompressed_pool
end

#Generate a random pool without replacement. Specify the randomer length, the size of the pool, and the available nucleotides.
#It can create ambiguous sequences, but it is possible that unique sequences will appear multiple times.
function generate_random_pool(;randomerlen=6,randpoolsize=100,nucleotides=dna"ACGT")
    randpool=[]
    randomer=dna"-"^randomerlen
    for i = 1:randomerlen
        bases=shuffle(nucleotides)
        randomer[i]=bases[1];
    end
    push!(randpool,randomer)
    while length(randpool)<randpoolsize
        randomer=dna"-"^randomerlen
        for i = 1:randomerlen
            bases=shuffle(nucleotides)
            randomer[i]=bases[1];
        end
        check=true
        for j = 1:length(randpool)
            if randomer==randpool[j]
                check=false
                break
            end
        end
        if check
            push!(randpool,randomer)
        end
        if length(randpool)==length(nucleotides)^randomerlen
            break
        end
    end

    return randpool
end

# Given a pool of sequences all of length n, find all of the unambiguous sequences of length n that don't appear in that pool.
function valid_candidates(sites)
    lengths=length.(sites)
    if all(lengths[1].==lengths)
        n=length(sites[1])
        nsites=length(sites)
        candidates=generate_random_pool(;randomerlen=n,randpoolsize=4^n);
        for i=1:nsites
            candidates=candidates[findall(j -> !iscompatible(j,sites[i]),candidates)]
        end
        return candidates
    else
        println("provide an array of oligos of the same length")
        return nothing
    end
end

# Given a list of valid candidates, determine the action space of those candidates. For each base position, return a list of compatible nucleotides.
# If for example, the first postion does not contain an A or a T in any of the valid candidates, then the action space for that base position is dna"GCS".
#returns bases, an array of dna sequences containing all of the valid actions for the agent to take at each position.
function define_action_space(candidates)
    if length(candidates)==0
        bases=[dna"" for i in 1:6]
        return bases
    else
        randomerlen=length(candidates[1]);
        ncands=length(candidates)
        bases=[];
        for i=1:randomerlen
            states=dna"-"^ncands;
            nucleotides=dna"----" # A list of nucleotides not present at a given position in the oligos, called "states" here
            #println("TEST: $nucleotides")
            for j=1:ncands
                states[j]=candidates[j][i];
            end
            if !occursin(dna"A",states)
                nucleotides[1]=DNA_A
            else
                nucleotides[1]=DNA_Gap
            end
            if !occursin(dna"T",states)
                nucleotides[2]=DNA_T
            else
                nucleotides[2]=DNA_Gap
            end
            if !occursin(dna"C",states)
                nucleotides[3]=DNA_C
            else
            nucleotides[3]=DNA_Gap
            end
            if !occursin(dna"G",states)
                nucleotides[4]=DNA_G
            else
                nucleotides[4]=DNA_Gap
            end
            #println("TEST: $states")
            push!(bases,omit_codes(nucleotides))
            #println("Test: $nucleotides")
        end
    end
    return bases
end

#given a list of base codes, return all of the base codes that do not contain any of those codes
function omit_codes(bases::LongSequence{DNAAlphabet{4}})
    testcodes=dna"AGCTMRWSYKVHDBN"
    testcodes_copy=deepcopy(testcodes)
    indexarray=[]
    for i= 1:length(testcodes)
        count=0;
        for nucleotide in bases
            if iscompatible(nucleotide,testcodes[i])
                count+=1;
            end
        end
        if count >0
            push!(indexarray,i)
        end
    end
    check=sort(indexarray,rev=true)
    for index in check
        deleteat!(testcodes_copy,index)
    end

    return testcodes_copy
end

####################################################################
# Oligo Simulation Functions
####################################################################

"""
    isvalid(seq, sites)

Check if any of the sites appear in the sequence. Also check for palindromes by calling is palindrom

"""
function isvalid(seq, sites)
    if is_palindrome(seq)
        return false
    end
    for site in sites
        if occursin(site, seq)
            return false
        end
    end
    return true
end

"""
    find_valid_randomer(prefix, bases, sites)

Randomly select bases to add the prefix, ensuring no site
appears in the randomer.
"""
function find_valid_randomer(prefix, bases, sites)
    n = length(bases)
    pre_n = length(prefix)
    randomer = prefix * (dna"-" ^ n)
    for i = 1:n
        candidates=shuffle(bases[i])
        for cand in candidates
            randomer[pre_n+i] = cand
            if isvalid(randomer[1:pre_n+i], sites)
                break
            else
                randomer[pre_n+i] = DNA_Gap
            end
        end
        if randomer[pre_n+i] == DNA_Gap
            break
        end
    end
    return randomer
end

############################################################
#The follwoing functions are used as part of the various reward systems. They are called by each reward function when applicable
###########################################################

function measure_hitcounts(randomer::LongSequence{DNAAlphabet{4}},mRNAs)
    mRNAnum=length(mRNAs);
    hitcounts=zeros(mRNAnum)
    for i=1:mRNAnum
        counter=0;
        indexcounter=0
        index=approxsearchindex(mRNAs[i][1:end],randomer,0);
        indexcounter+=index;
        if index > 0
            counter +=1
        end
        while index >0
            index=approxsearchindex(mRNAs[i][indexcounter+1:end],randomer,0);
            indexcounter+=index
            if index==0
                break
            else
                counter+=1;
            end
        end
        hitcounts[i]=counter
    end
    return hitcounts
end

function mRNA_uniform_coverage_function_factory(randomer::LongSequence{DNAAlphabet{4}},mRNAs,hitcounts;uniformity_weight=1,total_hits_weight=1,kwargs...)
    hitcounts+=measure_hitcounts(randomer,mRNAs)
    function mRNA_uniform_coverage(randomer,mRNAs)
        mRNAnum=length(mRNAs);
        DULQ_score=0;
        new_hitcounts=measure_hitcounts(randomer,mRNAs)
        scored_hits=hitcounts+new_hitcounts;
        if sum(scored_hits)==0
            final_score=0
        else
            sorted_hits=sort(scored_hits);
            Davg=mean(sorted_hits)
            lq=round(Integer,mRNAnum/4);
            Dlq=mean(sorted_hits[1:lq])
            DULQ_score=Dlq/Davg;
            filtered_hits=filter(x -> x > 0, scored_hits)
            final_score=length(filtered_hits)+total_hits_weight*sum(log10.(filtered_hits))+uniformity_weight*mRNAnum*DULQ_score
        end
        return final_score
    end
    function mRNA_uniform_coverage_analysis(randomer,mRNAs)
        mRNAnum=length(mRNAs);
        DULQ_score=0;
        new_hitcounts=measure_hitcounts(randomer,mRNAs)
        scored_hits=hitcounts+new_hitcounts;
        if sum(scored_hits)==0
            final_score=0
        else
            sorted_hits=sort(scored_hits);
            Davg=mean(sorted_hits)
            lq=round(Integer,mRNAnum/4);
            Dlq=mean(sorted_hits[1:lq])
            DULQ_score=Dlq/Davg;
            filtered_hits=filter(x -> x > 0, scored_hits)
            final_score=length(filtered_hits)+total_hits_weight*sum(log10.(filtered_hits))+uniformity_weight*mRNAnum*DULQ_score
        end
        return final_score,hitcounts,DULQ_score, total_hits_weight, uniformity_weight
    end
    return mRNA_uniform_coverage,mRNA_uniform_coverage_analysis,hitcounts
end

function measure_hit_positions(randomer::LongSequence{DNAAlphabet{4}},mRNAs)
    mRNAnum=length(mRNAs)
    hit_positions=[[] for _ =1:mRNAnum]
    for i=1:mRNAnum
        indexcounter=0
        index=approxsearchindex(mRNAs[i][1:end],randomer,0);
        indexcounter+=index;
        if index > 0
            push!(hit_positions[i],indexcounter)
        end
        while index >0
            index=approxsearchindex(mRNAs[i][indexcounter+1:end],randomer,0);
            indexcounter+=index
            if index==0
                break
            else
                push!(hit_positions[i],indexcounter)
            end
        end
    end
    return hit_positions
end

function calculate_gap_distances(hit_indicies,genelen,randomerlen)
    dists=[]
    if length(hit_indicies)>0
        hit_indicies=sort(hit_indicies)
        push!(dists,hit_indicies[1])
        if length(hit_indicies)>1
            for i=2:length(hit_indicies)
                push!(dists,(hit_indicies[i]-hit_indicies[i-1]))
            end
        end
        push!(dists,(genelen-hit_indicies[end]))
    end
    return dists
end

#hit_positions= [[] for _ in 1:mRNAnum]
function mRNA_within_gene_uniformity_function_factory(randomer::LongSequence{DNAAlphabet{4}},mRNAs,hit_positions;uniformity_weight=1,total_hits_weight=1,kwargs...)
    hit_positions=vcat.(hit_positions,measure_hit_positions(randomer,mRNAs));
    function mRNA_within_gene_uniformity(randomer,mRNAs)
        mRNAnum=length(mRNAs)
        gene_lengths=length.(mRNAs)
        new_hit_positions=measure_hit_positions(randomer,mRNAs)
        scored_hit_positions=vcat.(hit_positions,new_hit_positions)
        randomerlen=length(randomer)
        if sum(length.(scored_hit_positions))==0
            final_score=0
        else
            gap_distances=calculate_gap_distances.(scored_hit_positions,gene_lengths,randomerlen)
            uniformity_score=zeros(mRNAnum)# "gap_standard_deviations_normalized" normalized by gene length
            for i =1:mRNAnum
                if length(gap_distances[i])==0
                    uniformity_score[i]=0
                else
                    sorted_gap_distances=sort(gap_distances);
                    Davg=mean(sorted_gap_distances)
                    lq=round(Integer,length(gap_distances)/4);
                    Dlq=mean(sorted_gap_distances[1:lq])
                    uniformity_score[i]=Dlq/Davg;
                end
            end
            filtered_hits=filter(x -> length(x)>0,scored_hit_positions)
            num_hits=length.(filtered_hits)
            final_score=length(filtered_hits)+ total_hits_weight*sum(log10.(num_hits))+ uniformity_weight*sum(uniformity_score)
        end
        return final_score
    end
    function mRNA_within_gene_uniformity_analysis(randomer,mRNAs)
        mRNAnum=length(mRNAs)
        gene_lengths=length.(mRNAs)
        new_hit_positions=measure_hit_positions(randomer,mRNAs)
        scored_hit_positions=vcat.(hit_positions,new_hit_positions)
        randomerlen=length(randomer)
        uniformity_score=zeros(mRNAnum)
        if sum(length.(scored_hit_positions))==0
            final_score=0
        else
            gap_distances=calculate_gap_distances.(scored_hit_positions,gene_lengths,randomerlen)
            # "gap_standard_deviations_normalized" normalized by gene length
            for i =1:mRNAnum
                if length(gap_distances[i])==0
                    uniformity_score[i]=0
                else
                    sorted_gap_distances=sort(gap_distances);
                    Davg=mean(sorted_gap+distances)
                    lq=round(Integer,length(gap_distances)/4);
                    Dlq=mean(sorted_gap_distances[1:lq])
                    uniformity_score[i]=Dlq/Davg;
                end
            end
            filtered_hits=filter(x -> length(x)>0,scored_hit_positions)
            num_hits=length.(filtered_hits)
            final_score=length(filtered_hits)+ total_hits_weight*sum(log10.(num_hits))+ uniformity_weight*sum(uniformity_score)
        end
        return final_score,hit_positions,mean(uniformity_score),total_hits_weight,uniformity_weight
    end
    return mRNA_within_gene_uniformity,mRNA_within_gene_uniformity_analysis, hit_positions
end

function calculate_DULQ_score(array)
    if length(array)==0
        DULQ_score=0
    else
        sorted_array=sort(array)
        grand_mean=mean(sorted_array)
        lq_idx=max(round(Integer,length(sorted_array)/4),1);
        lq_mean=mean(sorted_array[1:lq_idx])
        DULQ_score=lq_mean/grand_mean;
    end
    return DULQ_score
end

function mRNA_all_rewards_function_factory(randomer::LongSequence{DNAAlphabet{4}},mRNAs,hitcounts,hit_positions;genes_hit_weight=1,intrauniformity_weight=0,total_hits_weight=0,interuniformity_weight=0,kwargs...)
    if intrauniformity_weight==0  # we will only track the hit positions if the intrauniformity score matters. This will save unnecessary computation
        hitcounts+=measure_hitcounts(randomer,mRNAs)

        function mRNA_all_rewards_A(randomer,mRNAs)
            mRNAnum=length(mRNAs);
            interuniformity_score=0;
            new_hitcounts=measure_hitcounts(randomer,mRNAs)
            scored_hits=hitcounts+new_hitcounts;
            if sum(scored_hits)==0
                final_score=0
            else
                filtered_hits=filter(x -> x > 0, scored_hits)
                interuniformity_score=calculate_DULQ_score(filtered_hits)
                final_score=genes_hit_weight*length(filtered_hits)+total_hits_weight*sum(filtered_hits)+interuniformity_weight*mRNAnum*interuniformity_score
            end
            return final_score
        end
        function mRNA_all_rewards_analysis_A(randomer,mRNAs)
            mRNAnum=length(mRNAs);
            interuniformity_score=0;
            new_hitcounts=measure_hitcounts(randomer,mRNAs)
            scored_hits=hitcounts+new_hitcounts;
            genes_hit=0
            total_hits=0
            interuniformity_score=0
            intrauniformity_score=0
            if sum(scored_hits)==0
                final_score=0
                genes_hit=0
                total_hits=0
                interuniformity_score=0
                intrauniformity_score=0;
            else
                filtered_hits=filter(x -> x > 0, scored_hits)
                interuniformity_score=calculate_DULQ_score(filtered_hits)
                genes_hit=length(filtered_hits)
                total_hits=sum(filtered_hits)
                intrauniformity_score= -1
                final_score=genes_hit_weight*length(filtered_hits)+total_hits_weight*sum(filtered_hits)+interuniformity_weight*mRNAnum*interuniformity_score
            end
            return final_score,hitcounts,hit_positions,genes_hit,total_hits,interuniformity_score,mean(intrauniformity_score),genes_hit_weight,total_hits_weight,interuniformity_weight,intrauniformity_weight
        end
        mRNA_all_rewards=mRNA_all_rewards_A
        mRNA_all_rewards_analysis=mRNA_all_rewards_analysis_A

    else ## intrauniformity_weight >0, keep track of hit positions
        hit_positions=vcat.(hit_positions,measure_hit_positions(randomer,mRNAs));
        hitcounts=length.(hit_positions)
        function mRNA_all_rewards_B(randomer,mRNAs)
            mRNAnum=length(mRNAs)
            gene_lengths=length.(mRNAs)
            new_hit_positions=measure_hit_positions(randomer,mRNAs)
            scored_hit_positions=vcat.(hit_positions,new_hit_positions)
            randomerlen=length(randomer)
            if sum(length.(scored_hit_positions))==0
                final_score=0
            else
                gap_distances=calculate_gap_distances.(scored_hit_positions,gene_lengths,randomerlen)
                gap_distances=filter(x -> length(x)>0,gap_distances)
                intrauniformity_scores=calculate_DULQ_score.(gap_distances)
                filtered_hit_positions=filter(x -> length(x)>0,scored_hit_positions)
                filtered_hitcounts=(length.(filtered_hit_positions));
                interuniformity_score=calculate_DULQ_score(filtered_hitcounts)
                filtered_hits=length.(filtered_hit_positions)
                final_score=genes_hit_weight*length(filtered_hits)+ total_hits_weight*sum(filtered_hits)+ interuniformity_weight*mRNAnum*interuniformity_score+intrauniformity_weight*mRNAnum*mean(intrauniformity_scores)
            end
            return final_score
        end
        function mRNA_all_rewards_analysis_B(randomer,mRNAs)
            mRNAnum=length(mRNAs)
            gene_lengths=length.(mRNAs)
            new_hit_positions=measure_hit_positions(randomer,mRNAs)
            scored_hit_positions=vcat.(hit_positions,new_hit_positions)
            randomerlen=length(randomer)
            genes_hit=0
            total_hits=0
            interuniformity_score=0
            intrauniformity_scores=0
            if sum(length.(scored_hit_positions))==0
                final_score=0
                genes_hit=0
                total_hits=0
                interuniformity_score=0
                intrauniformity_scores=0
            else
                gap_distances=calculate_gap_distances.(scored_hit_positions,gene_lengths,randomerlen)
                gap_distances=filter(x -> length(x)>0,gap_distances)
                intrauniformity_scores=calculate_DULQ_score.(gap_distances)
                filtered_hit_positions=filter(x -> length(x)>0,scored_hit_positions)
                filtered_hitcounts=(length.(filtered_hit_positions));
                interuniformity_score=calculate_DULQ_score(filtered_hitcounts)
                filtered_hits=length.(filtered_hit_positions)
                genes_hit=length(filtered_hits)
                total_hits=sum(filtered_hits)
                final_score=genes_hit_weight*length(filtered_hits)+ total_hits_weight*sum(filtered_hits)+ interuniformity_weight*mRNAnum*interuniformity_score+intrauniformity_weight*mRNAnum*mean(intrauniformity_scores)
            end
            return final_score,hitcounts,hit_positions,genes_hit,total_hits,interuniformity_score,mean(intrauniformity_scores),genes_hit_weight,total_hits_weight,interuniformity_weight,intrauniformity_weight
        end
        mRNA_all_rewards=mRNA_all_rewards_B
        mRNA_all_rewards_analysis=mRNA_all_rewards_analysis_B
    end
    return mRNA_all_rewards,mRNA_all_rewards_analysis,hitcounts,hit_positions
end
#################################################
# Simulation function. Call Reward functions and simulate a valid randomer
##################################################
function simulate(prefix,bases,sites,mRNAs,reward_function; nsims=1000, horizon=length(bases), kwargs...)
    rewards=zeros(nsims)
    true_horizon = min(horizon, length(bases))
    n=length(prefix)+true_horizon;
    sim_randomer_list=[dna"-"^n for _ in 1:nsims];

    Threads.@threads for i = 1:nsims   #This is where multithreading is incorporated
        sim_randomer_list[i]=find_valid_randomer(prefix,bases[1:true_horizon],sites);
        rewards[i]= reward_function(sim_randomer_list[i],mRNAs)
    end
    max_sim_score,max_index=findmax(rewards);
    max_sim_randomer=sim_randomer_list[max_index];
    return mean(rewards),max_sim_randomer,max_sim_score

end

########################################
# Rollout Function
#######################################
function oligo_rollout(bases,sites, mRNA; reward_function=mRNA_uniform_coverage, kwargs...)
    n=length(bases)
    randomer= dna"-" ^ n
    # Find the longest blocking site and use this as the horizon. (they should all be equal)
    max_len = map(length, sites) |> maximum
    candidates=valid_candidates(sites);
    best_sim_randomer=dna"-"^n
    best_sim_score=0;
    for i= 1:n
        best_base= DNA_Gap
        best_score=0;
        start = max(i - max_len, 1) # only need to look back max_len when checking the oligos
        bases=define_action_space(candidates[findall(j -> iscompatible(j[start:i-1],randomer[start:i-1]),candidates)])
        if length(bases[i])==1
            best_base=bases[i][1];
        else
            for base in bases[i]
                if ~isvalid(randomer[start:i-1] * base, sites)
                    # adding this bases creates a restriction site; skip it
                    continue
                end
                score,max_sim_randomer,max_sim_score= simulate(randomer[start:i-1] * base, bases[i+1:end], sites, mRNA,reward_function;kwargs...)
                # Update use this base if:
                #   1. the mean randomer mRNA coverage is higher than the previous best, OR
                #   2. the randomer mRNA coverage is tied with the previous best, but the new base has
                #      higher individual degeneracy
                if score > best_score || (score == best_score && degeneracy(base) > degeneracy(best_base))
                    best_score = score
                    best_base = base
                end
                ## Update the overall best simulation result if
                #   1. the best globally simulated randomer mRNA score is better than the previous best
                #   2. the best globally simulated randomer mRNA score is tied but the randomer has higher degeneracy
                if max_sim_score >best_sim_score || (max_sim_score==best_sim_score && degeneracy(max_sim_randomer)>degeneracy(best_sim_randomer))
                    best_sim_score=max_sim_score;
                    best_sim_randomer=max_sim_randomer
                end
            end
        end
        if best_base == DNA_Gap
            # The sequence has terminated; there are no bases that can be added without making
            # a restriction site.
            break
        else
            randomer[i] = best_base

        end
        #println("randomer: $randomer")
    end
    final_randomer_score=reward_function(randomer,mRNA)
    if best_sim_score >final_randomer_score
        randomer=best_sim_randomer
        println("Triggered Best Simulated Randomer")
    end
    return randomer
end

####################################
#Functions for Reading Genome Data CSV Files Imported from R
######################################

## Read in the list of sequences you want to exclude in the designed primers.
# Function finds every unique hexamer in the list of sequences and stores it as an array
#. For example, use list of rRNA and tRNA sequences. This funciton takes care of orientation. We end up using the reverse complement since we are finding primers that won't hybridize
## Example
#blocking_sites=read_blocking_sites(CSV.read("./SMU_UA159_rRNA_tRNA.csv"))
function read_blocking_sites(data; randomerlen=6)
    blocking_sites=[];
    for i =1:nrow(data)
        seq=reverse_complement(LongDNASeq(string(data[i,:Sequence])));
        for j=1:length(seq)-(randomerlen-1)
            if isvalid(seq[j:j+(randomerlen-1)],blocking_sites)
                push!(blocking_sites,seq[j:j+(randomerlen-1)])
            end
        end
    end
    return blocking_sites
end
# RNAs=read_rRNA_tRNA_list(CSV.read("./SMU_UA159_rRNA_tRNA.csv"))
#This function Does not hash up the rRNAs and tRNAs like read_blocking_sites does. That is the only difference
function read_rRNA_tRNA_list(data)
    RNAs=[]
    for i=1:nrow(data)
        seq=reverse_complement(LongDNASeq(string(data[i,:Sequence])));
        if isvalid(seq,RNAs)
            push!(RNAs,seq)
        end
    end
    return RNAs
end
## Read in the list of sequences you want the designed primers to hybridize to
# Function takes care of the orientation, we want to end up using the reverse complmement.
## Example
# mRNAs=read_mRNA_list(CSV.read("./SMU_UA159_Genes.csv"))
function read_mRNA_list(data)
    mRNAs=[]
    for i=1:nrow(data)
        seq=reverse_complement(LongDNASeq(string(data[i,:mRNA_seq])));
        if isvalid(seq,mRNAs)
            push!(mRNAs,seq)
        end
    end
    return mRNAs
end

#Finds the number of bases as well as the GC content of the given file
# Example
# bases,GC_content=calculate_transcriptome_stats(read_mRNA_list(CSV.read("./SMU_UA159_Genes.csv")))
function calculate_transcriptome_stats(sequences)
    bases=0;  #number of bases in the file
    GC_count=0 # Count of G's and C's in file
    for i=1:length(sequences)
        bases += length(sequences[i]);
        comp=composition(sequences[i]);
        GC_count += comp[DNA_G] +comp[DNA_C];
    end
    GC_content= GC_count/bases*100;
    return bases,GC_content
end

#################################################
#### Extra Analysis Tools
#################################################

## Finds the edit distance distribution for a given designed primer to all of the blocking sites
function base_error_dist(randomer,sites)
    error_dist=[]
    for i in 0:length(randomer)
        count=0
        found_sites=[]
        for j in 1:length(sites)
            valid=approxsearchindex(randomer,sites[j],i)
            if valid == 1
                count += valid
                push!(found_sites,j)
            end
        end
        deleteat!(sites,found_sites)
        push!(error_dist,count)
    end

    return error_dist
end

function ID_genes_missing(data,genes_missing)
    missing_ID=[];
    for i in genes_missing
        push!(missing_ID,data[i,:New_ID])
    end
    return missing_ID
end

## Calculates the distribution number of randomers that hit each gene at least once. Not very useful for most cases
function coverage_distribution(randomers, mRNA)
    coverage_dist=zeros(length(randomers)+1)
     for seq in mRNA
         count=0;
         for randomer in randomers
             if occursin(randomer, seq)
                 count+=1;
             end
         end
         coverage_dist[count+1]+=1;

     end
     return coverage_dist
end

# Calcualtes the number of times each gene is hit by all randomers in the pool. Useful for comparing the "between gene uniformity".
function gene_coverage_counts(randomers,mRNAs)
    gene_hits=zeros(length(mRNAs))
    randomerlen=length(randomers[1]);
    for i in 1:length(mRNAs)
        counter=0;
        for j in 1:length(mRNAs[i])-randomerlen
            for k in 1:length(randomers)
                if occursin(randomers[k],mRNAs[i][j:j+randomerlen-1])
                    counter+=1
                end
            end
        end
        gene_hits[i]=counter
    end
    data=DataFrame(counts=gene_hits)
    return data

end

 ########################################################################
 # MAIN RUNNING FUNCTION
 ########################################################################

 #Quick Test#
 # data=run_NSR_RL(pool_size=3,nsims=50)
 #CSV.write("NSR_RL_S_mutans_4_12_21",data)

 #Input the Species name, the length of primer, and the size of the pool.
function run_NSR_RL(;species="S_mutans",randomerlen=6,pool_size=25,all_bases = dna"AGCTMRWSYKVHDBN",kwargs...)

    #Intialization
    randomers=[]
    mRNAfile=string("./Genome Data/",species,"_Genes.csv")
    rRNAfile=string("./Genome Data/",species,"_rRNA_tRNA.csv")
    mRNA_DF=CSV.read(mRNAfile,DataFrame);
    rRNA_DF=CSV.read(rRNAfile,DataFrame);
    mRNAs=read_mRNA_list(mRNA_DF);
    rRNAs=read_rRNA_tRNA_list(rRNA_DF);
    numgenes=length(mRNAs)
    blocking_sites=read_blocking_sites(rRNA_DF,randomerlen=randomerlen);
    mRNAbp,mRNA_GC=calculate_transcriptome_stats(mRNAs);
    rRNAbp,rRNA_GC=calculate_transcriptome_stats(rRNAs);
    transcriptome_size=mRNAbp+rRNAbp;
    all_ns = [all_bases for _ in 1:randomerlen];
    times=[]
    num_blocked_sites=zeros(pool_size)
    println("Initialization Successful")
    ##Initialize the reward-associated data structure and reward function here. ex hitcounts=zeros(numgenes), hit_positions=[[] for _ in 1:numgenes]
    hitcounts=zeros(numgenes)
    hit_positions=[[] for _ in 1:numgenes]
    #Call the appropriate reward function factory
    reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(dna"-"^randomerlen,mRNAs,hitcounts,hit_positions;kwargs...)
    ##
    for i = 1:pool_size
        num_blocked_sites[i]=length(blocking_sites)
        randomerstats= @timed oligo_rollout(all_ns,blocking_sites,mRNAs; reward_function=reward_function,kwargs...); #create new randomer and add to pool
        push!(randomers,randomerstats.value)
        push!(times,randomerstats.time)
        push!(blocking_sites,LongDNASeq(randomers[i]));
        push!(blocking_sites,LongDNASeq(reverse_complement(randomers[i]))); ## avoid adding randomers that are reverse complements of each other
        blocking_sites=unique(decompress_pool(blocking_sites));
        reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(randomers[i],mRNAs,hitcounts,hit_positions; kwargs...) #officially add hit positions from the new randomer to the reward function
        #println("Iteration $i")
    end
    final_score,hitcounts,hit_positions,genes_hit,total_hits,interuniformity_score,intrauniformity_score,genes_hit_weight,total_hits_weight,interuniformity_weight,intrauniformity_weight =analysis_function(dna"-"^randomerlen,mRNAs)
    genes_missing=findall(x -> x==0,hitcounts)
    missing_IDs=join(ID_genes_missing(mRNA_DF,genes_missing),",")
    totaldeg=pooled_degeneracy(randomers)
    total_time=sum(times)
    string_time=join(times,",")
    string_num_blocked_sites=join(num_blocked_sites,",")
    string_randomers=join(dnaarray2stringarray(randomers),",")
    data=DataFrame(poolsize=pool_size, primer_length=randomerlen, genes_hit_weight=genes_hit_weight, total_hits_weight=total_hits_weight,interuniformity_weight=interuniformity_weight,intrauniformity_weight=intrauniformity_weight,species=species,num_genes=numgenes,
        transcriptome_size=transcriptome_size,mRNA_GC=mRNA_GC,blocked_GC=rRNA_GC,randomers=string_randomers,degeneracy=totaldeg,num_genes_hit=genes_hit,
        genes_missing=missing_IDs,total_hits=total_hits,interuniformity_score=interuniformity_score,intrauniformity_score=intrauniformity_score,individual_times=string_time,time=total_time,num_blocked_sites=string_num_blocked_sites,final_score=final_score)
    return data
end

function benchmark_NSR_RL(;species=["S_mutans"],genes_hit_weights=[1],total_hits_weights=[1e-4],interuniformity_weights=[1],intrauniformity_weights=[1],kwargs...)
    DF1=DataFrame()
    a=length(genes_hit_weights)
    b=length(total_hits_weights)
    c=length(interuniformity_weights)
    d=length(intrauniformity_weights)
    e=length(species)
    if a != b & a != c & a != d & a !=e
        error("Provide weights for equal numbers of experiments.")
    end
    nexps=a;
    for i = 1:nexps
        data=run_NSR_RL(;species=species[i],genes_hit_weight=genes_hit_weights[i],total_hits_weight=total_hits_weights[i],interuniformity_weight=interuniformity_weights[i],intrauniformity_weight=intrauniformity_weights[i],kwargs...)
        DF1=vcat(DF1,data);
    end
    return DF1
end

##Code for analyzing previously made primer pools. stores the same data that run_NSR_RL stores
function analyze_NSR_Pool(randomers;species="S_mutans",mRNAfile="./SMU_UA159_Genes.csv",rRNAfile="./SMU_UA159_rRNA_tRNA.csv",kwargs...)
    mRNAfile=string("./Genome Data/",species,"_Genes.csv")
    rRNAfile=string("./Genome Data/",species,"_rRNA_tRNA.csv")
    mRNA_DF=CSV.read(mRNAfile,DataFrame);
    rRNA_DF=CSV.read(rRNAfile,DataFrame);
    mRNAs=read_mRNA_list(mRNA_DF);
    rRNAs=read_rRNA_tRNA_list(rRNA_DF);
    numgenes=length(mRNAs)
    blocking_sites=read_blocking_sites(rRNA_DF,randomerlen=randomerlen);
    randomerlen=length(randomers[1]);
    mRNAbp,mRNA_GC=calculate_transcriptome_stats(mRNAs);
    rRNAbp,rRNA_GC=calculate_transcriptome_stats(rRNAs);
    transcriptome_size=mRNAbp+rRNAbp;
    hitcounts=zeros(numgenes)
    hit_positions=[[] for _ in 1:numgenes]
    reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(dna"-"^randomerlen,mRNAs,hitcounts,hit_positions;genes_hit_weight=1,total_hits_weight=1,interuniformity_weight=1,intrauniformity_weight=1)
    for i=1:length(randomers)
        reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(randomers[i],mRNAs,hitcounts,hit_positions;genes_hit_weight=1,total_hits_weight=1,interuniformity_weight=1,intrauniformity_weight=1); #officially add hit positions from the new randomer
    end
    final_score,hitcounts,hit_positions,genes_hit,total_hits,interuniformity_score,intrauniformity_score,genes_hit_weight,total_hits_weight,interuniformity_weight,intrauniformity_weight =analysis_function(dna"-"^randomerlen,mRNAs)
    genes_missing=findall(x->x==0,hitcounts);
    missing_IDs=join(ID_genes_missing(mRNAdata,genes_missing),",")
    totaldeg=pooled_degeneracy(randomers)
    total_time=0
    string_time="N/A"
    string_randomers=join(dnaarray2stringarray(randomers),",")
    genes_hit_weight=0
    total_hits_weight=0
    interuniformity_weight=0
    intrauniformity_weight=0
    pool_size=length(randomers);
    data=DataFrame(poolsize=pool_size, primer_length=randomerlen, genes_hit_weight=genes_hit_weight, total_hits_weight=total_hits_weight,interuniformity_weight=interuniformity_weight,intrauniformity_weight=intrauniformity_weight,species=species,num_genes=numgenes,
        transcriptome_size=transcriptome_size,mRNA_GC=mRNA_GC,blocked_GC=rRNA_GC,randomers=string_randomers,degeneracy=totaldeg,num_genes_hit=genes_hit,
        genes_missing=missing_IDs,total_hits=total_hits,interuniformity_score=interuniformity_score,intrauniformity_score=intrauniformity_score,individual_times=string_time,time=total_time)
    return data
end

function analyze_all_genomes_NSR_Pools(species_list)
    DF1=DataFrame()
    n=length(species_list)
    for i = 1:n
        randomers=stringarray2dnaarray(CSV.read(string("./Brute Force Pools V2 Compressed/",species_list[i],"_NSR_Brute_Force_Compressed.csv"),DataFrame)[:Sequence])
        mRNAfile=string("./Genome Data/",species_list[i],"_Genes.csv")
        rRNAfile=string("./Genome Data/",species_list[i],"_rRNA_tRNA.csv")
        data=analyze_NSR_Pool(randomers;species=species_list[i],mRNAfile=mRNAfile,rRNAfile=rRNAfile)
        DF1=vcat(DF1,data);
        pct_complete=i/n*100;
        println("$pct_complete % Complete")
        CSV.write("./Brute_Force_Pool_Compressed_Analysis_all_rewards.csv",DF1)
    end

end
function analyze_NSR_Pool_gene_coverage(species_name)
    randomers=stringarray2dnaarray(CSV.read(string("./Brute Force Pools V2/",species_name,"_NSR_Brute_Force_all_randomers.csv"),DataFrame)[:Sequence])
    mRNAfile=string("./Genome Data/",species_name,"_Genes.csv")
    data=gene_coverage_counts(randomers,mRNAfile);
    outputfile=string("./Brute Force Pools V2 Gene Coverage/",species_name,"_Gene_Coverage.csv")
    CSV.write(outputfile,data)
end

function analyze_cumulative_NSR_Pool(randomers,species_name,kwargs...)
    mRNAfile=string("./Genome Data/",species_name,"_Genes.csv")
    mRNAdata=CSV.read(mRNAfile,DataFrame);
    mRNAs=read_mRNA_list(mRNAdata);
    nrandomers=length(randomers);
    randomerlen=length(randomers[1]);
    cumulative_hits=zeros(nrandomers)
    cumulative_score=zeros(nrandomers)
    cumulative_interuniformity_score=zeros(nrandomers)
    cumulative_intrauniformity_score=zeros(nrandomers)
    randomer_num=1:nrandomers;
    str_randomers=dnaarray2stringarray(randomers);
    cumulative_genes_hit=zeros(nrandomers);
    hitcounts=zeros(length(mRNAs))
    hit_positions=[[] for _ in 1:length(mRNAs)]
    reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(dna"-"^randomerlen,mRNAs,hitcounts,hit_positions;genes_hit_weight=1,total_hits_weight=1,interuniformity_weight=1,intrauniformity_weight=1)
    for i=1:nrandomers
        reward_function,analysis_function,hitcounts,hit_positions=mRNA_all_rewards_function_factory(randomers[i],mRNAs,hitcounts,hit_positions;genes_hit_weight=1,total_hits_weight=1,interuniformity_weight=1,intrauniformity_weight=1); #officially add hit positions from the new randomer
        final_score,hitcounts,hit_positions,genes_hit,total_hits,interuniformity_score,intrauniformity_score,genes_hit_weight,total_hits_weight,interuniformity_weight,intrauniformity_weight=analysis_function(dna"-"^randomerlen,mRNAs)
        cumulative_score[i]=final_score
        cumulative_hits[i]=total_hits
        cumulative_interuniformity_score[i]=interuniformity_score
        cumulative_genes_hit[i]=genes_hit
        cumulative_intrauniformity_score[i]=intrauniformity_score
    end
    data=DataFrame(run=randomer_num,randomer=str_randomers,cumulative_genes_hit=cumulative_genes_hit,cumulative_hits=cumulative_hits,cumulative_interuniformity_score=cumulative_interuniformity_score,cumulative_intrauniformity_score=cumulative_intrauniformity_score,cumulative_score=cumulative_score)
    return data
end

function cumulative_runtime()
    data=CSV.read("./S_mutans_NSR_Refactor_NSR_RL_All_Rewards_BF_Comparison_6_21_21.csv",DataFrame)
    hitcountstime=parse.(Float64,(split(data[2,19],",")))
    hit_pos_time=parse.(Float64,(split(data[4,19],",")))
    run=1:length(hitcountstime)
    cum_hitcountstime=zeros(length(hitcountstime))
    cum_hit_pos_time=zeros(length(hit_pos_time))
    for i =1:length(hit_pos_time)
        cum_hitcountstime[i]=sum(hitcountstime[1:i])
        cum_hit_pos_time[i]=sum(hit_pos_time[1:i])
    end
    df=DataFrame(run=run,without_intra=cum_hitcountstime,with_intra=cum_hit_pos_time)
    CSV.write("NSR_pool_size_timing_6_23_21.csv",df)
    end

###################################
# Run this code to re-create the experiments for Figure 3.
###################################

# For Panels C-F
s_mutans_data=benchmark_NSR_RL(;pool_size=100,nsims=100,randomerlen=6,species=["S_mutans"],genes_hit_weights=[1],total_hits_weights=[1e-4],interuniformity_weights=[1],intrauniformity_weights=[1])
CSV.write("./S_mutans_benchmark.csv",s_mutans_data)
s_mutans_randomers=stringarray2dnaarray(split(s_mutans_data[:randomers][1],","))
s_mutans_cumulative_data=analyze_cumulative_NSR_Pool(s_mutans_randomers,"S_mutans")
CSV.write("./S_mutans_benchmark_cumulative.csv",s_mutans_cumulative_data)

# For Panel G
species_list=CSV.read("./Genome Data/Species_List.csv",DataFrame)[:Species]
weight_vector=repeat([1],length(species_list))
all_species_data=benchmark_NSR_RL(;pool_size=30,randomerlen=6,nsims=100,species=species_list,genes_hit_weights=weight_vector,total_hits_weights=weight_vector,interuniformity_weights=weight_vector,intrauniformity_weights=weights_vector)
CSV.write("./All_species_benchmark.csv",all_species_data)

# For Panel H
intrauniformity_timing_data=benchmark_NSR_RL(;pool_size=100,randomerlen=6,nsims=100,species=["S_mutans","S_mutans"],genes_hit_weights=[1 1],total_hits_weights=[1 1],interuniformity_weights=[1 1],intrauniformity_weights=[1 0])
CSV.write("./Intrauniformity_timing_benchmark.csv",intrauniformity_timing_data)
num_oligos=length(stringarray2dnaarray(split(intrauniformity_timing_data[:randomers][1],",")))
run=1:num_oligos
intra0=parse.(Float64,split(intrauniformity_timing_data[:individual_times][1],","))
intra1=parse.(Float64,split(intrauniformity_timing_data[:individual_times][2],","))
cumulative_intra0=zeros(num_oligos)
cumulative_intra1=zeros(num_oligos)
for i = 1:num_oligos
    cumulative_intra0[i]=sum(intra0[1:i])
    cumulative_intra1[i]=sum(intra1[1:i])
end
intrauniformity_timing_cumulative=DataFrame(run=run,with_intra=cumulative_intra1,without_intra=cumulative_intra0)
CSV.write("./Intrauniformity_timing_benchmark_cumulative.csv",intrauniformity_timing_cumulative)


text:
using Random

function count_elt(map_el::MapElites)
    return length(collect(skipmissing(map_el.solutions)))
end

function add_to_map(map_el::MapElites,coordinates,solution::Individual,performance::Array{Float64})
    if (!haskey(map_el.solutions.data,Tuple(coordinates)))
        map_el.solutions.data[Tuple(coordinates)] = solution
        map_el.performances.data[Tuple(coordinates)] = performance
    elseif (map_el.performances.data[Tuple(coordinates)] < performance)
        map_el.solutions.data[Tuple(coordinates)] = solution
        map_el.performances.data[Tuple(coordinates)] = performance
    end
end

function select_random(map_el::MapElites)
    available_ind = collect(skipmissing(map_el.solutions))
    nb_ind = length(available_ind)
    rand_idx = Random.rand(1:nb_ind)
    return available_ind[rand_idx]
end

function select_most_promising(map_el::MapElites)
    available_ind = collect(skipmissing(map_el.solutions))
    available_perf = collect(skipmissing(map_el.performances))
    idx_best = argmax(available_perf)
    return available_ind[idx_best]
end

ext: jl


text:
using NeuralNetworkEconomics

runweave= "weave" ∈ ARGS
runnotebook= "notebook" ∈ ARGS

src=ARGS[1]

if runweave
  println("weaving markdown for $src")
  using Weave
  weave(src,out_path="md",
        cache=:refresh, cache_path="weavecache",
        doctype="github", mod=Main,
        args=Dict("md" => true))
end

if runnotebook
  println("weaving notebook for $src")
  using Weave
  notebook(src, out_path=joinpath(pwd(),"build"), nbconvert_options="--allow-errors")
end

if (isfile("build/temp.md"))
  rm("build/temp.md")
end

max_stars_repo_name: schrimpf/GradientBoostingNotes

max_issues_repo_name: schrimpf/GradientBoostingNotes

max_forks_repo_name: schrimpf/GradientBoostingNotes


text:
import Base: *, /, +, -, ^, sin, cos, exp, log, abs, abs2, sign, tanh, sqrt
import Base: sum, dropdims, transpose, adjoint, minimum, maximum, getindex, setindex!, reshape
import LinearAlgebra: mul!
import Statistics: mean
import Espresso: ExGraph, ExNode, matchingex, rewrite
import Espresso
using Requires

include("fwddecl.jl")
include("utils.jl")
include("tape.jl")
include("tracked.jl")
include("ops.jl")
include("forward/scalar.jl")
include("forward/tensor.jl")
include("backward/scalar.jl")
include("backward/tensor.jl")
include("macro.jl")
include("compile.jl")
include("devices.jl")
include("grad.jl")
include("update.jl")
include("helpers.jl")

function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("cuda.jl")
end


text:
"""
`module  JAC.MultipoleMoment`
    ... a submodel of JAC that contains all methods for computing (electric) dipole D_z amplitudes
        between two bound-state levels.
"""
module MultipoleMoment

    using Printf, JAC, ..Basics, ..Defaults, ..ManyElectron, ..Radial, ..InteractionStrength

    """
    `MultipoleMoment.amplitude(mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                                   grid::Radial.Grid; display::Bool=false)`
        ... to compute the multipole-moment amplitude
            <alpha_f J_f || Q^(M) (omega, gauge) || alpha_i J_i>  for the given final and initial level. All the amplitudes are
            calculated by means of Johnson's (2007) multipole-transition amplitudes from which also all constrains w.r.t the
            allowed gauges are inherited. A value::ComplexF64 is returned.
    """
    function amplitude(mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                       grid::Radial.Grid; display::Bool=false)

        amplitude = MultipoleMoment.transitionAmplitude(mp, gauge, omega, finalLevel, initialLevel, grid; display=false)
        #
        if  display
            sa = @sprintf("%.5e", amplitude.re) * "  " * @sprintf("%.5e", amplitude.im)
            println("    < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || Q^($mp) ($omega a.u., $gauge) ||" *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "    Multipole-moment amplitude:     " *
                                  " < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || Q^($mp) ($omega a.u., $gauge) ||" *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

    """
    `MultipoleMoment.dipoleAmplitude(finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)`
         ... to compute the dipole amplitude   <(alpha_f J_f, kappa) J_i || D || alpha_i J_i>  for the given final
            and initial level. A value::ComplexF64 is returned.
    """
    function dipoleAmplitude(finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)

        if     finalLevel.parity == initialLevel.parity
            amplitude = 0.0 + 0.0 * im
        else
            nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
            printstyled("Compute dipole matrix of dimension $nf x $ni in the final- and initial-state bases " *
                        "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
            matrix = zeros(ComplexF64, nf, ni)
            #
            for  r = 1:nf
                for  s = 1:ni
                    ##x wa = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                    # Calculate the spin-angular coefficients
                    if  Defaults.saRatip()
                        waR = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                        wa  = waR
                    end
                    if  Defaults.saGG()
                        subshellList = initialLevel.basis.subshells
                        opa = SpinAngular.OneParticleOperator(1, plus, true)
                        waG = SpinAngular.computeCoefficients(opa, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s], subshellList)
                        wa  = waG
                    end
                    if  Defaults.saRatip() && Defaults.saGG() && true
                        if  length(waR) != 0     println("\n>> Angular coeffients from GRASP/MCT   = $waR ")    end
                        if  length(waG) != 0     println(  ">> Angular coeffients from SpinAngular = $waG ")    end
                    end
                    #
                    for  coeff in wa
                        ja = Basics.subshell_2j(finalLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(initialLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = JAC.InteractionStrength.dipole(finalLevel.basis.orbitals[coeff.a], initialLevel.basis.orbitals[coeff.b], grid)
                        matrix[r,s] = matrix[r,s] + coeff.T * tamp
                    end
                end
            end
            printstyled("done. \n", color=:light_green)
            amplitude::ComplexF64 = transpose(finalLevel.mc) * matrix * initialLevel.mc
        end
        #
        if  display
            sa = @sprintf("%.5e", amplitude.re) * "  " * @sprintf("%.5e", amplitude.im)
            println("    < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || D ||" *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "    Dipole amplitude:  < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || D ||" *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

    """
    `JAC.MultipoleMoment.transitionAmplitude(mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                                             grid::Radial.Grid; display::Bool=false)`
        ... to compute Johnson's (2007) multipole-transition amplitude  <alpha_f J_f || T^(M, absorption) (omega, gauge) || alpha_i J_i>
            for the given final and initial level. Allowed gauges are {Velocity, Length} for the electric-multipole transition
            amplitudes and Magnetic for the magnetic-multipole transition amplitudes, respectively. A value::ComplexF64 is
            returned.
    """
    function transitionAmplitude(mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                                 grid::Radial.Grid; display::Bool=false)
        nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
        printstyled("Compute multipole-moment transition matrix of dimension $nf x $ni in the final- and initial-state bases " *
                    "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
        matrix = zeros(ComplexF64, nf, ni)
        #
        for  r = 1:nf
            for  s = 1:ni
                wa = compute("angular coefficients: 1-p, Grasp92", 0, mp.L, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                for  coeff in wa
                    ##x ja = Basics.subshell_2j(finalLevel.basis.orbitals[coeff.a].subshell)
                    ##x jb = Basics.subshell_2j(initialLevel.basis.orbitals[coeff.b].subshell)
                    tamp  = JAC.InteractionStrength.multipoleTransition(mp, gauge, omega, finalLevel.basis.orbitals[coeff.a],
                                                                        initialLevel.basis.orbitals[coeff.b], grid)
                    matrix[r,s] = matrix[r,s] + coeff.T * tamp
                end
            end
        end
        printstyled("done. \n", color=:light_green)
        amplitude = transpose(finalLevel.mc) * matrix * initialLevel.mc
        #
        if  display
            sa = @sprintf("%.5e", amplitude.re) * "  " * @sprintf("%.5e", amplitude.im)
            println("    < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] ||" *
                    " T^($mp, absorption) ($omega a.u., $gauge) ||" *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "    Multipole-transition amplitude:      < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] ||" *
                                  " T^($mp, absorption) ($omega a.u., $gauge) ||" *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

end # module


text:
using DashBootstrapComponents

checklist = html_div([
    dbc_checklist(
        id = "checklist-selected-style",
        options = [
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
            Dict("label" => "Option 3", "value" => 3),
        ],
        label_checked_style = Dict("color" => "red"),
        input_checked_style = Dict(
            "backgroundColor" => "#fa7268",
            "borderColor" => "#ea6258",
        ),
    ),
    html_hr(),
    dbc_radioitems(
        id = "radio-selected-style",
        options = [
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
            Dict("label" => "Option 3", "value" => 3),
        ],
        labelCheckedClassName = "text-success",
        inputCheckedClassName = "border border-success bg-success",
    ),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
include("thermalgeneration_constructor.jl")
include("branch_constructor.jl")
include("renewablegeneration_constructor.jl")
include("load_constructor.jl")
include("storage_constructor.jl")


text:
module JuliaAPI

using Compat
using Compat.REPL

function eval_str(code::AbstractString;
                  scope::Module = Main,
                  filename::AbstractString = "string",
                  auto_jlwrap = true,
                  force_jlwrap = false)
    return include_string(scope, code, filename)
end

@static if VERSION < v"0.7-"
    getattr(obj, name) = getfield(obj, Symbol(name))
    function setattr(obj, name, value)
        setfield!(obj, Symbol(name), value)
        return nothing
    end
else
    getattr(obj, name) = getproperty(obj, Symbol(name))
    function setattr(obj, name, value)
        setproperty!(obj, Symbol(name), value)
        return nothing
    end
end

function setattr(mod::Module, name, value)
    Base.eval(mod, :($(Symbol(name)) = $value))
    return nothing
end

@static if VERSION < v"0.7-"
    completions(_a...; __k...) = String[]
else
    function completions(string, pos, context_module = Main)
        ret, _, should_complete =
            REPL.completions(string, pos, context_module)
        if should_complete
            return map(REPL.completion_text, ret)
        else
            return String[]
        end
    end
end

end  # module


text:
export MSet

mutable struct MSet{T} <: AbstractSet{T}
  dict::Dict{T, Int}
  MSet{T}() where {T} = new(Dict{T,Int}())
  MSet{T}(itr) where {T} = union!(new(Dict{T,Int}()), itr)
end

MSet() = MSet{Any}()
MSet(itr) = MSet{eltype(itr)}(itr)

Base.eltype(::Type{MSet{T}}) where {T} = T
Base.similar(s::MSet{T}) where {T} = MSet{T}()
Base.similar(s::MSet, T::Type) = MSet{T}()

#TODO: compact printing, remove trailing , ... the works...
function Base.show(io::IO, s::MSet)
    print(io,"MSet")
    if isempty(s)
        print(io,"{",eltype(s),"}()")
        return
    end
    print(io,"(")
    for (k,v) = s.dict
      if v > 1
        print(io, "$k : $v, ") # ugly
      else
        print(io, "$k, ") # ugly
      end
    end
    print(io,")")
end

Base.isempty(s::MSet) = isempty(s.dict)
Base.length(s::MSet)  = BigInt(sum(values(s.dict)))
Base.IteratorSize(s::MSet) = Base.HasLength()
Base.IteratorEltype(s::MSet) = Base.HasEltype()
Base.eltype(s::MSet{T}) where {T} = T
Base.in(x, s::MSet) = haskey(s.dict, x)
function Base.push!(s::MSet, x, mult::Int = 1)
  if haskey(s.dict, x)
    s.dict[x] += mult
  else
    s.dict[x] = mult
  end
end

function Base.pop!(s::MSet, x)
  s.dict[x] -= 1
  if s.dict[x] == 0
    delete!(s.dict, x)
  end
  return x
end
Base.pop!(s::MSet, x, deflt) = x in s ? pop!(s, x) : deflt
Base.pop!(s::MSet) = (val = iterate(s.dict)[1][1]; pop!(s, val))
Base.delete!(s::MSet, x) = (delete!(s.dict, x); s)

Base.copy(s::MSet) = union!(similar(s), s)

function Base.iterate(s::MSet)
  I = iterate(s.dict)
  I === nothing && return I
  return I[1][1], (I[1], I[2], 1)
end

function Base.iterate(s::MSet, state)
  if state[3] < state[1][2]
    return state[1][1], (state[1], state[2], state[3]+1)
  else
    I = iterate(s.dict, state[2])
    I === nothing && return I
    val, st = I
    return (val[1], (val, st, 1))
  end
end

Base.union(s::MSet) = copy(s)
function Base.union(s::MSet, sets::Set...)
    u = MSet{join_eltype(s, sets...)}()
    union!(u,s)
    for t in sets
        union!(u,t)
    end
    return u
end
Base.union!(s::MSet, xs) = (for x=xs; push!(s,x); end; s)
Base.union!(s::MSet, xs::AbstractArray) = (for x=xs; push!(s,x); end; s)

function multiplicities(s::MSet)
  return values(s.dict)
end

function Base.unique(s::MSet)
  return collect(keys(s.dict))
end

############################################
# subsets iterator
############################################

struct MSubSetItr{T}
  b::Array{T, 1}
  m::Array{Int, 1}
  length::Int
end

function subsets(s::MSet{T}) where T
  # subsets are represented by integers in a weird base
  # the distinct elements are b0...bn with mult mi
  # subset (bi, ni) -> sum ni gi where gi = prod (mj+1)
  b = collect(unique(s))
  m = [s.dict[x] for x = b]
  #= not needed for the iterator
  g = [1]
  for i=2:length(b)
    push!(g, g[end]*(m[i]+1))
  end
  =#
  return MSubSetItr{T}(b, m, length(m) == 0 ? 1 : prod(x+1 for x=m))
end

function int_to_elt(M::MSubSetItr{T}, i::Int) where T
  s = MSet{T}()
  for j=1:length(M.b)
    k = i % (M.m[j]+1)
    for l=1:k
      push!(s, M.b[j])
    end
    i = div(i-k, M.m[j]+1)
  end
  return s
end

function Base.iterate(M::MSubSetItr)
  return int_to_elt(M, 0), 0
end

function Base.iterate(M::MSubSetItr, st::Int)
  st += 1
  st >= M.length && return nothing
  return int_to_elt(M, st), st
end

function Base.length(M::MSubSetItr)
  return M.length
end

Base.IteratorSize(M::MSubSetItr) = Base.HasLength()
Base.IteratorEltype(M::MSubSetItr) = Base.HasEltype()
Base.eltype(M::MSubSetItr{T}) where {T} = MSet{T}

function Base.show(io::IO, M::MSubSetItr)
  println(io, "subset iterator of length $(M.length) for $(M.b) with multiplicities $(M.m)")
end

#... to be completed from base/Set.jl ...


text:
module AVDevice
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(avdevice_dir, f)

  using ..AVUtil
  using ..AVFormat
  using ..AVCodecs
  include(w("LIBAVDEVICE.jl"))
end


text:
"""
    Arf <: Real
"""
struct Arf <: Real
    arf::arf_struct
    prec::Int

    Arf(; prec::Integer = DEFAULT_PRECISION[]) = new(arf_struct(), prec)

    Arf(x::Union{UInt,Int}; prec::Integer = DEFAULT_PRECISION[]) = new(arf_struct(x), prec)
end

"""
    Mag <: Real
"""
struct Mag <: Real
    mag::mag_struct

    Mag() = new(mag_struct())

    Mag(x::mag_struct) = new(mag_struct(x))

    Mag(x::Union{Mag,Arf}) = new(mag_struct(cstruct(x)))
end

"""
    Arb <: Real
"""
struct Arb <: Real
    arb::arb_struct
    prec::Int

    Arb(; prec::Integer = DEFAULT_PRECISION[]) = new(arb_struct(), prec)
end

"""
    Acb <: Number
"""
struct Acb <: Number
    acb::acb_struct
    prec::Int

    Acb(; prec::Integer = DEFAULT_PRECISION[]) = new(acb_struct(), prec)
end

# Refs are in reverse order to model their possible depencies
"""
    AcbRef <: Number
"""
struct AcbRef <: Number
    acb_ptr::Ptr{acb_struct}
    prec::Int
    parent::Union{Nothing,acb_vec_struct,acb_mat_struct}
end

"""
    ArbRef <: Real
"""
struct ArbRef <: Real
    arb_ptr::Ptr{arb_struct}
    prec::Int
    parent::Union{acb_struct,AcbRef,arb_vec_struct,arb_mat_struct}
end

"""
    ArfRef <: Real
"""
struct ArfRef <: Real
    arf_ptr::Ptr{arf_struct}
    prec::Int
    parent::Union{arb_struct,ArbRef}
end

"""
    MagRef <: Real
"""
struct MagRef <: Real
    mag_ptr::Ptr{mag_struct}
    parent::Union{arb_struct,ArbRef}
end

"""
    ArbPoly
"""
struct ArbPoly
    arb_poly::arb_poly_struct
    prec::Int

    ArbPoly(; prec::Integer = DEFAULT_PRECISION[]) = new(arb_poly_struct(), prec)
end

"""
    ArbSeries <: Number
"""
struct ArbSeries <: Number
    arb_poly::arb_poly_struct
    degree::Int
    prec::Int

    ArbSeries(degree::Integer; prec::Integer = DEFAULT_PRECISION[]) =
        new(arb_poly_struct(), degree, prec)
end

"""
    AcbPoly
"""
struct AcbPoly
    acb_poly::acb_poly_struct
    prec::Int

    AcbPoly(; prec::Integer = DEFAULT_PRECISION[]) = new(acb_poly_struct(), prec)
end

"""
    AcbSeries <: Number
"""
struct AcbSeries <: Number
    acb_poly::acb_poly_struct
    degree::Int
    prec::Int

    AcbSeries(degree::Integer; prec::Integer = DEFAULT_PRECISION[]) =
        new(acb_poly_struct(), degree, prec)
end

"""
    ArbVector <: DenseVector{Arb}
"""
struct ArbVector <: DenseVector{Arb}
    arb_vec::arb_vec_struct
    prec::Int
end
ArbVector(n::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    ArbVector(arb_vec_struct(n), prec)

"""
    AcbVector <: DenseVector{Acb}
"""
struct AcbVector <: DenseVector{Acb}
    acb_vec::acb_vec_struct
    prec::Int
end
AcbVector(n::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    AcbVector(acb_vec_struct(n), prec)

"""
    ArbMatrix <: DenseMatrix{Arb}
"""
struct ArbMatrix <: DenseMatrix{Arb}
    arb_mat::arb_mat_struct
    prec::Int
end
ArbMatrix(r::Integer, c::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    ArbMatrix(arb_mat_struct(r, c), prec)

"""
    AcbMatrix <: DenseMatrix{Acb}
"""
struct AcbMatrix <: DenseMatrix{Acb}
    acb_mat::acb_mat_struct
    prec::Int
end
AcbMatrix(r::Integer, c::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    AcbMatrix(acb_mat_struct(r, c), prec)

"""
    ArbRefVector <: DenseMatrix{ArbRef}
"""
struct ArbRefVector <: DenseVector{ArbRef}
    arb_vec::arb_vec_struct
    prec::Int
end
ArbRefVector(n::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    ArbRefVector(arb_vec_struct(n), prec)

"""
    AcbRefVector <: DenseMatrix{AcbRef}
"""
struct AcbRefVector <: DenseVector{AcbRef}
    acb_vec::acb_vec_struct
    prec::Int
end
AcbRefVector(n::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    AcbRefVector(acb_vec_struct(n), prec)

"""
    ArbRefMatrix <: DenseMatrix{ArbRef}
"""
struct ArbRefMatrix <: DenseMatrix{ArbRef}
    arb_mat::arb_mat_struct
    prec::Int
end
ArbRefMatrix(r::Integer, c::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    ArbRefMatrix(arb_mat_struct(r, c), prec)

"""
    AcbRefMatrix <: DenseMatrix{AcbRef}
"""
struct AcbRefMatrix <: DenseMatrix{AcbRef}
    acb_mat::acb_mat_struct
    prec::Int
end
AcbRefMatrix(r::Integer, c::Integer; prec::Integer = DEFAULT_PRECISION[]) =
    AcbRefMatrix(acb_mat_struct(r, c), prec)

# conversion between ref and non-ref arrays.
for T in [:Arb, :Acb], A in [:Vector, :Matrix]
    TA = Symbol(T, A)
    TRefA = Symbol(T, :Ref, A)
    @eval begin
        $TRefA(M::$TA) = $TRefA(cstruct(M), precision(M))
        $TA(M::$TRefA) = $TA(cstruct(M), precision(M))
    end
end

for (T, prefix) in (
    (Mag, :mag),
    (Arf, :arf),
    (Arb, :arb),
    (Acb, :acb),
    (Union{ArbVector,ArbRefVector}, :arb_vec),
    (Union{AcbVector,AcbRefVector}, :acb_vec),
    (Union{ArbMatrix,ArbRefMatrix}, :arb_mat),
    (Union{AcbMatrix,AcbRefMatrix}, :acb_mat),
    (Union{ArbPoly,ArbSeries}, :arb_poly),
    (Union{AcbPoly,AcbSeries}, :acb_poly),
)
    arbstruct = Symbol(prefix, :_struct)
    @eval begin
        cprefix(::Type{<:$T}) = $(QuoteNode(Symbol(prefix)))
        cstruct(x::$T) = getfield(x, cprefix($T))
        cstruct(x::$arbstruct) = x
        cstructtype(::Type{<:$T}) = $arbstruct
        Base.convert(::Type{$arbstruct}, x::$T) = cstruct(x)
    end
end

# handle Ref types
for prefix in [:mag, :arf, :arb, :acb]
    T = Symbol(uppercasefirst(string(prefix)))
    TRef = Symbol(T, :Ref)
    TStruct = Symbol(prefix, :_struct)
    TPtr = Symbol(prefix, :_ptr)
    @eval begin
        cprefix(::Type{$TRef}) = $(QuoteNode(prefix))
        cstructtype(::Type{$TRef}) = $TStruct
        cstruct(x::$TRef) = x.$TPtr
        Base.convert(::Type{Ptr{$TStruct}}, x::$TRef) = cstruct(x)
        Base.cconvert(::Type{Ref{$TStruct}}, x::$TRef) = cstruct(x)

        _nonreftype(::Type{<:Union{$T,$TRef}}) = $T

        parentstruct(x::$T) = cstruct(x)
        parentstruct(x::$TRef) = x
        parentstruct(x::$TStruct) = x
    end
end

const MagOrRef = Union{Mag,MagRef}
const ArfOrRef = Union{Arf,ArfRef}
const ArbOrRef = Union{Arb,ArbRef}
const AcbOrRef = Union{Acb,AcbRef}
const MagLike = Union{Mag,MagRef,cstructtype(Mag),Ptr{cstructtype(Mag)}}
const ArfLike = Union{Arf,ArfRef,cstructtype(Arf),Ptr{cstructtype(Arf)}}
const ArbLike = Union{Arb,ArbRef,cstructtype(Arb),Ptr{cstructtype(Arb)}}
const AcbLike = Union{Acb,AcbRef,cstructtype(Acb),Ptr{cstructtype(Acb)}}
const ArbVectorLike = Union{ArbVector,ArbRefVector,cstructtype(ArbVector)}
const AcbVectorLike = Union{AcbVector,AcbRefVector,cstructtype(AcbVector)}
const ArbMatrixLike = Union{ArbMatrix,ArbRefMatrix,cstructtype(ArbMatrix)}
const AcbMatrixLike = Union{AcbMatrix,AcbRefMatrix,cstructtype(AcbMatrix)}
const ArbPolyLike = Union{ArbPoly,ArbSeries,cstructtype(ArbPoly)}
const AcbPolyLike = Union{AcbPoly,AcbSeries,cstructtype(AcbPoly)}
const ArbTypes = Union{
    Mag,
    MagRef,
    Arf,
    ArfRef,
    Arb,
    ArbRef,
    Acb,
    AcbRef,
    ArbVector,
    ArbRefVector,
    AcbVector,
    AcbRefVector,
    ArbMatrix,
    ArbRefMatrix,
    AcbMatrix,
    AcbRefMatrix,
    ArbPoly,
    ArbSeries,
    AcbPoly,
    AcbSeries,
}

Base.setindex!(x::Union{Mag,MagRef,Arf,ArfRef,Arb,ArbRef,Acb,AcbRef}, z::Number) =
    set!(x, z)
Base.setindex!(x::MagOrRef, z::Ptr{mag_struct}) = set!(x, z)
Base.setindex!(x::ArfOrRef, z::Ptr{arf_struct}) = set!(x, z)
Base.setindex!(x::ArbOrRef, z::Ptr{arb_struct}) = set!(x, z)
Base.setindex!(x::AcbOrRef, z::Ptr{acb_struct}) = set!(x, z)


text:
############################
#### Homework No 1 #########
#### Ahmed El Fatmaoui #####
############################

#cd("/Users/ahmedelfatmaoui/Documents/MetricsIII/HOW1")
#---------------
#-----Q1--------
#---------------

using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Distributions, Kronecker

function q1()
    #a)
         Random.seed!(1234)

         A = rand(Uniform(-5, 10),10,7)
         B = rand(Normal(-2, 15),10,7)
         C = [A[1:5,1:5] B[1:5,6:7]]
         D = (A .<= 0) .* A
    #d-g)
         elm = size(A)
         unqElm = unique(D)
         E = reshape(B,70,1)
         E2 = vec(B)
         F =  reshape([A B], 10, 7, 2)
         F = permutedims(F, [3, 1, 2])
         G = kronecker(B,C)
        # G2 = kronecker(C,F)
    #g-k)
    save("matrixpractice.jld", "A", A, "B", B, "C", C,"D", D, "E", E, "F", F, "G", G )
    save("firstmatrix.jld", "A", A, "B", B, "C", C,"D", D )
    CSV.write("Cmatrix.csv",DataFrame(C,:auto))
    CSV.write("Dmatrix.dat",DataFrame(D,:auto), delim = ' ')

     return A, B, C, D

end

A, B, C, D = q1()

#---------------
#-----Q2--------
#---------------

function q2(A, B, C)
    #a)
  AB = [A[i,j]*B[i,j] for i in 1:10, j in 1:7]

  #b)
  Cprime =  [C[i,j] for i in 1:5, j in 1:7 if -5 <= C[i,j] <= 5]

  C[-5 .<= C .<= 5]

  #c)
  N, K, T = 15169, 6, 5

  X = zeros(N,K,T);

  x5 = rand(Binomial(20,0.6),N)
  x6 = rand(Binomial(20,0.5),N)
      for k in 1:T

          X[:,1,k] = ones(N)

          if rand(1)[1] <= 75*(6-k)/5
              X[:,2,k] = ones(N)
          end

          X[:,3,k] = rand(Normal(15 + k - 1, 5*(k-1)),N)

          X[:,4,k] = rand(Normal(pi*(6-k)/3, 1/exp(1)),N)

          # let columns 1, 5 and 6 remain stationary over time.
          X[:,5, k] = x5
          X[:,6, k] = x6

      end

  # d)
  β = zeros(K,T);
  β[1,:] = [1:0.25:2;];
  β[2,:] = [log(j) for j in 1:T];
  β[3,:] = [-sqrt(j) for j in 1:T];
  β[4,:] = [exp(j) - exp(j+1) for j in 1:T];
  β[5,:] = [j for j in 1:T];
  β[6,:] = [j/3 for j in 1:T];

  β

  # e)

 Y = [X[:,:,t] * β[:,t] + rand(Normal(0,0.36),N,1)  for t in 1:T]

end

q2(A,B,C)

#---------------
#-----Q3--------
#---------------
df = CSV.read("nlsw88.csv", DataFrame);
function q3()
    # a)
    df = CSV.read("nlsw88.csv", DataFrame);
    # ismissing(df)
    # [print(string(":", x, ", ")) for x in names(df)]

        # b)
        NotMarrPer = 1 - sum(df[:,"married"])/size(df)[1];
        GradPer = sum(df[:,"collgrad"])/size(df)[1];

        println("Percent not married = ", round(NotMarrPer, digits = 4) * 100,
        "\nCollege graduate percentage = ", round(GradPer, digits =4 ) * 100)

            # c)  note: tabulate command is undefined
            # source of freq. tables:
            # https://stackoverflow.com/questions/34654489/rs-table-function-in-julia-for-dataframes

            prop(freqtable(df.race))

                # d)

                intrq = (describe(df,:q75)[:,2] - describe(df,:q25)[:,2])
                summarystats = [describe(df,:nmissing, :mean, :median, :std, :min, :max, :nunique) intrq]
                rename(summarystats, :x1 => :itrq)

                    # e)
                    ## Cross tab between industry and occupation
                    freqtable(df.industry, df.occupation)

                        # f)
                        # source: https://dataframes.juliadata.org/stable/man/split_apply_combine/

                        q1 = df[:, [:race, :occupation, :wage]];
                        gdf = groupby(q1, [:race, :occupation]);
                        combine(gdf, :wage => mean)

end

q3()

#---------------
#-----Q4--------
#---------------
function q4()
        fmat = load("firstmatrix.jld")
        #keys(fmat)
        #fmat["B"]

        # a-e)
        function matrixops(A, B)
            # element by element product of inputs, A and b
            # input matrices product
            # input matrices sum
                    if(size(A) != size(B))
                    println("inputs must have the same size")
                    else
                        return println(
                            "A .* B = ",  A .* B,
                            "\n\nA' * B = ", A' * B,
                            "\n\nA + B = ",  A + B
                        )
                    end
        end

        matrixops(A, B)

        # f)
        # we get the error message back
        matrixops(C, D)

        # g)

        matrixops(
            convert(Array,df.ttl_exp ),
            convert(Array,df.wage )
        )
end
    q4()


text:
struct NormalMesh{U,D1,T} <: AbstractMesh{U,D1,T}
    vertices::Vector{SVector{U,T}}
    faces::Vector{SVector{D1,Int}}
    normals::Vector{Int}
end


text:
# The branching factors for the various INode levels, as determined by
# _config_for_size(N). A NodeConfig remains valid until the number of dictionary
# entries N reaches next_config_size - changing configurations is handled by
# insert!(dict,h,k,v).
struct NodeConfig
    inode_capacity::Vector{Int}
    next_config_size::Int
end

function _config_for_size2(N::Int)::NodeConfig
    if N < 16
        return NodeConfig(Int64[], 16)
    elseif N < 65536
        # For small dictionaries, use inode_capacity=4, and height increasing with N.
        num_bits = ceil(log2(N+1))
        leaf_bits = min(num_bits,4)
        num_inodes = max(1, convert(Int, ceil((num_bits-leaf_bits)/2)))
        inode_capacity = [4 for i in 1:num_inodes]
        next_config_size = 2^(4+2*(num_inodes))
    else
        # For large dictionaries, we use height 7, and inode capacity increasing with
        # N.  We partition ceil(log2(N)) bits among 7 levels, and to grow capacity as
        # N increases we add a single bit to one level.

        # We will have 1/7 of the bits at the leaf node, so 6/7 of the bits will be
        # partitioned among the inodes.  We start by giving all levels the same number
        # of bits, then take care of the remainder by putting one extra bit at deeper
        # levels.
        num_bits = ceil(log2(N+1))
        leaf_bits = max(4, convert(Int, floor(num_bits/7)))
        inode_bits = num_bits - leaf_bits
        m = convert(Int, floor(inode_bits/6))
        inode_capacity = [2^m for i in 1:6]
        have_bits = m*6
        for k=6:-1:1
            if have_bits < inode_bits
                inode_capacity[k] *= 2
                @dassert1 inode_capacity[k] <= 256   # N <= 2^56
                have_bits += 1
            else
                break
            end
        end
        next_config_size = convert(Int, 2^(floor(log2(N))+1))
    end
    shifts = [0 for i in 1:6]
    return NodeConfig(inode_capacity, next_config_size)
end

const config1 = _config_for_size2(1)
const config16 = _config_for_size2(16)
const config64 = _config_for_size2(64)
const config256 = _config_for_size2(256)

function _config_for_size(N::Int)::NodeConfig
    # For N < 1024, use one of the pre-allocated instances config1/config16/
    # config64/config256, to reduce space use by small dictionaries.
    if N < 64
        if N < 16
            @dassert1 config1.next_config_size == 16
            return config1
        else
            @dassert1 config16.next_config_size == 64
            return config16
        end
    elseif N < 1024
        if N < 256
            @dassert1 config64.next_config_size == 256
            return config64
        else
            @dassert1 config256.next_config_size == 1024
            return config256
        end
    else
        return _config_for_size2(N)
    end
end


text:
export plot_1D_solutions_branch

"""Calculate distance between a given state and a stable branch"""
function _closest_branch_index(res::Result,state::Vector{Float64},index::Int64)
    #search only among stable solutions
    physical     = filter_solutions.(
                            filter_solutions.(
                            res.solutions, res.classes["physical"]),res.classes["stable"])

    steadystates = reduce(hcat,physical[index]) #change structure for easier indexing
    distances    = vec(sum(abs2.(steadystates .- state),dims=1))
    return argmin(replace!(distances, NaN=>Inf))
end

"""Return the indexes and values following stable branches along a 1D sweep.
When a no stable solutions are found (e.g. in a bifurcation), the next stable solution is calculated by time evolving the previous solution (quench).

Keyword arguments
- `y`:  Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
- `sweep`: Direction for the sweeping of solutions. A `right` (`left`) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
- `tf`: time to reach steady
- `ϵ`: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)
"""
function follow_branch(starting_branch::Int64,res::Result; y="sqrt(u1^2+v1^2)",sweep="right",tf=10000,ϵ=1E-4)
    sweep_directions = ["left", "right"]
    sweep ∈ sweep_directions  || error("Only the following (1D) sweeping directions are allowed:  ", sweep_directions)

    #filter stable solutions
    Y  = transform_solutions(res, y)
    Ys = filter_solutions.(filter_solutions.(Y, res.classes["physical"]),res.classes["stable"]);
    Ys = real.(reduce(hcat,Ys)) #facilitate indexing. Ys is an array (length(sweep))x(#solutions)

    followed_branch    = zeros(Int64,length(Y)) #followed branch indexes
    followed_branch[1] = starting_branch

    if sweep == "left"
        Ys = reverse(Ys,dims=2)
    end

    p1 = collect(keys(res.swept_parameters))[1] #parameter values

    for i in 2:size(Ys)[2]
        s = Ys[followed_branch[i-1],i] #solution amplitude in the current branch and current parameter index
        if !isnan(s) #the solution is not unstable or unphysical
            followed_branch[i] = followed_branch[i-1]
        else #bifurcation found
            if sweep == "right"
                next_index = i
            elseif sweep == "left" #reverse order of iterators
                next_index = length(Y)-i+1
            end

            # create a synthetic starting point out of an unphysical solution: quench and time evolve
            #the actual solution is complex there, i.e. non physical. Take real part for the quench.
            sol_dict = get_single_solution(res, branch=followed_branch[i-1], index= next_index)
            print("bifurcation @ ", p1 ," = ", real(sol_dict[p1])," ")
            sol_dict   = Dict(zip(keys(sol_dict),real.(values(sol_dict)) .+ 0.0im .+ ϵ*rand(length(values(sol_dict)))))
            problem_t  = TimeEvolution.ODEProblem(res.problem.eom, steady_solution=sol_dict, timespan=(0,tf))
            res_t      = TimeEvolution.solve(problem_t,saveat=tf)

            followed_branch[i] = _closest_branch_index(res,res_t.u[end],next_index) #closest branch to final state

            print("switched branch ", followed_branch[i-1] ," -> ", followed_branch[i],"\n")
        end
    end

    if sweep == "left"
        Ys = reverse(Ys,dims=2)
        followed_branch = reverse(followed_branch)
    end

    return followed_branch,Ys
end

"""1D plot with the followed branch highlighted"""
function plot_1D_solutions_branch(starting_branch::Int64,res::Result; x::String, y::String, sweep="right",tf=10000,ϵ=1E-4, xscale="linear",yscale="linear",plot_only=["physical"],marker_classification="stable",filename=nothing,kwargs...)
    plot_1D_solutions(res; x=x, y=y, xscale=xscale,yscale=yscale,
                        plot_only=plot_only,marker_classification=marker_classification,filename=filename,kwargs...)

    followed_branch,Ys = follow_branch(starting_branch,res,y=y,sweep=sweep,tf=tf,ϵ=ϵ)
    if sweep == "left"
        m = "<"
    elseif sweep == "right"
        m = ">"
    end

    Y_followed = [Ys[branch,param_idx] for (param_idx,branch) in enumerate(followed_branch)]

    lines = plot(transform_solutions(res, x),Y_followed,c="k",marker=m; kwargs...)

    #extract plotted data and return it
    xdata,ydata = [line.get_xdata() for line in lines], [line.get_ydata() for line in lines]
    markers = [line.get_marker() for line in lines]
    save_dict= Dict(string(x) => xdata,string(y)=>ydata,"markers"=>markers)

    !isnothing(filename) ? JLD2.save(_jld2_name(filename), save_dict) : nothing
    return save_dict
end


text:
using HMCExamples

HMCExamples.main_rbc_2_joint(ARGS)

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
@testset "bounded_by" begin
    demo_problems = demo_problem.(1:length(demo_problems_definitions))

    @testset "test problems" begin
        i = 1
        for (f, (a, b), _, fmax, _, fabsmax) in demo_problems
            # We check that the bound is satisfied for a value
            # slightly larger than the maximum but not satisfied for a
            # value slightly smaller. What slightly means depends on
            # the degree we use, sqrt(eps(fmax)) for the default and
            # abs(fmax) / 1000 for degree = -1.
            @test ArbExtras.bounded_by(f, a, b, ubound(fmax + sqrt(eps(fmax))))
            @test !ArbExtras.bounded_by(f, a, b, lbound(fmax - sqrt(eps(fmax))))

            @test ArbExtras.bounded_by(f, a, b, ubound(fmax + abs(fmax) / 100), degree = -1)
            @test !ArbExtras.bounded_by(
                f,
                a,
                b,
                lbound(fmax - abs(fmax) / 100),
                degree = -1,
            )

            @test ArbExtras.bounded_by(
                f,
                a,
                b,
                ubound(fabsmax + sqrt(eps(fabsmax))),
                abs_value = true,
            )
            @test !ArbExtras.bounded_by(
                f,
                a,
                b,
                lbound(fabsmax - sqrt(eps(fabsmax))),
                abs_value = true,
            )

            @test ArbExtras.bounded_by(
                f,
                a,
                b,
                ubound(fabsmax + abs(fabsmax) / 100),
                abs_value = true,
                degree = -1,
            )
            @test !ArbExtras.bounded_by(
                f,
                a,
                b,
                lbound(fabsmax - abs(fabsmax) / 100),
                abs_value = true,
                degree = -1,
            )
        end
    end

    @testset "special cases" begin
        # Thin interval
        @test ArbExtras.bounded_by(x -> 2x, Arf(2), Arf(2), Arf(4))
        @test !ArbExtras.bounded_by(x -> 2x, Arf(2), Arf(2), Arf(3))
        @test ArbExtras.bounded_by(x -> -2x, Arf(2), Arf(2), Arf(4), abs_value = true)
        @test !ArbExtras.bounded_by(x -> -2x, Arf(2), Arf(2), Arf(3), abs_value = true)

        # Non-valid endpoints
        @test_throws ArgumentError ArbExtras.bounded_by(identity, Arf(2), Arf(1), Arf(0))

        # Threading enabled
        @test ArbExtras.bounded_by(x -> 2x, Arf(2), Arf(3), Arf(6.1), threaded = true)
        @test !ArbExtras.bounded_by(x -> 2x, Arf(2), Arf(3), Arf(5.9), threaded = true)

        # verbose = true
        @test_logs (:warn, "reached maximum depth 20") match_mode = :any ArbExtras.bounded_by(
            acos,
            Arf(0),
            Arf(0.5),
            ubound(Arb(π) / 2),
            degree = -1,
            verbose = true,
        )

        @test_logs (:warn, "reached maximum number of evaluations 101 >= 100") match_mode =
            :any ArbExtras.bounded_by(
            sin ∘ exp,
            Arf(0),
            Arf(8),
            Arf(0.999),
            degree = -1,
            maxevals = 100,
            verbose = true,
        )
    end

end


text:
function iterate!(state::State, idx)
  for I in CartesianIndices(state.res)
    val = state.op(state.field[I], state.deltas[I])
    if (abs2(val) <= 65536)
      state.fieldIterations[I] = idx
      state.field[I] = val
    end
  end
end

function step!(state::State)
  for i = 1:state.iterations
    iterate!(state, i)
  end
end

max_stars_repo_name: benchislett/Skunkworks

max_issues_repo_name: benchislett/Skunkworks

max_forks_repo_name: benchislett/Skunkworks


text:
using PkgSRA
using Plots, Random, Distributions, Interpolations
pyplot()
Random.seed!(11)
using BSON: @save
include("../scripts/paper_settings.jl");
# include("../../utils/sindy_turing_utils.jl")
# include("../../utils/main_algorithm_utils.jl")
# include("../../utils/sindy_statistics_utils.jl")

################################################################################
#####
##### Van der Pol example: get data
#####
include(EXAMPLE_FOLDERNAME*"example_vanDerPol.jl")

# Define the multivariate forcing function
num_ctr = 3;
    U_starts = rand(1, num_ctr) .* tspan[2]/2
    U_widths = 0.2;
    # U_widths = 0;
    amplitude = 2.0
my_U_func_time2(t) = U_func_time(t, u0,
                        U_widths, U_starts,
                        F_dim=1,
                        amplitude=amplitude)

# Get data
sol = solve_vdp_system(U_func_time=my_U_func_time2)
dat = Array(sol)
numerical_grad = numerical_derivative(dat, ts)
true_grad = core_dyn_true(dat)

U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time2(t)
end

## Also get baseline true/ideal cases
# Uncontrolled
dat_raw = Array(solve_vdp_system())
numerical_grad_raw = numerical_derivative(dat_raw, ts)
true_grad_raw = zeros(size(dat))

# Intialize truth object
this_truth = sra_truth_object(dat_raw, true_grad, U_true, core_dyn_true)

#####
##### Build SRA object and analyze
#####
# Initialize
this_model = sra_stateful_object(ts, tspan, dat, u0, numerical_grad)
prams = this_model.parameters
prams.sindyc_ensemble_parameters[:selection_criterion] =
    sindy_cross_validate;
# Several changes are required because there are only 2 variables here
# p.sindyc_ensemble_parameters[:variable_names] = ["x", "y"];
prams.variable_names = ["x", "y"];
prams.sindy_terms_list = Iterators.product(1:3, 1:3)
prams.sindy_library["cross_terms"] = [2, 3] # Also include cubic terms

fit_first_model(this_model, 30);
print_true_equations(this_truth)
print_current_equations(this_model)

#################################################################
### Iterate
calculate_subsampled_ind(this_model);
all_models = fit_model(this_model);

print_true_equations(this_truth)
print_current_equations(this_model)

# plot_subsampled_points(this_model)
# plot_subsampled_simulation(this_model, 2)
# plot_subsampled_derivatives(this_model, 1)
# plot_residual(this_model)

# plot_subsampled_points_and_control(this_model, this_truth)

#####
##### Save Van der Pol data
#####
this_dat_name = "dat_library_of_examples_vdp_"

save_for_plotting(this_model, this_truth, this_dat_name)


text:
function _read_bash(src::String;
        verbose = true,
        cmdkwargs...
    )
    cmd = Cmd(`bash -c $(src)`; cmdkwargs...)
    out = read(cmd, String)
    verbose && println(out)
    return out
end

function _spawn_bash(src::String; cmdkwargs...)
    cmd = Cmd(`bash -c $(src)`; cmdkwargs...)
    proc = run(cmd; wait = false)
    return _try_getpid(proc)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using LoggingExtras
using Test
using Test: collect_test_logs, TestLogger
using Dates

using Base.CoreLogging
using Base.CoreLogging: BelowMinLevel, Debug, Info, Warn, Error

@testset "Tee" begin
    @testset "contructor" begin
        @testset "mixed types" begin
            @test TeeLogger(TestLogger(), NullLogger()) isa TeeLogger
        end

        @testset "errors if given nonloggers" begin
            @test_throws Exception TeeLogger(stdout, stderr)
        end
    end
    @testset "basic use with compositional levels" begin
        testlogger_info = TestLogger(min_level=Info)
        testlogger_warn = TestLogger(min_level=Warn)

        with_logger(TeeLogger(testlogger_warn, testlogger_info)) do
            @info "info1"
            @warn "warn1"
            @info "info2"
        end
        @test length(testlogger_info.logs) == 3
        @test length(testlogger_warn.logs) == 1
    end
end

@testset "File" begin
    mktempdir() do dir
        for (filepath, sink) in [
                (f = joinpath(dir, "log"); (f, f)), # Filepath
                (f = joinpath(dir, "log_io"); (f, open(f, "w"))), # IOStream
            ]
            with_logger(FileLogger(sink)) do
                @info "first"
                @warn "second"
                @info "third"
            end
            sink isa IOStream && close(sink)
            logtext = read(filepath, String)
            @test occursin("first", logtext)
            @test occursin("second", logtext)
            @test occursin("third", logtext)
        end
    end
end

@testset "Active Filter" begin
    testlogger = TestLogger()
    yodawg_filter(logargs) = startswith(logargs.message, "Yo Dawg!")

    filtered_logger = ActiveFilteredLogger(yodawg_filter, testlogger)

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
    end
    @test length(testlogger.logs) == 2
end

@testset "Early Filter" begin
    testlogger = TestLogger()
    filtered_logger = EarlyFilteredLogger(testlogger) do logargs
        logargs.level == Info  # Only exactly Info, nothing more nothing less
    end

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
    end
    @test length(testlogger.logs) == 3
end

@testset "MinLevel Filter" begin
    testlogger = TestLogger()
    filtered_logger = MinLevelLogger(testlogger, Warn)

    with_logger(filtered_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
        @error "MISTAKES WERE MADE"
    end
    @test length(testlogger.logs) == 2
end

@testset "Transformer" begin
    testlogger = TestLogger(min_level=Error)
    transformer_logger = TransformerLogger(testlogger) do log_msg
        # We are going to transform all warnings into errors
        if log_msg.level == Warn
            return merge(log_msg, (; level=Error))
        else
            return log_msg
        end
    end

    with_logger(transformer_logger) do
        @info "info1"
        @warn "Yo Dawg! It is a warning"
        @info "info2"
        @info "Yo Dawg! It's all good"
        @info "info 3"
        @error "MISTAKES WERE MADE"
    end
    @test length(testlogger.logs) == 2
end

@testset "DatetimeRotatingFileLogger" begin
    mktempdir() do dir
        drfl_min = DatetimeRotatingFileLogger(dir, raw"\m\i\n-YYYY-mm-dd-HH-MM.\l\o\g")
        drfl_hour = DatetimeRotatingFileLogger(dir, raw"\h\o\u\r-YYYY-mm-dd-HH.\l\o\g")
        func = (io, args) -> println(io, reverse(args.message))
        drfl_fmt = DatetimeRotatingFileLogger(func, dir, raw"\f\m\t-YYYY-mm-dd-HH-MM.\l\o\g")
        callback_record = []
        callback(f) = filesize(f) > 0 && push!(callback_record, f)
        drfl_cb = DatetimeRotatingFileLogger(dir, raw"\c\o\m\p-YYYY-mm-dd-HH-MM.\l\o\g";
                                             rotation_callback=callback)

        sink = TeeLogger(drfl_min, drfl_hour, drfl_fmt, drfl_cb)
        with_logger(sink) do
            # Make sure to trigger one minute-level rotation and no hour-level
            # rotation by sleeping until HH:MM:55
            n = now()
            sleeptime = mod((55 - second(n)), 60)
            minute(n) == 59 && (sleeptime += 60)
            sleep(sleeptime)
            @info "first"
            @info "second"
            sleep(10) # Should rotate to next minute
            @info("third")
        end

        # Drop anything that's not a .log file or empty
        files = sort(map(f -> joinpath(dir, f), readdir(dir)))
        files = filter(f -> endswith(f, ".log") && filesize(f) > 0, files)
        min_files = filter(f -> startswith(basename(f), "min-"), files)
        @test length(min_files) == 2

        hour_files = filter(f -> startswith(basename(f), "hour-"), files)
        @test length(hour_files) == 1

        fmt_files = filter(f -> startswith(basename(f), "fmt-"), files)
        @test length(fmt_files) == 2

        # Two files exist, but just one have been rotated
        @test length(callback_record) == 1
        @test occursin(r"comp-\d{4}(-\d{2}){4}\.log$", callback_record[1])

        min1_data = String(read(min_files[1]))
        @test occursin("first", min1_data)
        @test occursin("second", min1_data)
        min2_data = String(read(min_files[2]))
        @test occursin("third", min2_data)

        min_data = String(read(hour_files[1]))
        @test occursin("first", min_data)
        @test occursin("second", min_data)
        @test occursin("third", min_data)

        fmt_data = String(read(fmt_files[1]))
        @test occursin("tsrif", fmt_data)
        @test occursin("dnoces", fmt_data)
        fmt_data = String(read(fmt_files[2]))
        @test occursin("driht", fmt_data)

        # Sub-minute resolution not allowed
        @test_throws(ArgumentError("rotating the logger with sub-minute resolution not supported"),
                     DatetimeRotatingFileLogger(dir, "HH-MM-SS"))

        # Test constructors with pattern as a DateFormat
        l = DatetimeRotatingFileLogger(dir, raw"yyyy-mm-dd.\l\o\g")
        l1 = DatetimeRotatingFileLogger(dir, dateformat"yyyy-mm-dd.\l\o\g")
        l2 = DatetimeRotatingFileLogger(identity, dir, dateformat"yyyy-mm-dd.\l\o\g")
        @test l.filename_pattern == l1.filename_pattern == l2.filename_pattern
    end
end

@testset "FormatLogger" begin
    io = IOBuffer()
    logger = FormatLogger(io) do io, args
        # Put in some bogus sleep calls just to test that
        # log records writes in one go
        print(io, args.level)
        sleep(rand())
        print(io, ": ")
        sleep(rand())
        println(io, args.message)
    end
    with_logger(logger) do
        @sync begin
            @async @debug "debug message"
            @async @info "info message"
            @async @warn "warning message"
            @async @error "error message"
        end
    end
    str = String(take!(io))
    @test occursin(r"^Debug: debug message$"m, str)
    @test occursin(r"^Info: info message$"m, str)
    @test occursin(r"^Warn: warning message$"m, str)
    @test occursin(r"^Error: error message$"m, str)
    @test logger.always_flush
    # Test constructor with default io and kwarg
    logger = FormatLogger(x -> x; always_flush=false)
    @test logger.stream === stderr
    @test !logger.always_flush
end

@testset "Deprecations" begin
    testlogger = TestLogger(min_level=BelowMinLevel)

    @test_logs (:warn, r"deprecated") match_mode=:any begin
        demux_logger = DemuxLogger(testlogger)
        @test demux_logger isa TeeLogger
        @test Set(demux_logger.loggers) == Set([testlogger, global_logger()])
    end

    @test_logs (:warn, r"deprecated") match_mode=:any begin
        demux_logger = DemuxLogger(testlogger; include_current_global=true)
        @test demux_logger isa TeeLogger
        @test Set(demux_logger.loggers) == Set([testlogger, global_logger()])
    end

    @test_logs (:warn, r"deprecated") match_mode=:any begin
        demux_logger = DemuxLogger(testlogger; include_current_global=false)
        @test demux_logger isa TeeLogger
        @test Set(demux_logger.loggers) == Set([testlogger])
    end
end