text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences, FASTX, XLSX, CPUTime

#########
#The following functions are the base set of tools that this program uses

# 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

function is_palindrome(randomer)
    return occursin(randomer,reverse_complement(randomer))
end
"""
    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

"""
    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
function iscompatible_sequence(seqa,seqb)
    n=length(seqa)
    count=0
    for i=1:n
        if iscompatible(seqa[i],seqb[i])
            count+=1;
        end
    end
    return count==n
end

function decompress_oligo(comp_oligo)
    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(comp_oligo;uselog2=false))
    randomerlen=length(comp_oligo)
    blockarray=zeros(Int128,randomerlen+1)
    blockarray[1]=nseqs
    degsarray=zeros(Int128,randomerlen)
    repsarray=zeros(Int128,randomerlen+1)
    repsarray[1]=1
    uncomp_oligos=[dna"-"^randomerlen for _ in 1:nseqs]
    for i = 1:randomerlen
        nucleotides=basedict[comp_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
            uncomp_oligos[m][i]=bases[m];
        end
    end
    return uncomp_oligos
end
function decompress_pool(compressed_pool)
    uncomppool=decompress_oligo(compressed_pool[1])
    for i = 2:length(compressed_pool)
        subset=decompress_oligo(compressed_pool[i])
        for i in subset
            push!(uncomppool,i)
        end
    end
    return uncomppool
end

#####################
#Define the action space
#####################
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

function valid_candidates(sites)
    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_sequence(j,sites[i]),candidates)]
    end
    return candidates
end

function define_action_space(sites)
    if length(sites)==0
        bases=[dna"" for i in 1:6]
        return bases
    else
        randomerlen=length(sites[1]);
        nsites=length(sites)
        bases=[];
        for i=1:randomerlen
            states=dna"-"^nsites;
            nucleotides=dna"----"
            #println("TEST: $nucleotides")
            for j=1:nsites
                states[j]=sites[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
function omit_codes(bases)
    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

#######################
#The follwoing functions are used as part of the various reward systems
######################

#SCORING SYSTEM FOR HITTING ALL mRNAs. Not used in the weighted scoring system
function mRNA_coverage(randomer,mRNA)
    hits=0;
    for seq in mRNA
        if occursin(randomer,seq)
            hits+=1;
        end
    end
    return hits
end

#Weighted Scoring System for mRNA coverage
function mRNA_weighted(randomer,mRNA,weights)
score=0;
for i=1:length(mRNA)
    if occursin(randomer,mRNA[i])
        score+=weights[i];
    end
end
return score
end

# Keeps track of the number of times a given gene has been hit. It is used to compute the weights above
function mRNA_hitcounter(randomer,mRNA,hitcounts)
    for i=1:length(mRNA)
        if occursin(randomer,mRNA[i])
            hitcounts[i]+=1;
        end
    end
    return hitcounts
end
function mRNA_uniform_coverage(randomer,mRNAs;hitcounts=zeros(length(mRNAs)),uniformity_weight=1,total_hits_weight=1,overwrite_hitcounts=false, kwargs...)
    mRNAnum=length(mRNAs);
    gene_scores=zeros(mRNAnum);
    dist_score=0;
    randomerlen=length(randomer);
    new_hitcounts=zeros(mRNAnum)
    #println("hitcounts: $hitcounts")
    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

        new_hitcounts[i]=counter
        hit_score=hitcounts[i]+new_hitcounts[i]
        if hit_score==0
            gene_scores[i]=0;
        else
            gene_scores[i]=1+total_hits_weight*log10(hit_score)
        end
    end
    scored_hits=hitcounts+new_hitcounts;
    if sum(scored_hits)==0
        final_score=0
    else
        sorted_hits=sort(scored_hits);
        DUavg=mean(sorted_hits)
        lq=round(Integer,mRNAnum/4);
        DUlq=mean(sorted_hits[1:lq])
        dist_score=DUlq/DUavg;
        final_score=sum(gene_scores)+uniformity_weight*mRNAnum*dist_score;
    end
    if overwrite_hitcounts
        return final_score,scored_hits,dist_score, uniformity_weight,total_hits_weight
    else
        return final_score,hitcounts,dist_score,uniformity_weight,total_hits_weight
    end
end

function mRNA_multiobjective(randomer, mRNAs; indicies=[[] for i=1:10], b1=1,b2=1,overwrite_indicies=false, return_indicies=false,kwargs...)
    # If overwrite_indicies=false, the funciton will score the randomer as if the randomer were to be added, but will not actually add the randomer
    # If overwrite_indicies=true, then the function will add the randomer hit sites to the offical count
    # b1 = weight for uniformity score
    # b2 = weight for hits per unit length
    mRNAnum=length(mRNAs)
    gene_scores=zeros(mRNAnum);
    means=zeros(mRNAnum);
    stds=ones(mRNAnum);
    randomerlen=length(randomer);
    indicies_copy=deepcopy(indicies);
    for i=1:mRNAnum
        iteration=0
        index=approxsearchindex(mRNAs[i][1:end],randomer,0);
        if index > 0
            push!(indicies_copy[i],index)
        end
        iteration = iteration+1;
        while index > 0
            index=approxsearchindex(mRNAs[i][((indicies_copy[i][end]+1):end)],randomer,0);
            iteration=iteration+1;
            if index == 0
                break
            else
                push!(indicies_copy[i],index+indicies_copy[i][end])

            end
        end
        if length(indicies_copy[i])==0
            continue
        end
        dists=calculate_gap_distances(indicies_copy[i],length(mRNAs[i]),randomerlen)
        stds[i]=std(dists)/length(mRNAs[i])
        gene_scores[i]=1+b1*(1-stds[i])+b2*length(indicies_copy[i])/length(mRNAs[i]);
    end

    if overwrite_indicies && return_indicies
    return indicies_copy, b1, b2
    elseif ~overwrite_indicies &&  return_indicies
    return  indicies
    else  # by having return_indicies=false, you can just return the gene scores. This just makes the output of the function more flexible.
        return sum(gene_scores)
    end
end

function calculate_gap_distances(hit_indicies,genelen,randomerlen)
    dists=[]
    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]-randomerlen))
        end
    end
    push!(dists,(genelen-hit_indicies[end]-randomerlen))
    return dists
end

# Rollout reward systems
#3 factor scoring policy. This is the default reward system
function simulate_random_mRNA_uniform_coverage(prefix,bases,sites,mRNAs; nsims=1000, horizon=length(bases), kwargs...)
    mRNA_score=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);
        mRNA_score[i],b,c,d,e= mRNA_uniform_coverage(sim_randomer_list[i],mRNAs;overwrite_hitcounts=false, kwargs...)
    end
    max_sim_score,max_index=findmax(mRNA_score);
    max_sim_randomer=sim_randomer_list[max_index];
    return mean(mRNA_score),max_sim_randomer,max_sim_score

end
function simulate_random_mRNA_multiobjective(prefix,bases,sites,mRNAs; nsims=1000, horizon=length(bases), kwargs...)
    mRNA_score=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);
        mRNA_score[i]= mRNA_multiobjective(sim_randomer_list[i],mRNAs;return_indicies=false, kwargs...)
    end
    max_sim_score,max_index=findmax(mRNA_score);
    max_sim_randomer=sim_randomer_list[max_index];
    return mean(mRNA_score),max_sim_randomer,max_sim_score
end

#Weighted scoring policy. The score for hitting a given gene decreases according to the number of times it's been hit
function simulate_random_mRNA_weighted(prefix,bases,sites,mRNA; nsims=1000, horizon=length(bases))
    mRNA_score=zeros(nsims)
    true_horizon = min(horizon, length(bases))
    Threads.@threads for i = 1:nsims
        mRNA_score[i]= mRNA_weighted(find_valid_randomer(prefix,bases[1:true_horizon],sites),mRNA,weights)
    end
    return mean(mRNA_score)
end

# Unweighted scoring policy. Finds randomers that hit the most genes at least once
function simulate_random_mRNA(prefix,bases,sites,mRNA; nsims=1000, horizon=length(bases))
    mRNA_hits=zeros(nsims)
    true_horizon = min(horizon, length(bases))
    Threads.@threads for i = 1:nsims
        mRNA_hits[i]= mRNA_coverage(find_valid_randomer(prefix,bases[1:true_horizon],sites),mRNA)
    end
    return mean(mRNA_hits)
end

### MAIN FUNCTION
function rollout_rt(bases,sites, mRNA; simulate=simulate_random_mRNA_multiobjective, 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
    horizon = max_len
    candidates=valid_candidates(sites);
    checkpoint=length(candidates)
    best_sim_randomer=dna"-"^n
    best_sim_score=0;
    #println("Valid Candidates: $checkpoint")
    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_sequence(j[start:i-1],randomer[start:i-1]),candidates)])
        #println("Bases $bases")
        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;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,hitcounts,dist_score,uniformity_weight,total_hits_weight=mRNA_uniform_coverage(randomer,mRNA;kwargs...)
    if best_sim_score >final_randomer_score
        randomer=best_sim_randomer
        println("Triggered Best Simulated Randomer")
    end
    return randomer
end

## 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"))
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 Funcitons

## 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
# genes_hit,genes_missing,total_hits,uniformity_score=hit_stats(indicies,mRNAs,6)
function mRNA_multiobjective_hit_stats(indicies,mRNAs,randomerlen)
    genes_hit=0;
    genes_missing=[]
    total_hits=0;
    stds=ones(length(indicies))
    for i=1:length(indicies);
        if length(indicies[i])>0
            genes_hit+=1;
            dists=calculate_gap_distances(indicies[i],length(mRNAs[i]),randomerlen)
            stds[i]=std(dists)/length(mRNAs[i]);
        else
            push!(genes_missing,i)
        end
        total_hits+=length(indicies[i]);
    end
    return genes_hit,genes_missing,total_hits,mean(1 .- stds)
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.
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,mRNAfile)
    mRNAdata=CSV.read(mRNAfile,DataFrame);
    mRNAs=read_mRNA_list(mRNAdata)
    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
function dnaarray2stringarray(array)
    array2=[]
    for i=1:length(array)
        push!(array2,convert(String,array[i]))
    end
    return array2
end
function stringarray2dnaarray(array)
    array2=[]
    for i=1:length(array)
        push!(array2,LongDNASeq(array[i]))
    end
    return array2
end

function pooled_degeneracy(pool)
    deg=0;
    for randomer in pool
        deg +=degeneracy(randomer;uselog2=false);
    end
    return deg
end
 ########################################################################
 # #CODE FOR CREATING RANDOMER POOLS And Saving them to files
 ########################################################################

 #Quick Test#
 # data=run_NSR_RL(pool_size=3,nsims=50)
 #CSV.write("NSR_RL_S_mutans_4_12_21",data)
function run_NSR_RL(;species="S_mutans",randomerlen=6,pool_size=25,all_bases = dna"AGCTMRWSYKVHDBN",kwargs...)
    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];
    hitcounts=zeros(numgenes)
    times=[]
    num_blocked_sites=zeros(pool_size)
    println("Initialization Successful")
    for i = 1:pool_size
        if i > 1
            push!(blocking_sites,LongDNASeq(randomers[i-1]));
            push!(blocking_sites,LongDNASeq(reverse_complement(randomers[i-1]))) ; ## avoid adding randomers that are reverse complements of each other
            blocking_sites=unique(decompress_pool(blocking_sites));
        end
        num_blocked_sites[i]=length(blocking_sites)
        randomerstats= @timed rollout_rt(all_ns,blocking_sites,mRNAs; hitcounts=hitcounts, simulate=simulate_random_mRNA_uniform_coverage,kwargs...); #create new randomer and add to pool
        push!(randomers,randomerstats.value)
        push!(times,randomerstats.time)
        final_score,hitcounts,dist_score,uniformity_weight,total_hits_weight=mRNA_uniform_coverage(randomers[i],mRNAs; hitcounts=hitcounts, overwrite_hitcounts=true,kwargs...); #officially add hit positions from the new randomer
        #println("Iteration $i")
    end
    final_score,hitcounts,dist_score,uniformity_weight,total_hits_weight =mRNA_uniform_coverage(dna"-"^randomerlen,mRNAs; hitcounts=hitcounts, overwrite_hitcounts=false,kwargs...);
    genes_hit=length(findall(x ->x>0,hitcounts))
    genes_missing=findall(x -> x==0,hitcounts)
    total_hits=sum(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, uniformity_weight=uniformity_weight,hits_weight=total_hits_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,coverage_uniformity_score=dist_score,individual_times=string_time,time=total_time,num_blocked_sites=string_num_blocked_sites)
    return data
end

function benchmark_NSR_RL()
    DF1=DataFrame()
    nexps=1;
    uniformity_weights=[0]
    total_hits_weights=[100]
    species="S_mutans"
    filename="./S_mutans_NSR_Alt_Valid_Randomer_Test_Timing_6_10_21.csv"
    for i = 1:nexps
        data=run_NSR_RL(;species=species,pool_size=3,nsims=100,uniformity_weight=uniformity_weights[i],total_hits_weight=total_hits_weights[i])
        DF1=vcat(DF1,data);
        CSV.write(filename,DF1)
    end
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...)
    mRNAdata=CSV.read(mRNAfile,DataFrame)
    rRNAdata=CSV.read(rRNAfile,DataFrame)
    mRNAs=read_mRNA_list(mRNAdata);
    rRNAs=read_rRNA_tRNA_list(rRNAdata);
    numgenes=length(mRNAs)
    randomerlen=length(randomers[1]);
    blocking_sites=read_blocking_sites(rRNAdata,randomerlen=randomerlen);
    mRNAbp,mRNA_GC=calculate_transcriptome_stats(mRNAs);
    rRNAbp,rRNA_GC=calculate_transcriptome_stats(rRNAs);
    transcriptome_size=mRNAbp+rRNAbp;
    hitcounts=zeros(numgenes)
    dist_score=0
    for i=1:length(randomers)
        final_score,hitcounts,dist_score,uniformity_weight,total_hits_weight=mRNA_uniform_coverage(randomers[i],mRNAs; hitcounts=hitcounts, overwrite_hitcounts=true,kwargs...); #officially add hit positions from the new randomer
    end
    genes_hit=length(filter(x->x>0,hitcounts));
    genes_missing=findall(x->x==0,hitcounts);
    total_hits=sum(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),",")
    uniformity_weight=0
    hits_weight=0
    pool_size=length(randomers);
    data=DataFrame(poolsize=pool_size, primer_length=randomerlen, uniformity_weight=uniformity_weight,hits_weight=hits_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,uniformity_score=dist_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_V2.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;outputfile="./S_mutans_Cumulative_NSR_Pool.csv" ,kwargs...)
    mRNAfile=string("./Genome Data/",species_name,"_Genes.csv")
    mRNAdata=CSV.read(mRNAfile,DataFrame);
    mRNAs=read_mRNA_list(mRNAdata);
    nrandomers=length(randomers);
    cumulative_hits=zeros(nrandomers)
    cumulative_dist_score=zeros(nrandomers)
    hitcounts=zeros(length(mRNAs))
    randomer_num=1:nrandomers;
    str_randomers=dnaarray2stringarray(randomers);
    cumulative_genes_hit=zeros(nrandomers);
    for i=1:nrandomers
        final_score,hitcounts,dist_score,uniformity_weight,total_hits_weight=mRNA_uniform_coverage(randomers[i],mRNAs; hitcounts=hitcounts, overwrite_hitcounts=true,kwargs...);
        cumulative_hits[i]=sum(hitcounts)
        cumulative_dist_score[i]=dist_score;
        cumulative_genes_hit[i]=sum(hitcounts.>0)
    end
    data=DataFrame(run=randomer_num,randomer=str_randomers,cumulative_genes_hit=cumulative_genes_hit,cumulative_hits=cumulative_hits,cumulative_dist_score=cumulative_dist_score)
    CSV.write(outputfile,data)
end

#data=run_NSR_RL(;nsims=100,poolsize=1)
#CSV.write("NSR_RL_Prelims_Genome_Test_4_18_21.csv",data)
#speciesdata=CSV.read("NSR_RL_Prelims_Genome_Test_4_18_21.csv",DataFrame)
#SMUrandomers=stringarray2dnaarray(split(speciesdata[3,10],","));
#countdata=gene_coverage_counts(SMUrandomers,"./E_coli_Genes_4_15_21.csv")
#CSV.write("Ecoli_Coverage_Counts_4_18_21.csv",countdata)

#species_list=CSV.read("./Genome Data/Species_List.csv")[:Species]
#analyze_NSR_Pool_gene_coverage.(species_list)

#species_list=CSV.read("./Genome Data/Species_List.csv")[:Species]
#analyze_all_genomes_NSR_Pools(species_list)
#randomers=stringarray2dnaarray(CSV.read("./Brute Force Pools V2 Compressed/S_mutans_NSR_Brute_Force_Compressed.csv",DataFrame)[:Sequence])
#countdata=gene_coverage_counts(randomers,"./Genome Data/S_mutans_Genes.csv")
#CSV.write("S_mutans_NSR_Comparison_2_gene_coverage_5-19-21.csv",countdata)
#data=CSV.read("S_mutans_Add_Best_Sim_Option_Test_2.csv",DataFrame)
#randomers=stringarray2dnaarray(split(data[2,10],","));
#analyze_cumulative_NSR_Pool(randomers,"S_mutans";outputfile="./S_mutans_Cumulative_NSR_Pool_Best_Sim_Option_6_8_21_Uniformity.csv")
#data=CSV.read("S_mutans_NSR_Num_Blocked_Sites_Test_6_8_21.csv",DataFrame)
#randomers=stringarray2dnaarray(split(data[1,10],","))
#randomers=filter(x -> degeneracy(x;uselog2=false)>0,randomers)
#randomers=decompress_pool(randomers)
#rcomps=reverse_complement.(randomers)
#blocking_sites=read_blocking_sites(CSV.read("./Genome Data/S_mutans_rRNA_tRNA.csv",DataFrame));
#sites=unique(vcat(blocking_sites,rcomps,randomers))
#benchmark_NSR_RL()


text:
module TestSetup

using JNeuron

include("import_test.jl")
include("neuron_test.jl")
include("extra_test.jl")
include("visualize_test.jl")

end

ext: jl

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
using Math
using Documenter

DocMeta.setdocmeta!(Math, :DocTestSetup, :(using Math); recursive=true)

makedocs(;
    modules=[Math],
    authors="liuxu <liuxu89@mail.ustc.edu.cn> and contributors",
    repo="https://github.com/liuxu89/Math/blob/{commit}{path}#{line}",
    sitename="Math",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://liuxu89.github.io/Math",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "积分" => "integral.md",
        "未分类" => "tmp.md",
    ],
)

deploydocs(;
    repo="github.com/liuxu89/Math",
    devbranch="main",
)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Stage-decomposition evaluation #
# ============================== #
function evaluate_decision(structure::StageDecompositionStructure, decision::AbstractVector)
    # Evalaute decision stage-wise
    cᵀx = _eval_first_stage(structure, decision)
    𝔼Q = _eval_second_stages(structure, decision, objective_sense(structure.first_stage))
    # Return evaluation result
    return cᵀx + 𝔼Q
end

function statistically_evaluate_decision(structure::StageDecompositionStructure, decision::AbstractVector)
    # Evalaute decision stage-wise
    cᵀx = _eval_first_stage(structure, decision)
    𝔼Q, σ² = _stat_eval_second_stages(structure, decision, objective_sense(structure.first_stage))
    return cᵀx + 𝔼Q, sqrt(σ²)
end

function _eval_first_stage(structure::StageDecompositionStructure, decision::AbstractVector)
    # Update decisions (checks handled by first-stage model)
    take_decisions!(structure.first_stage,
                    all_decision_variables(structure.first_stage, 1),
                    decision)
    # Optimize first_stage model
    optimize!(structure.first_stage)
    # Switch on return status
    status = termination_status(structure.first_stage)
    result = if status in AcceptableTermination
        result = objective_value(structure.first_stage)
    else
        result = if status == MOI.INFEASIBLE
            result = objective_sense(structure.first_stage) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            result = objective_sense(structure.first_stage) == MOI.MAX_SENSE ? Inf : -Inf
        else
            error("First-stage model could not be solved, returned status: $status")
        end
    end
    # Revert back to untaken decisions
    untake_decisions!(structure.first_stage, all_decision_variables(structure.first_stage, 1))
    # Return evaluation result
    return result
end

function _eval_second_stages(structure::StageDecompositionStructure{2,1,Tuple{SP}},
                             decision::AbstractVector,
                             sense::MOI.OptimizationSense) where SP <: ScenarioProblems
    update_known_decisions!(structure.decisions[2], decision)
    map(subprob -> update_known_decisions!(subprob), subproblems(structure))
    return outcome_mean(subproblems(structure), probability.(scenarios(structure)), sense)
end
function _eval_second_stages(structure::StageDecompositionStructure{2,1,Tuple{SP}},
                             decision::AbstractVector,
                             sense::MOI.OptimizationSense) where SP <: DistributedScenarioProblems
    Qs = Vector{Float64}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async Qs[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                decision,
                sense) do sp, d, x, sense
                    scenarioproblems = fetch(sp)
                    decisions = fetch(d)
                    num_scenarios(scenarioproblems) == 0 && return 0.0
                    update_known_decisions!(decisions, x)
                    map(subprob -> update_known_decisions!(subprob), subproblems(scenarioproblems))
                    return outcome_mean(subproblems(scenarioproblems),
                                        probability.(scenarios(scenarioproblems)),
                                        sense)
                end
        end
    end
    return sum(Qs)
end

function _stat_eval_second_stages(structure::StageDecompositionStructure{2,1,Tuple{SP}},
                                  decision::AbstractVector,
                                  sense::MOI.OptimizationSense) where SP <: ScenarioProblems
    update_known_decisions!(structure.decisions[2], decision)
    map(subprob -> update_known_decisions!(subprob), subproblems(structure))
    return welford(subproblems(structure), probability.(scenarios(structure)), sense)
end
function _stat_eval_second_stages(structure::StageDecompositionStructure{2,1,Tuple{SP}},
                                  decision::AbstractVector,
                                  sense::MOI.OptimizationSense) where SP <: DistributedScenarioProblems
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                decision,
                sense) do sp, d, x, sense
                    scenarioproblems = fetch(sp)
                    decisions = fetch(d)
                    num_scenarios(scenarioproblems) == 0 && return zero(eltype(x)), zero(eltype(x)), zero(eltype(x)), zero(Int)
                    update_known_decisions!(decisions, x)
                    map(subprob -> update_known_decisions!(subprob), subproblems(scenarioproblems))
                    return welford(subproblems(scenarioproblems), probability.(scenarios(scenarioproblems)), sense)
                end
        end
    end
    𝔼Q, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼Q, σ²
end


text:
using Test

using CUDA

# generate some data on the CPU
dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = a .+ b

# upload to the GPU
d_a = CuArray(a)
d_b = CuArray(b)

# compute using array abstractions
d_c = d_a .+ d_b
@test c ≈ Array(d_c)

# compute using a kernel
function vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    c[i] = a[i] + b[i]
    return
end
len = prod(dims)
@cuda threads=len vadd(d_a, d_b, d_c)
@test c ≈ Array(d_c)


text:
# Programs used to visualize output, profile code or perform debugging tasks, as opposed to unit testing
module Diagnostics

using ..OpticSim
using ..OpticSim: replprint
using ..OpticSim.Vis
using ..OpticSim.Optimizable
using LinearAlgebra
using StaticArrays

import Unitful
import Plots

function testbigfloat()
    λ = 550 * Unitful.u"nm"
    temp = 20 * Unitful.u"°C"

    TOLERANCE = 1e-8

    function printintersections(a::Type{T}) where {T<:Real}
        println("* Intersections for type $T *")
        r1 = OpticalRay(Vector{T}([0.0, 0.0, 1.0]), Vector{T}([0.0, 0.0, -1.0]), one(T), T(0.55))
        r2 = OpticalRay(Vector{T}([2.0, 2.0, 1.0]), Vector{T}([0.0, 0.0, -1.0]), one(T), T(0.55))
        r3 = OpticalRay(Vector{T}([5.0, 5.0, 1.0]), Vector{T}([0.0, 0.0, -1.0]), one(T), T(0.55))
        r4 = OpticalRay(Vector{T}([0.0, -5.0, 1.0]), Vector{T}([0.0, sind(5.0), -cosd(5.0)]), one(T), T(0.55))
        r5 = OpticalRay(Vector{T}([-5.0, -5.0, 1.0]), Vector{T}([sind(5.0), sind(2.0), -cosd(2.0) * cosd(5.0)]), one(T), T(0.55))

        a = TestData.doubleconvex(T, temperature = temp)
        for r in (r1, r2, r3, r4, r5)
            println(point(OpticSim.intersection(OpticSim.trace(a, r, test = true))))
        end
        println()
    end

    printintersections(Float64)
    printintersections(BigFloat)

    a = TestData.doubleconvex(BigFloat)
    r2 = OpticalRay(Vector{BigFloat}([2.0, 2.0, 1.0]), Vector{BigFloat}([0.0, 0.0, -1.0]), one(BigFloat), BigFloat(0.55))
    r3 = OpticalRay(Vector{BigFloat}([5.0, 5.0, 1.0]), Vector{BigFloat}([0.0, 0.0, -1.0]), one(BigFloat), BigFloat(0.55))
    r4 = OpticalRay(Vector{BigFloat}([0.0, -5.0, 1.0]), Vector{BigFloat}([0.0, sind(5.0), -cosd(5.0)]), one(BigFloat), BigFloat(0.55))
    r5 = OpticalRay(Vector{BigFloat}([-5.0, -5.0, 1.0]), Vector{BigFloat}([sind(5.0), sind(2.0), -cosd(2.0) * cosd(5.0)]), one(BigFloat), BigFloat(0.55))
    @assert isapprox(point(OpticSim.intersection(OpticSim.trace(a, r2, test = true))), [-0.06191521590711035, -0.06191521590711035, -67.8], atol = TOLERANCE)
    @assert isapprox(point(OpticSim.intersection(OpticSim.trace(a, r3, test = true))), [-0.2491105067897657, -0.2491105067897657, -67.8], atol = TOLERANCE)
    @assert isapprox(point(OpticSim.intersection(OpticSim.trace(a, r4, test = true))), [0.0, 5.639876913179362, -67.8], atol = TOLERANCE)
    @assert isapprox(point(OpticSim.intersection(OpticSim.trace(a, r5, test = true))), [5.75170097290395, 2.504152441922817, -67.8], atol = TOLERANCE)
end

function vistest(sys::OpticalSystem{Float64}; kwargs...)
    λ = 0.550
    r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ)
    r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ)
    raygen = RayListSource(r1, r2, r3, r4, r5)
    Vis.drawtracerays(sys, raygenerator = raygen, trackallrays = true, test = true; kwargs...)
    for (i, r) in enumerate(raygen)
        t = OpticSim.trace(sys, r, test = true)
        if t !== nothing
            println("=== RAY $i ===")
            println(OpticSim.point(t))
            println(OpticSim.pathlength(t))
            println()
        end
    end
end

function visualizerefraction()
    nₛ = SVector{3,Float64}([0.0, 0.0, 1.0])
    c = 1 / sqrt(2.0)
    r = Ray([0.0, c, c], [0.0, -c, -c])
    rray = OpticSim.refractedray(1.0, 1.5, nₛ, direction(r))
    Vis.draw(Ray([0.0, 0.0, 0.0], Array{Float64,1}(nₛ)), color = :black)
    Vis.draw!(r, color = :green)
    Vis.draw!(Ray([0.0, 0.0, 0.0], Array{Float64,1}(rray)), color = :red)
    Vis.draw!(Ray([0.0, 0.0, 0.0], Array{Float64,1}(-nₛ)), color = :yellow)

    temp = OpticSim.refractedray(1.0, 1.5, nₛ, direction(r))
    temp = [0.0, temp[2], -temp[3]]
    r = Ray([0.0, -temp[2], -temp[3]], temp)
    rray = OpticSim.refractedray(1.5, 1.0, nₛ, direction(r))
    Vis.draw!(Ray([0.0, 0.0, 0.0], Array{Float64,1}(nₛ)), color = :black)
    Vis.draw!(r, color = :green)
    Vis.draw!(Ray([0.0, 0.0, 0.0], Array{Float64,1}(rray)), color = :red)
    Vis.draw!(Ray([0.0, 0.0, 0.0], Array{Float64,1}(-nₛ)), color = :yellow)
end

function plotreflectedvsrefractedpower()
    lens = OpticSim.SphericalLens(OpticSim.GlassCat.SCHOTT.BAK50, 0.0, Inf64, Inf64, 5.0, 10.0)
    reflectpow = Array{Float64,1}(undef, 0)
    refractpow = Array{Float64,1}(undef, 0)
    green = 500 * Unitful.u"nm"
    glass = OpticSim.GlassCat.SCHOTT.BAK50

    incidentindex = OpticSim.GlassCat.index(glass, green)
    transmittedindex = OpticSim.GlassCat.index(OpticSim.GlassCat.Air, green)

    for θ in 0.0:0.01:(π / 2.0)
        origin = [0.0, tan(θ), 1.0]
        dir = -origin
        r = Ray(origin, dir)

        intsct = OpticSim.closestintersection(OpticSim.surfaceintersection(lens(), r))

        nml = SVector{3,Float64}(0.0, 0.0, 1.0)

        rdir = direction(r)
        (nᵢ, nₜ) = OpticSim.mᵢandmₜ(incidentindex, transmittedindex, nml, r)
        reflected = OpticSim.reflectedray(nml, rdir)
        refracted = OpticSim.refractedray(nᵢ, nₜ, nml, rdir)

        (sinθᵢ, sinθₜ) = OpticSim.snell(nml, direction(r), nᵢ, nₜ)
        (powᵣ, powₜ) = OpticSim.fresnel(nᵢ, nₜ, sinθᵢ, sinθₜ)
        push!(reflectpow, powᵣ)
        push!(refractpow, powₜ)
    end

    Plots.plot((0.0:0.01:(π / 2.0)), [reflectpow, refractpow], label = ["reflected" "refracted"])
end

### OPTIMIZATION TESTING

using FiniteDifferences
using ForwardDiff
using DataFrames
using Optim
using JuMP
using Ipopt
using Zygote
using NLopt

doubleconvexprescription() = DataFrame(Surface = [:Object, 1, 2, :Image], Radius = [(Inf64), 60.0, -60.0, (Inf64)], Thickness = [(Inf64), (10.0), (77.8), missing], Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing], SemiDiameter = [(Inf64), (9.0), (9.0), (15.0)])

function doubleconvex(a::AbstractVector{T}; detpix::Int = 100) where {T<:Real}
    frontradius = a[1]
    rearradius = a[2]
    #! format: off
    AxisymmetricOpticalSystem{T}(DataFrame(
        Surface = [:Object, 1, 2, :Image],
        Radius = [T(Inf64), frontradius, rearradius, T(Inf64)],
        Conic = [missing, -1.0, 1.0, missing],
        Thickness = [T(Inf64), T(10.0), T(77.8), missing],
        Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, missing],
        SemiDiameter = [T(Inf64), T(9.0), T(9.0), T(15.0)],
    ), detpix, detpix, T, temperature = OpticSim.GlassCat.TEMP_REF_UNITFUL, pressure = OpticSim.GlassCat.PRESSURE_REF)
    #! format: on
end

function RMS_spot_size(a::AbstractVector{T}, b::AxisymmetricOpticalSystem{T}, samples::Int = 3) where {T}
    # RMSE spot size
    lens = Optimizable.updateoptimizationvariables(b, a)
    field = HexapolarField(lens, collimated = true, samples = samples)
    error = zero(T)
    hits = 0
    for r in field
        traceres = OpticSim.trace(lens, r, test = true)
        if traceres !== nothing
            hitpoint = point(traceres)
            if abs(hitpoint[1]) > eps(T) && abs(hitpoint[2]) > eps(T)
                dist_to_axis = hitpoint[1]^2 + hitpoint[2]^2
                error += dist_to_axis
            end
            hits += 1
        end
    end
    if hits > 0
        error = sqrt(error / hits)
    end
    return error
end

function testfinitedifferences()
    lens = Examples.doubleconvex(60.0, -60.0)
    start = Optimizable.optimizationvariables(lens)

    u = 60.0
    println("time to compute objective function $(@time RMS_spot_size([u,-60.0],lens))")
    fdm = central_fdm(3, 1)
    f1(u) = RMS_spot_size([u, -60.0], lens)
    fu = 0.0
    for i in 1:100
        fu = fdm(f1, 60.0)
    end
    return fu
end
export testfinitedifferences

function testoptimization(; lens = Examples.cooketriplet(), constrained = false, algo = constrained ? IPNewton() : LBFGS(), chunk_size = 1, samples = 3)
    @info "NaN safe: $(ForwardDiff.NANSAFE_MODE_ENABLED)"

    start, lower, upper = Optimizable.optimizationvariables(lens)
    optimobjective = (arg) -> RMS_spot_size(arg, lens, samples)

    if chunk_size == 0
        gcfg = ForwardDiff.GradientConfig(optimobjective, start)
        hcfg = ForwardDiff.HessianConfig(optimobjective, start)
    else
        gcfg = ForwardDiff.GradientConfig(optimobjective, start, ForwardDiff.Chunk{chunk_size}())
        hcfg = ForwardDiff.HessianConfig(optimobjective, start, ForwardDiff.Chunk{chunk_size}())
        @info "$(length(start)) optimization variables, chunk size = $chunk_size"
    end

    g! = (G, x) -> ForwardDiff.gradient!(G, optimobjective, x, gcfg)
    h! = (H, x) -> ForwardDiff.hessian!(H, optimobjective, x, hcfg)

    @info "Starting objective function $(optimobjective(start))"
    if constrained && !all(isinf.(abs.(lower))) && !all(isinf.(upper))
        if !(algo isa Optim.ConstrainedOptimizer)
            @error "This algorithm can't handle constraints, use IPNewton"
            return
        end
        @info "Constraints: $lower, $upper"
        df = TwiceDifferentiable(optimobjective, g!, h!, start)
        dfc = TwiceDifferentiableConstraints(lower, upper)
        res = optimize(df, dfc, start, algo, Optim.Options(show_trace = true, iterations = 100, allow_f_increases = true))
    else
        if constrained
            @warn "No constraints to apply, using unconstrained optimization"
        end
        res = optimize(optimobjective, g!, h!, start, algo, Optim.Options(show_trace = true, iterations = 100, allow_f_increases = true))
    end

    # println("== Computing gradient")
    # println("First call:")
    # @time g!(zeros(length(start)), start)
    # println("Second call:")
    # @time g!(zeros(length(start)), start)

    # println("== Computing hessian")
    # println("First call:")
    # @time h!(zeros(length(start), length(start)), start)
    # println("Second call:")
    # @time h!(zeros(length(start), length(start)), start)

    # println()

    final = Optim.minimizer(res)
    newlens = Optimizable.updateoptimizationvariables(lens, final)
    @info "Result: $final"

    field = HexapolarField(newlens, collimated = true, samples = samples)
    Vis.drawtraceimage(newlens, raygenerator = field)
    Vis.drawtracerays(newlens, raygenerator = field, trackallrays = true, test = true)
end

function testnlopt()
    lens = Examples.doubleconvex(60.0, -60.0)
    # lens = Examples.doubleconvex()
    # lens = Examples.telephoto(6,.5)
    # Vis.drawtraceimage(lens)
    start = OpticSim.optimizationvariables(lens)

    optimobjective = (arg) -> RMS_spot_size(arg, lens)
    println("starting objective function $(optimobjective(start))")

    model = Model(NLopt.optimizer)
    @variable(model, rad1, start = 60.0)
    @variable(model, rad2, start = -60.0)
    register(model, :f, 2, f, autodiff = true)
    @NLconstraint(model, 30.0 <= rad1 <= 85.0)
    @NLconstraint(model, -100.0 <= rad1 <= -55.0)
    @NLobjective(model, Min, f(rad1, rad2))
    optimize!(model)

    println("rad1 $(value(rad1)) rad2 $(value(rad2))")

end

function testJuMP()
    function innerfunc(radius1::T, radius2::T) where {T<:Real}

        lens = Examples.doubleconvex(radius1, radius2)

        field = HexapolarField(lens, collimated = true, samples = 10)
        error = zero(T)
        hits = 0
        for r in field
            traceres = OpticSim.trace(lens, r, test = true)

            if !(nothing === traceres)
                hitpoint = point(traceres)
                if abs(hitpoint[1]) > eps(T) && abs(hitpoint[2]) > eps(T)
                    dist_to_axis = sqrt(hitpoint[1]^2 + hitpoint[2]^2)
                    error += dist_to_axis
                end
                hits += 1
            end
        end

        # if isa(radius1,ForwardDiff.Dual)
        #     println("radius1 $(realpart(radius1)) radius2 $(realpart(radius2))")
        # end

        error /= hits
        # println("error in my code $error")

        # Vis.drawtracerays(doubleconvex(radiu1,radius2),trackallrays = true, test = true)

        return error
    end

    function f(radius1::T, radius2::T) where {T<:Real}
        innerfunc(radius1, radius2)
        # try
        #     innerfunc(radius1, radius2)
        # catch err
        #     return zero(T)
        # end
    end

    # Vis.drawtracerays(doubleconvex(60.0,-60.0), trackallrays = true, test = true)

    model = Model(Ipopt.Optimizer)
    set_time_limit_sec(model, 10.0)
    register(model, :f, 2, f, autodiff = true)
    @variable(model, 10.0 <= rad1 <= 85.0, start = 60.0)
    @variable(model, -150 <= rad2 <= 20.0, start = -10.0)
    @NLobjective(model, Min, f(rad1, rad2))
    optimize!(model)

    println("rad1 $(value(rad1)) rad2 $(value(rad2))")
    Vis.drawtracerays(Examples.doubleconvex(value(rad1), value(rad2)), trackallrays = true, test = true)
end

function simpletest()
    lens = Examples.cooketriplet()
    # lens = Examples.doubleconvex(60.0,-60.0)
    testgenericoptimization(lens, RMS_spot_size)
end

function RMS_spot_size(lens, x::T...) where {T}
    lens = Optimizable.updateoptimizationvariables(lens, collect(x))
    field = HexapolarField(lens, collimated = true, samples = 10)
    error = zero(T)
    hits = 0
    for r in field
        traceres = OpticSim.trace(lens, r, test = true)

        if !(nothing === traceres)
            hitpoint = point(traceres)
            if abs(hitpoint[1]) > eps(T) && abs(hitpoint[2]) > eps(T)
                dist_to_axis = sqrt(hitpoint[1]^2 + hitpoint[2]^2)
                error += dist_to_axis
            end
            hits += 1
        end
    end

    error /= hits

    return error
end

function testZygoteGradient(lens::S, objective::Function) where {S<:OpticalSystem}
    vars = Optimizable.optimizationvariables(lens)
    gradient = objective'(vars)
    println(gradient)
end

function testgenericoptimization(lens::S, objective::Function) where {S<:OpticalSystem}
    vars = Optimizable.optimizationvariables(lens)
    function newobjective(a)
        objective(lens, a...)
    end

    numvars = length(vars)
    model = Model(Ipopt.Optimizer)
    set_time_limit_sec(model, 10.0)
    register(model, :newobjective, numvars, newobjective, autodiff = true)
    # @variable(model, bounds[i][1] <= x[i=1:numvars] <= bounds[i][2],start = vars)
    @variable(model, x[i = 1:numvars], start = vars[i])

    @NLobjective(model, Min, newobjective(x...))
    # d = NLPEvaluator(model)
    # MOI.initialize(d, [:Grad])
    # println("objective $(MOI.eval_objective(d, vars))")

    # ∇f = zeros(length(vars))
    # println(vars)
    # MOI.eval_objective_gradient(d, ∇f, vars)
    # println("gradient $∇f")
    # println("starting forward diff")
    # time = @time ForwardDiff.gradient(newobjective,vars)
    # println("forward diff time1 $time")
    # time = @time ForwardDiff.gradient(newobjective,vars)
    # println("forward diff time2 $time")
    println("starting reverse diff")
    time = @time ReverseDiff.gradient(newobjective, vars)
    println("reverse diff time1 $time")
    time = @time ReverseDiff.gradient(newobjective, vars)
    println("reverse diff time2 $time")

    println("ending forward diff")
    # fdm = central_fdm(2,1)
    # replprint( @benchmark(FiniteDifferences.grad($fdm,$newobjective,$vars...)))
    # println("finite difference gradient $(FiniteDifferences.grad(fdm,newobjective,vars...))")
    # optimize!(model)
    lens = Optimizable.updateoptimizationvariables(lens, value.(x))

    Vis.drawtracerays(lens, trackallrays = true, test = true)
end

function testIpopt()

    model = Model(Ipopt.Optimizer)

    @variable(model, x1, start = 1)
    @variable(model, x2, start = 5)
    @variable(model, x3, start = 5)
    @variable(model, x4, start = 1)

    f(x1, x2, x3, x4) = x1 * x4 * (x1 + x2 + x3) + x3

    register(model, :f, 4, f, autodiff = true)
    @NLobjective(model, Min, f(x1, x2, x3, x4))
    @NLconstraint(model, x1 * x2 * x3 * x4 >= 25)
    @NLconstraint(model, x1^2 + x2^2 + x3^2 + x4^2 == 40)

    optimize!(model)

    println("x1 $(value(x1)) x2 $(value(x2)) x3 $(value(x3)) x4 $(value(x4))")
end

function testtypesignature()
    println("here")
    a = OpticSim.Sphere(1.0)
    b = csgintersection(csgunion(a, a), csgunion(a, a))
end

function testoptimizationvariables()
    lens = Examples.cooketriplet()
    vars = OpticSim.optimizationvariables(lens)
    vars .= [Float64(i) for i in 1:length(vars)]
    newlens = OpticSim.updateoptimizationvariables(lens, vars)
    println("original lens")
    show(lens)
    println("updated lens")
    show(newlens)
end

end #module
export Diagnostics


text:
module ExtrudedArrays

using StaticArrays
using Base: RefValue
using Base.Broadcast: Broadcasted, Extruded
using Base.Broadcast: newindexer
using Swizzles.WrapperArrays
using Swizzles.ShallowArrays
using Swizzles.ArrayifiedArrays
using Swizzles.Virtuals
using Swizzles.Properties
using Swizzles: combinetuple
using Swizzles: maptuple
using LinearAlgebra

using Base: tail

export ExtrudedArray
export keeps, kept
export Keep, Extrude
export keep, extrude
export lift_keeps, freeze_extrudes

struct Keep end
struct Extrude end
const keep = Keep()
const extrude = Extrude()

Base.:|(::Keep, ::Keep) = keep
Base.:|(::Keep, ::Extrude) = keep
Base.:|(::Keep, ::Bool) = keep

Base.:|(::Extrude, ::Keep) = keep
Base.:|(::Extrude, ::Extrude) = extrude
Base.:|(::Extrude, k::Bool) = k

Base.:|(::Bool, ::Keep) = keep
Base.:|(k::Bool, ::Extrude) = k

kept(::Keep) = true
kept(::Extrude) = false
kept(k::Bool) = k

struct ExtrudedArray{T, N, Arg<:AbstractArray{T, N}, Keeps<:Tuple{Vararg{Any, N}}} <: ShallowArray{T, N, Arg}
    arg::Arg
    keeps::Keeps
    function ExtrudedArray{T, N, Arg, Keeps}(arg::Arg, keeps::Keeps) where {T, N, Arg, Keeps}
        @assert keeps isa Tuple{Vararg{Any, N}}
        @assert all(ntuple(n -> kept(keeps[n]) || size(arg, n) == 1, N)) "$(size(arg)) $(keeps)"
        return new{T, N, Arg, Keeps}(arg, keeps)
    end
end

function ExtrudedArray(arg)
    return ExtrudedArray(arg, maptuple(k-> kept(k) ? keep : extrude, keeps(arg)...))
end
function ExtrudedArray(arg::ExtrudedArray)
    return ExtrudedArray(arg.arg, maptuple(k-> kept(k) ? keep : extrude, keeps(arg)...))
end
function ExtrudedArray(arg, keeps)
    return ExtrudedArray{eltype(arg), ndims(arg), typeof(arg), typeof(keeps)}(arg, keeps)
end
function ExtrudedArray(arg::ExtrudedArray, keeps)
    return ExtrudedArray(arg.arg, keeps)
end

@inline Base.parent(arr::ExtrudedArray) = arr.arg
@inline WrapperArrays.adopt(arr::ExtrudedArray, arg::AbstractArray) = ExtrudedArray{eltype(arg), ndims(arr), typeof(arg), typeof(arr.keeps)}(arg, arr.keeps)

#keeps returns a tuple where each element of the tuple specifies whether the
#corresponding dimension is intended to have size 1, possibly using the traits
#keep and extrude (instances of the singleton types `keep` and `extrude`).
keeps(x) = ntuple(ndims(x)) do n
    return size(x, n) != 1
end
keeps(::Tuple{}) = (keep,)
keeps(::Tuple{Any}) = (extrude,)
keeps(::Tuple) = (keep,)
function keeps(arr::Adjoint)
    arg_keeps = keeps(parent(arr))
    if length(arg_keeps) == 0
        return (extrude, extrude)
    elseif length(arg_keeps) == 1
        return (extrude, arg_keeps[1])
    elseif length(arg_keeps) == 2
        return (arg_keeps[2], arg_keeps[1])
    else
        throw(ArgumentError("TODO"))
    end
end
function keeps(arr::Transpose)
    arg_keeps = keeps(parent(arr))
    if length(arg_keeps) == 0
        return (extrude, extrude)
    elseif length(arg_keeps) == 1
        return (extrude, arg_keeps[1])
    elseif length(arg_keeps) == 2
        return (arg_keeps[2], arg_keeps[1])
    else
        throw(ArgumentError("TODO"))
    end
end
keeps(arr::ExtrudedArray) = arr.keeps
keeps(ext::Extruded) = ext.keeps
keeps(bc::Broadcasted) = combinetuple(|, map(keeps, bc.args)...)
keeps(arr::ArrayifiedArray) = keeps(arr.arg)
keeps(arr::StaticArray) = ntuple(ndims(arr)) do n
    return size(arr, n) == 1 ? keep : extrude
end

import Base.Broadcast

@inline Base.Broadcast.extrude(x::Extruded) = x
Base.@propagate_inbounds Base.Broadcast.newindex(arg::ExtrudedArray, I::CartesianIndex) = CartesianIndex(Base.Broadcast._newindex(I.I, keeps(arg), maptuple(first, axes(arg)...)))
Base.@propagate_inbounds Base.Broadcast.newindex(arg::ExtrudedArray, I::Integer) = CartesianIndex(Base.Broadcast._newindex((I,), keeps(arg), maptuple(first, axes(arg)...)))
@inline Base.Broadcast.newindexer(A::ExtrudedArray) = (keeps(A), maptuple(first, axes(A)...))
@inline Base.Broadcast.newindex(i::Integer, ::Tuple{Keep}, idefault) = i
@inline Base.Broadcast.newindex(i::Integer, ::Tuple{Extrude}, idefault) = idefault[i]
@inline Base.Broadcast._newindex(I, keeps::Tuple{Keep, Vararg{Any}}, Idefault) =
    (I[1], Base.Broadcast._newindex(tail(I), tail(keeps), tail(Idefault))...)
@inline Base.Broadcast._newindex(I, keeps::Tuple{Extrude, Vararg{Any}}, Idefault) =
    (Idefault[1], Base.Broadcast._newindex(tail(I), tail(keeps), tail(Idefault))...)

static_extrude(x) = Extruded(x, maptuple(k-> kept(k) ? keep : extrude, keeps(x)...), maptuple(first, axes(x)...))

#Add stable Extrudes to all internal broadcast expressions so that the broadcast_getindex does not need dynamic checks.
_freeze_extrudes(bc::Broadcasted{Style}) where {Style} = Broadcasted{Style}(bc.f, maptuple(_freeze_extrudes, bc.args...), bc.axes)
_freeze_extrudes(ext::Extruded) = static_extrude(freeze_extrudes(ext.x)) #FIXME avoid recomputation in redundant passes
_freeze_extrudes(x) = static_extrude(freeze_extrudes(x))
freeze_extrudes(bc::Broadcasted) = _freeze_extrudes(bc)
freeze_extrudes(x) = x
function freeze_extrudes(arr::AbstractArray)
    if iswrapper(arr)
        adopt(freeze_extrudes(parent(arr)), arr)
    else
        arr
    end
end
function freeze_extrudes(arr::ArrayifiedArray)
    return arrayify(freeze_extrudes(arr.arg))
end

#Add ExtrudedArrays to a broadcast/lazyarray expression so that the keeps of every recognizeable node can be determined from the types.
lift_keeps(bc::Broadcasted{Style}) where {Style} = Broadcasted{Style}(bc.f, maptuple(lift_keeps, bc.args...), bc.axes)
lift_keeps(ext::Extruded) = static_extrude(lift_keeps(ext.x))
lift_keeps(arr::ExtrudedArray) = ExtrudedArray(lift_keeps(ext.x))
lift_keeps(arr::ArrayifiedArray) = arrayify(lift_keeps(arr.arg))
function lift_keeps(x)
    if ndims(x) == 0
        return x
    else
        return ExtrudedArray(arrayify(x))
    end
end
lift_keeps(tu::Tuple) = tu
function lift_keeps(arr::Adjoint{T}) where {T}
    arg = lift_keeps(parent(arr))
    return Adjoint{T, typeof(arg)}(arg)
end
function lift_keeps(arr::Transpose{T}) where {T}
    arg = lift_keeps(parent(arr))
    return Transpose{T, typeof(arg)}(arg)
end
function lift_keeps(arr::AbstractArray)
    if ndims(arr) == 0
        return arr
    else
        ExtrudedArray(arr)
    end
end

function Virtuals.virtual(root, ::Type{Keep})
    return keep
end
function Virtuals.virtual(root, ::Type{Extrude})
    return extrude
end
function Virtuals.virtualize(root, ::Type{ExtrudedArray{T, N, Arg, Keeps}}) where {T, N, Arg, Keeps}
    arg = virtualize(:($root.arg), Arg)
    keeps = virtualize(:($root.keeps), Keeps)
    return ExtrudedArray{T, N, typeof(arg), typeof(keeps)}(arg, keeps)
end

Base.:|(::Keep, ::Virtual{Bool}) = keep
Base.:|(::Virtual{Bool}, ::Keep) = keep
Base.:|(::Extrude, y::Virtual{Bool}) = y
Base.:|(x::Virtual{Bool}, ::Extrude) = x
Base.:|(x::Virtual{Bool}, ::Virtual{Bool}) = Virtual{Bool}(:($(x.ex) | $(y.ex)))

keeps(::Virtual{Tuple{}}) = (keep,)
keeps(::Virtual{Tuple{Any}}) = (extrude,)
keeps(::Virtual{<:Tuple}) = (keep,)

end


text:
quantile_exp(lam::Number, rand1::Number)::Number = -log(1 - rand1) / lam;

function gausscopulagen2(t::Int, Σ::Matrix{Float64}, mcBaseData::AbstractMonteCarloConfiguration)
    z = rand(mcBaseData.rng, MvNormal(Σ), t)
    for i = 1:size(Σ, 2)
        d = Normal(0, sqrt(Σ[i, i]))
        @views z[i, :] .= cdf.(d, z[i, :])
    end
    return z
end

function gausscopulagen2!(z, Σ::Matrix{Float64}, mcBaseData::AbstractMonteCarloConfiguration)
    rand!(mcBaseData.rng, MvNormal(Σ), z)
    for i = 1:size(Σ, 2)
        d = Normal(0, sqrt(Σ[i, i]))
        @views z[i, :] .= cdf.(d, z[i, :])
    end
    return z
end


text:
highlighter_exists(config::CthulhuConfig) =
    Sys.which(config.highlighter.exec[1]) !== nothing

__init__() = CONFIG.enable_highlighter = highlighter_exists(CONFIG)

function highlight(io, x, lexer, config::CthulhuConfig)
    config.enable_highlighter || return print(io, x)
    if !highlighter_exists(config)
        @warn "Highlighter command $(config.highlighter.exec[1]) does not exist."
        return print(io, x)
    end
    cmd = `$(config.highlighter) $lexer`
    open(pipeline(cmd; stdout=io, stderr=stderr), "w") do io
        print(io, x)
    end
end

function cthulhu_llvm(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    @static if VERSION >= v"1.5.0-DEV.393"
        dump = InteractiveUtils._dump_function_linfo_llvm(
            mi, params.world,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            #=dump_module=# false,
            optimize, debuginfo ? :source : :none, Base.CodegenParams())
    else
        dump = InteractiveUtils._dump_function_linfo(
            mi, params.world, #=native=# false,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            #=dump_module=# false, #=syntax=# config.asm_syntax,
            optimize, debuginfo ? :source : :none)
    end
    highlight(io, dump, "llvm", config)
end

function cthulhu_native(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    @static if VERSION >= v"1.5.0-DEV.393"
        dump = InteractiveUtils._dump_function_linfo_native(
            mi, params.world,
            #=wrapper=# false, #=syntax=# config.asm_syntax,
            debuginfo ? :source : :none)
    else
        dump = InteractiveUtils._dump_function_linfo(
            mi, params.world, #=native=# true,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            #=dump_module=# false, #=syntax=# config.asm_syntax,
            optimize, debuginfo ? :source : :none)
    end
    highlight(io, dump, "asm", config)
end

function cthulhu_ast(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    ast = definition(Expr, meth)
    if ast!==nothing
        if !config.pretty_ast
            dump(io, ast; maxdepth=typemax(Int))
        else
            show(io, ast)
            # Meta.show_sexpr(io, ast)
            # Could even highlight the above as some kind-of LISP
        end
    else
        @info "Could not retrieve AST. AST display requires Revise.jl to be loaded." meth
    end
end

function cthulhu_source(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    src, line = definition(String, meth)
    highlight(io, src, "julia", config)
end

cthulhu_warntype(args...) = cthulhu_warntype(stdout, args...)
function cthulhu_warntype(io::IO, src, rettype, debuginfo)
    if VERSION < v"1.1.0-DEV.762"
    elseif VERSION < v"1.2.0-DEV.229"
        lineprinter = Base.IRShow.debuginfo[debuginfo]
    else
        debuginfo = Base.IRShow.debuginfo(debuginfo)
        lineprinter = Base.IRShow.__debuginfo[debuginfo]
    end

    lambda_io::IOContext = stdout
    if src.slotnames !== nothing
        lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES =>  Base.sourceinfo_slotnames(src))
    end
    print(io, "Body")
    InteractiveUtils.warntype_type_printer(io, rettype, true)
    println(io)
    if VERSION < v"1.1.0-DEV.762"
        Base.IRShow.show_ir(lambda_io, src, InteractiveUtils.warntype_type_printer)
    else
        Base.IRShow.show_ir(lambda_io, src, lineprinter(src), InteractiveUtils.warntype_type_printer)
    end
    return nothing
end

function cthulu_typed(io::IO, debuginfo_key, CI, rettype, mi, iswarn)
    println()
    println("│ ─ $(string(Callsite(-1, MICallInfo(mi, rettype))))")

    if iswarn
        cthulhu_warntype(stdout, CI, rettype, debuginfo_key)
    elseif VERSION >= v"1.1.0-DEV.762"
        show(stdout, CI, debuginfo = debuginfo_key)
    else
        display(CI=>rt)
    end
    println()
end

# These are standard code views that don't need any special handling,
# This namedtuple maps toggle::Symbol to function
const codeviews = (;
    llvm=cthulhu_llvm,
    native=cthulhu_native,
    ast=cthulhu_ast,
    source=cthulhu_source,
)


text:
# Variational                       #
# --------------------------------- #
@rule MvNormalMeanPrecision(:Λ, Marginalisation) (q_out::Any, q_μ::Any) = begin
    m_out, v_out   = mean_cov(q_out)
    m_mean, v_mean = mean_cov(q_μ)

    df = ndims(q_μ) + 2
    S  = Matrix(Hermitian(cholinv(v_mean + v_out + (m_mean - m_out) * (m_mean - m_out)')))

    return Wishart(df, S)
end

@rule MvNormalMeanPrecision(:Λ, Marginalisation) (q_out_μ::Any, ) = begin
    m_out_μ, v_out_μ   = mean_cov(q_out_μ)

    d = div(ndims(q_out_μ), 2)

    mdiff = @views m_out_μ[1:d] - m_out_μ[d+1:end]
    vdiff = @views v_out_μ[1:d,1:d] - v_out_μ[1:d,d+1:end] - v_out_μ[d+1:end, 1:d] + v_out_μ[d+1:end,d+1:end]
    S     = Matrix(Hermitian(cholinv(vdiff + mdiff * mdiff')))

    return Wishart(d + 2, S)
end


text:
using ConstrainedDynamics
using ConstrainedDynamics: LVᵀmat, Rotational, Translational, ∂2g∂posaa, ∂2g∂posab, ∂2g∂posba, ∂2g∂posbb, g
using ForwardDiff
using Rotations
using Rotations: params
using LinearAlgebra

function transtest()
    xa = srand(3)
    qa = rand(UnitQuaternion)
    xb = srand(3)
    qb = rand(UnitQuaternion)
    z = [xa;params(qa);xb;params(qb)]

    axis = srand(3)
    vert = (srand(3),srand(3))

    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))
    trans = Translational{Float64,2}(body1, body2; axis = axis, p1 = vert[1], p2 = vert[2])[1]

    XXaa, XQaa, QXaa, QQaa = ∂2g∂posaa(trans, xa, qa, xb, qb)
    XXab, XQab, QXab, QQab = ∂2g∂posab(trans, xa, qa, xb, qb)
    XXba, XQba, QXba, QQba = ∂2g∂posba(trans, xa, qa, xb, qb)
    XXbb, XQbb, QXbb, QQbb = ∂2g∂posbb(trans, xa, qa, xb, qb)

    G = [
        XXaa XQaa XXab XQab
        QXaa QQaa QXab QQab
        XXba XQba XXbb XQbb
        QXba QQba QXbb QQbb
    ]

    function firstderivative(z)
        E = diagm(ones(3))
        Z34 = zeros(3,4)
        Z43 = zeros(4,3)
        postmultmat = [
            E 0*I 0*I 0*I
            Z43 LVᵀmat(UnitQuaternion(z[4:7]...,false)) Z43 Z43
            0*I 0*I E 0*I
            Z43 Z43 Z43 LVᵀmat(UnitQuaternion(z[11:14]...,false))
        ]
        ForwardDiff.jacobian(x->g(trans,x[1:3],UnitQuaternion(x[4:7]...,false),x[8:10],UnitQuaternion(x[11:14]...,false)),z)*postmultmat
    end
    function secondderivative(z)
        ForwardDiff.jacobian(x->vec(firstderivative(x)),z)
    end
    @test isapprox(norm(G-secondderivative(z)), 0.0; atol = 1e-8)
end

function rottest()
    xa = srand(3)
    qa = rand(UnitQuaternion)
    xb = srand(3)
    qb = rand(UnitQuaternion)
    z = [xa;params(qa);xb;params(qb)]

    axis = srand(3)
    qoffset = rand(UnitQuaternion)

    body1 = Body(Box(1., 1., 1., 1.))
    body2 = Body(Box(1., 1., 1., 1.))
    rot = Rotational{Float64,2}(body1, body2; axis = axis, qoffset = qoffset)[1]

    XXaa, XQaa, QXaa, QQaa = ∂2g∂posaa(rot, xa, qa, xb, qb)
    XXab, XQab, QXab, QQab = ∂2g∂posab(rot, xa, qa, xb, qb)
    XXba, XQba, QXba, QQba = ∂2g∂posba(rot, xa, qa, xb, qb)
    XXbb, XQbb, QXbb, QQbb = ∂2g∂posbb(rot, xa, qa, xb, qb)

    G = [
        XXaa XQaa XXab XQab
        QXaa QQaa QXab QQab
        XXba XQba XXbb XQbb
        QXba QQba QXbb QQbb
    ]

    function firstderivative(z)
        E = diagm(ones(3))
        Z34 = zeros(3,4)
        Z43 = zeros(4,3)
        postmultmat = [
            E 0*I 0*I 0*I
            Z43 LVᵀmat(UnitQuaternion(z[4:7]...,false)) Z43 Z43
            0*I 0*I E 0*I
            Z43 Z43 Z43 LVᵀmat(UnitQuaternion(z[11:14]...,false))
        ]
        ForwardDiff.jacobian(x->g(rot,x[1:3],UnitQuaternion(x[4:7]...,false),x[8:10],UnitQuaternion(x[11:14]...,false)),z)*postmultmat
    end
    function secondderivative(z)
        ForwardDiff.jacobian(x->vec(firstderivative(x)),z)
    end
    @test isapprox(norm(G-secondderivative(z)), 0.0; atol = 1e-8)
end

for i=1:10
    transtest()
    rottest()
end


text:
using PyPlot, DataFrames, Query
pygui(true)
csvraw = joinpath("data", "Statistikdatabasen_2017-12-23 16_43_51.csv")
csvform = joinpath("data", "influi0916.csv")

if !(ispath(csvform))
	run(pipeline(`./sos.awk`, stdin=csvraw, stdout=csvform))
end

influi0916 = readtable(csvform)
influi0916a = @from i in influi0916 begin
		@group i by @NT(Ålder=i.Ålder, Kön=i.Kön) into g
		@select {Ålder=g.key.Ålder,
			Kön=g.key.Kön,
			Incidens=mean(g..Incidens)}
		@collect DataFrame
	end

for sex in ["Kvinnor"; "Män"]
	influi0916as = influi0916a[influi0916a[:Kön].==sex, :]
	plot(influi0916as[:Ålder], log(influi0916as[:Incidens]./10^5),
		"-o", label = "$sex")
end
grid(1)
xlabel("Ålder")
ylabel("log(incidens)")
title("Sjukhusvårdade influensa/lunginflammation Sverige 2009\u201316")
legend()

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# For generating symbols that work as dataframe column names; symbols
# generated by gensym() don't work for this.
global _rvnum = 0

function _make_rvname(name)
    global _rvnum += 1
    return Symbol("$(name)!$(_rvnum)")
end

function _make_dims(args)
    dims = Vector{Any}()
    for arg in args
        # println("Arg: $arg")

        if @capture(arg, i_Int)  # scalar (must be integer)
            # println("arg is Int")
            dim = i

        elseif @capture(arg, first_Int:last_)   # last can be an int or 'end', which is converted to 0
            # println("arg is range")
            last = last == :end ? 0 : last
            dim = :($first:$last)

        elseif @capture(arg, first_Int:step_:last_)
            # println("arg is step range")
            last = last == :end ? 0 : last
            dim = :($first:$step:$last)

        elseif @capture(arg, s_Symbol)
            if arg == :(:)
                # println("arg is Colon")
                dim = Colon()
            else
                # println("arg is Symbol")
                dim = :($(QuoteNode(s)))
            end

        elseif @capture(arg, s_String)
            dim = s

        elseif isa(arg, Expr) && arg.head == :tuple  # tuple of Strings/Symbols (@capture didn't work...)
            argtype = typeof(arg.args[1])            # ensure all are same type as first element
            if ! isempty(filter(s -> typeof(s) != argtype, arg.args))
                error("A parameter dimension tuple must all String or all Symbol (got $arg)")
            end
            dim = :(convert(Vector{$argtype}, $(arg.args)))

        else
            error("Unrecognized stochastic parameter specification: $arg")
        end
        push!(dims, dim)
    end
    # println("dims = $dims")
    return dims
end

macro defmcs(expr)
    let # to make vars local to each macro invocation
        local _rvs        = []
        local _corrs      = []
        local _transforms = []
        local _saves      = []

        # distilled into a function since it's called from two branches below
        function saverv(rvname, distname, distargs)
            expr = :(RandomVariable($(QuoteNode(rvname)), $distname($(distargs...))))
            push!(_rvs, esc(expr))
        end

        @capture(expr, elements__)
        for elt in elements
            # Meta.show_sexpr(elt)
            # println("")
            # e.g.,  rv(name1) = Normal(10, 3)
            if @capture(elt, rv(rvname_) = distname_(distargs__))
                saverv(rvname, distname, distargs)

            elseif @capture(elt, save(vars__))
                for var in vars
                    # println("var: $var")
                    if @capture(var, comp_.datum_)
                        expr = :($(QuoteNode(comp)), $(QuoteNode(datum)))
                        push!(_saves, esc(expr))
                    else
                        error("Save arg spec must be of the form comp_name.datum_name; got ($var)")
                    end
                end

            # handle vector of distributions
            elseif @capture(elt, extvar_ = [items__])
                for pair in items
                    if (@capture(pair, [dims__] => distname_(distargs__)) ||
                        @capture(pair,     dim_ => distname_(distargs__)))

                        dims = _make_dims(dims === nothing ? [dim] : dims)

                        rvname = _make_rvname(extvar)
                        saverv(rvname, distname, distargs)

                        expr = :(TransformSpec($(QuoteNode(extvar)), :(=), $(QuoteNode(rvname)), [$(dims...)]))
                        push!(_transforms, esc(expr))
                    end
                end

            # e.g., name1:name2 = 0.7
            elseif @capture(elt, name1_:name2_ = value_)
                expr = :(CorrelationSpec($(QuoteNode(name1)), $(QuoteNode(name2)), $value))
                push!(_corrs, esc(expr))

            # e.g., ext_var5[2010:2050, :] *= name2
            # A bug in Macrotools prevents this shorter expression from working:
            # elseif @capture(elt, ((extvar_  = rvname_Symbol) |
            #                       (extvar_ += rvname_Symbol) |
            #                       (extvar_ *= rvname_Symbol) |
            #                       (extvar_  = distname_(distargs__)) |
            #                       (extvar_ += distname_(distargs__)) |
            #                       (extvar_ *= distname_(distargs__))))
            elseif (@capture(elt, extvar_  = rvname_Symbol) ||
                    @capture(elt, extvar_ += rvname_Symbol) ||
                    @capture(elt, extvar_ *= rvname_Symbol) ||
                    @capture(elt, extvar_  = distname_(distargs__)) ||
                    @capture(elt, extvar_ += distname_(distargs__)) ||
                    @capture(elt, extvar_ *= distname_(distargs__)))

                # For "anonymous" RVs, e.g., ext_var2[2010:2100, :] *= Uniform(0.8, 1.2), we
                # gensym a name based on the external var name and process it as a named RV.
                if rvname === nothing
                    param_name = @capture(extvar, name_[args__]) ? name : extvar
                    rvname = _make_rvname(param_name)
                    saverv(rvname, distname, distargs)
                end

                op = elt.head
                if @capture(extvar, name_[args__])
                    # println("Ref:  $name, $args")
                    # Meta.show_sexpr(extvar)
                    # println("")

                    # if extvar.head == :ref, extvar.args must be one of:
                    # - a scalar value, e.g., name[2050] => (:ref, :name, 2050)
                    #   convert to tuple of dimension specifiers (:name, 2050)
                    # - a slice expression, e.g., name[2010:2050] => (:ref, :name, (:(:), 2010, 2050))
                    #   convert to (:name, 2010:2050) [convert it to actual UnitRange instance]
                    # - a tuple of symbols, e.g., name[(US, CHI)] => (:ref, :name, (:tuple, :US, :CHI))
                    #   convert to (:name, (:US, :CHI))
                    # - combinations of these, e.g., name[2010:2050, (US, CHI)] => (:ref, :name, (:(:), 2010, 2050), (:tuple, :US, :CHI))
                    #   convert to (:name, 2010:2050, (:US, :CHI))

                    dims = _make_dims(args)
                    expr = :(TransformSpec($(QuoteNode(name)), $(QuoteNode(op)), $(QuoteNode(rvname)), [$(dims...)]))
                else
                    expr = :(TransformSpec($(QuoteNode(extvar)), $(QuoteNode(op)), $(QuoteNode(rvname))))
                end
                push!(_transforms, esc(expr))
            else
                error("Unrecognized expression '$elt' in @defmcs")
            end
        end
        return :(MonteCarloSimulation([$(_rvs...)],
                                      [$(_transforms...)],
                                      CorrelationSpec[$(_corrs...)],
                                      Tuple{Symbol, Symbol}[$(_saves...)]))
    end
end


text:
using BinDeps

isfile("deps.jl") && rm("deps.jl")

@BinDeps.setup
    libtriangle = library_dependency("libtriangle", aliases = ["libtriangle.dylib"], runtime = true)

    rootdir = BinDeps.depsdir(libtriangle)
    srcdir = joinpath(rootdir, "src")
    prefix = joinpath(rootdir, "usr")
    libdir = joinpath(prefix, "lib")
    headerdir = joinpath(prefix, "include")

    if Sys.iswindows()
        libfile = joinpath(libdir, "libtriangle.dll")
        arch = "x86"
        if Sys.WORD_SIZE == 64
            arch = "x64"
        end
        @build_steps begin
            FileRule(libfile, @build_steps begin
                 BinDeps.run(@build_steps begin
                    ChangeDirectory(srcdir)
                    `cmd /c compile.bat all $arch`
                    `cmd /c copy libtriangle.dll $libfile`
                    `cmd /c copy triangle.h $headerdir`
                    `cmd /c copy tricall.h $headerdir`
                    `cmd /c copy commondefine.h $headerdir`
                    `cmd /c compile.bat clean $arch`
            end) end) end

        provides(Binaries, URI(libfile), libtriangle)
    else
        libname = "libtriangle.so"
        if Sys.isapple()
            libname = "libtriangle.dylib"
        end
        libfile = joinpath(libdir, libname)
        provides(BinDeps.BuildProcess, (@build_steps begin
                    FileRule(libfile, @build_steps begin
                        BinDeps.ChangeDirectory(srcdir)
                        `make clean`
                        `make`
                        `cp libtriangle.so $libfile`
                        `cp triangle.h $headerdir/`
                        `cp tricall.h $headerdir/`
                        `cp commondefine.h $headerdir/`
                        `make clean`
                    end)
                end), libtriangle)
    end

@BinDeps.install Dict(:libtriangle => :libtriangle)


text:
using MultivariateOrthogonalPolynomials, ArrayLayouts, BandedMatrices, Test
import MultivariateOrthogonalPolynomials: ModalInterlace, ModalInterlaceLayout

@testset "ModalInterlace" begin
    ops = [brand(2,3,1,2), brand(1,2,1,1), brand(1,2,1,2)]
    A = ModalInterlace(ops, (3,5), (2,4))
    @test MemoryLayout(A) isa ModalInterlaceLayout
    @test A[[1,4],[1,4,11]] == ops[1]
    @test A[[2],[2,7]] == ops[2]
    @test A[[5],[5,12]] == A[[6],[6,13]] == ops[3]
end


text:
#!/usr/local/bin/julia
#
const PKGS = ["RDatasets","CSV","JLD","HDF5","Query","Match","PyPlot","GR","StatsBase",
	      "Distributions","GLM","HypothesisTests","KernelDensity","Clustering"]

println("Installing required packages ...")

using Pkg
ipks = Pkg.installed()
for p in PKGS
    println("> $p")
    !haskey(ipks,p) && Pkg.add(p)
end
println("Done!")


text:
using Documenter, DPMMSubClusters

push!(LOAD_PATH,"../src/")
makedocs(
    sitename = "DPMMSubClusters.jl",
    doctest  = false,
    pages    = [
        "index.md",
        "getting_started.md",
        "usage.md",
        "priors.md",
        "data_generation.md",
        "perf.md"
    ]
)

deploydocs(
    repo = "github.com/dinarior/DPMMSubClusters.jl.git",
)


text:
struct SMILESParseException <: Exception
    msg::String
    smiles::String
    charloc::Int
end
Base.showerror(io::IO, e::SMILESParseException) = print(io, "Failed to parse SMILES on character ", e.charloc, "!")

function ParseSMILES(s::AbstractString, calculate_implicit_hydrogens = true )
    MoleculeGraph = SimpleWeightedGraph();
    MolecularData = Element[]

    len = length(s)

    chainstack = Int16[]

    nextedgestart = 0
    weight = 1
    lastcursor = 0
    RingClosures = Dict{Int,Vector{Int}}()

    idx = 1
    while idx <= len
        moiety = nothing
        cursor = s[idx]

        if isdigit(cursor) || cursor == '%' #Handle rings
            if isdigit(cursor)
                # Single-digit ring id
                id = Int16(cursor - '0')
                idx = nextind(s, idx)
            else
                # Two-digit ring id
                id = parse(Int16, s[idx+1:idx+2])
                idx += 3
            end
            push!( MolecularData[end].ringID, id )
            if id in keys( RingClosures )
                push!( RingClosures[ id ], length(MolecularData) )
            else
                RingClosures[ id ] = [ length(MolecularData) ]
            end
        elseif isletter( cursor )
            symbollist = isuppercase( cursor ) ? aliphatics : aromatics
            aromatic = islowercase( cursor )
            symbol, idx = tryparseelement(s, idx, symbollist)
            if isa(symbol, Nothing)
                throw(SMILESParseException("Invalid SMILES. Unrecognized chemical symbol.", s, idx))
            end
            moiety = Element(symbol)
        elseif isspecialoperator(cursor)    #Handle operators
            if cursor == '['
                idxclose = findnext(']', s, idx)
                idxclose === nothing && throw(SMILESParseException("Invalid SMILES. Bracket is either empty, or does not end.", s, idx))
                moiety = ParseBracket(SubString(s, idx+1, idxclose-1))
                idx = nextind(s, idxclose)
            end
            if cursor == '('
                if lastcursor == 0
                    push!( chainstack, length( MolecularData )  )
                else
                    push!( chainstack, lastcursor  )
                end
                idx = nextind(s, idx)
            end
            if cursor == ')'
                isempty(chainstack) && throw(SMILESParseException("Invalid SMILES. Parenthesis/chain does not have a beginning.", s, idx))
                idx = nextind(s, idx)
                lastcursor = nextedgestart = pop!( chainstack )
            end
        elseif isbondoperator(cursor) #Handle bonds
            weight = bonds[ cursor ]
            idx = nextind(s, idx)
        end
        #New atom/moiety was parsed
        if isa( moiety, Element )
            lastcursor = 0
            SimpleWeightedGraphs.add_vertex!( MoleculeGraph )
            push!( MolecularData, moiety )
            l = length( MolecularData )
            if (l > 1)
                if nextedgestart == 0
                    SimpleWeightedGraphs.add_edge!(MoleculeGraph, l - 1, l, weight )
                else
                    SimpleWeightedGraphs.add_edge!(MoleculeGraph, nextedgestart, l, weight )
                    nextedgestart = 0
                end
                if weight > 1
                    weight = 1
                end
            end
        end
    end
    #Close rings!
    for (k, v) in RingClosures
        for id in 1 : ( length( v ) - 1 )
            if ( v[id+1] - v[id] ) != -1
                SimpleWeightedGraphs.add_edge!(MoleculeGraph, v[id], v[id+1])
            end
        end
    end
    # Add implicit hydrogens to graph structure
    if calculate_implicit_hydrogens
        for ( i, atom ) in enumerate( MolecularData )
            #Find all edges with this atom.
            implicitH = 0
            #Lookup valence
            if atom.symbol in keys( valence )
                Valence = valence[ atom.symbol ]
                BondedElectrons = sum( MoleculeGraph.weights[:,i].nzval )
                Aromaticity = (isa(atom.aromatic, Nothing) ? 0 : atom.aromatic)
                implicitH = Valence .- BondedElectrons .- Aromaticity
                #Look I'm not keeping track of non-bonding e- pairs
                #All I'm doing here is seeing if this is likely invalid structure...
                if length(implicitH) > 1
                    if all( implicitH .< 0 )
                        #Wiggle to find proper valence
                        @warn("Illegal number of implicit hydrogens in $(atom.symbol) (Atom #$i). Defaulting to 0 hydrogens.")
                        implicitH = 0
                    else
                        implicitH = implicitH[ implicitH .> 0 ]
                        implicitH = reduce(min, implicitH .> 0)
                    end
                end
            else
                #Valence not known to package :/
                @warn("The valence for $(atom.symbol) isn't known, implicit hydrogens not determined.")
            end
            # implicitH >= atom.explicithydrogens || error("more hydrogens were supplied than supported by the valence in $Sorig (at $S)")
            MolecularData[i].implicithydrogens = max(0, implicitH - atom.explicithydrogens)
        end
    end

    return MoleculeGraph, MolecularData
end


text:
# Automatically generated using Clang.jl wrap_c

# Skipping MacroDefinition: ROSIDL_GET_ACTION_TYPE_SUPPORT ( PkgName , Name ) ROSIDL_TYPESUPPORT_INTERFACE__ACTION_SYMBOL_NAME ( rosidl_typesupport_c , PkgName , action , Name ) ( )

const rosidl_service_typesupport_handle_function = Ptr{Cvoid}

struct rosidl_service_type_support_t
    typesupport_identifier::Cstring
    data::Ptr{Cvoid}
    func::rosidl_service_typesupport_handle_function
end

const rosidl_message_typesupport_handle_function = Ptr{Cvoid}

struct rosidl_message_type_support_t
    typesupport_identifier::Cstring
    data::Ptr{Cvoid}
    func::rosidl_message_typesupport_handle_function
end

struct rosidl_action_type_support_t
    goal_service_type_support::Ptr{rosidl_service_type_support_t}
    result_service_type_support::Ptr{rosidl_service_type_support_t}
    cancel_service_type_support::Ptr{rosidl_service_type_support_t}
    feedback_message_type_support::Ptr{rosidl_message_type_support_t}
    status_message_type_support::Ptr{rosidl_message_type_support_t}
end

@cenum(rosidl_runtime_c_message_initialization,
    ROSIDL_RUNTIME_C_MSG_INIT_ALL = 0,
    ROSIDL_RUNTIME_C_MSG_INIT_SKIP = 1,
    ROSIDL_RUNTIME_C_MSG_INIT_ZERO = 2,
    ROSIDL_RUNTIME_C_MSG_INIT_DEFAULTS_ONLY = 3,
)

# Skipping MacroDefinition: ROSIDL_GET_MSG_TYPE_SUPPORT ( PkgName , MsgSubfolder , MsgName ) ROSIDL_TYPESUPPORT_INTERFACE__MESSAGE_SYMBOL_NAME ( rosidl_typesupport_c , PkgName , MsgSubfolder , MsgName ) ( )
# Skipping MacroDefinition: ROSIDL_GENERATOR_C__PRIMITIVE_SEQUENCE ( STRUCT_NAME , TYPE_NAME ) typedef struct rosidl_generator_c__ ## STRUCT_NAME ## __Sequence { TYPE_NAME * data ; size_t size ; /*!< The number of valid items in data */ size_t capacity ; /*!< The number of allocated items in data */ } rosidl_generator_c__ ## STRUCT_NAME ## __Sequence ;

struct rosidl_generator_c__bool__Sequence
    data::Ptr{Bool}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__byte__Sequence
    data::Ptr{UInt8}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__char__Sequence
    data::Ptr{Cchar}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__float32__Sequence
    data::Ptr{Cfloat}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__float64__Sequence
    data::Ptr{Cdouble}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__int8__Sequence
    data::Ptr{Int8}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__uint8__Sequence
    data::Ptr{UInt8}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__int16__Sequence
    data::Ptr{Int16}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__uint16__Sequence
    data::Ptr{UInt16}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__int32__Sequence
    data::Ptr{Int32}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__uint32__Sequence
    data::Ptr{UInt32}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__int64__Sequence
    data::Ptr{Int64}
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__uint64__Sequence
    data::Ptr{UInt64}
    size::Csize_t
    capacity::Csize_t
end

# Skipping MacroDefinition: ROSIDL_GENERATOR_C__DECLARE_PRIMITIVE_SEQUENCE_FUNCTIONS ( STRUCT_NAME , TYPE_NAME ) ROSIDL_GENERATOR_C_PUBLIC bool rosidl_generator_c__ ## STRUCT_NAME ## __Sequence__init ( rosidl_generator_c__ ## STRUCT_NAME ## __Sequence * sequence , size_t size ) ; ROSIDL_GENERATOR_C_PUBLIC void rosidl_generator_c__ ## STRUCT_NAME ## __Sequence__fini ( rosidl_generator_c__ ## STRUCT_NAME ## __Sequence * sequence ) ;
# Skipping MacroDefinition: ROSIDL_GET_SRV_TYPE_SUPPORT ( PkgName , SrvSubfolder , SrvName ) ROSIDL_TYPESUPPORT_INTERFACE__SERVICE_SYMBOL_NAME ( rosidl_typesupport_c , PkgName , SrvSubfolder , SrvName ) ( )

struct rosidl_generator_c__String
    data::Cstring
    size::Csize_t
    capacity::Csize_t
end

struct rosidl_generator_c__String__Sequence
    data::Ptr{rosidl_generator_c__String}
    size::Csize_t
    capacity::Csize_t
end

# Skipping MacroDefinition: ROSIDL_GENERATOR_C_EXPORT __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: ROSIDL_GENERATOR_C_PUBLIC __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: ROSIDL_GENERATOR_C_LOCAL __attribute__ ( ( visibility ( "hidden" ) ) )


text:
module RocketToArrayOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: to_array()" begin

    println("Testing: operator to_array()")

    run_proxyshowcheck("ToArray", to_array())

    run_testset([
        (
            source = from(1:5) |> to_array(),
            values = @ts([ [ 1, 2, 3, 4, 5 ], c ]),
            source_type = Vector{Int}
        ),
        (
            source = from("Hello, world") |> to_array(),
            values = @ts([ ['H', 'e', 'l', 'l', 'o', ',', ' ', 'w', 'o', 'r', 'l', 'd'], c ]),
            source_type = Vector{Char},
        ),
        (
            source = of('a') |> to_array(),
            values = @ts([ [ 'a' ], c ]),
            source_type = Vector{Char}
        ),
        (
            source = completed() |> to_array(),
            values = @ts([ [], c ])
        ),
        (
            source = faulted(1) |> to_array(),
            values = @ts(e(1))
        ),
        (
            source = never() |> to_array(),
            values = @ts()
        )
    ])

end

end


text:
using DocumentDistances: findclosetoken, filtertokens
using DocumentDistances: FilteredNGramDocument, ngrams
using TextAnalysis

@testset "TextAnalysis interface" begin

@testset "findclosetoken" begin
    @test findclosetoken("president") == "president"
    @test findclosetoken("president.") == "president"
    @test findclosetoken("president-") == "president"
    @test findclosetoken("president,") == "president"

    @test findclosetoken("President") == "president"
    @test findclosetoken("President.") == "president"
    @test findclosetoken("President-") == "president"
    @test findclosetoken("President,") == "president"

    @test findclosetoken("ewirjp2or,") === nothing
end # @testset "findclosetoken"

@testset "filtertokens" begin
    res = filtertokens(["a", "President.", "Robert,"])
    @test sort(res) == sort(["a", "president", "robert"])
end # @testset "filtertokens"

@testset "FilteredNGramDocument from string and other docs" begin
    d = FilteredNGramDocument("a President, Robert, a Swedish programmer.")
    ng = ngrams(d)
    @test sort(collect(keys(ng))) == sort(AbstractString["a", ",", ".",
        "president", "robert", "swedish", "programmer"])
    @test ng["a"] == 2
    @test ng["president"] == 1
    @test ng["robert"] == 1
    @test ng["swedish"] == 1
    @test ng["programmer"] == 1
    @test ng["."] == 1
    @test ng[","] == 2

    d2 = FilteredNGramDocument(TokenDocument(["a", "President", "Robert", "a", "Swedish", "programmer"]))
    ng2 = ngrams(d2)
    @test sort(collect(keys(ng2))) == sort(AbstractString["a", "president", "robert", "swedish", "programmer"])
    @test ng2["a"] == 2

    ExampleFile = joinpath(dirname(@__FILE__()), "data", "example_with_strange_tokens.txt")
    d3 = FilteredNGramDocument(FileDocument(ExampleFile))
    @test sort(collect(keys(ngrams(d3)))) == sort(AbstractString["1", "2", "parl", "was", "modern", "1,2", "arne", ".", "best", "an", "success", "reported", "to", "sloan", "interpretative", "is", "mo", "and", "software", "sandor", "large-scale", "study", "methods", ":", "june", "development", "lea", "on", "abstract", "involvement", "adaptations", "coordination", "printed", "companies", "not", "guaranteed", "but", "for", "very", ",", "tory", "paper", "customer", "case", "co-located", "findings", "suit", "believed", "areas", "are", "done", "small", "at", "one", "imploring", "&", "it", "ebbe", "key", "nor", "the", "archi", "1977", "ro", "good", "3"])
end # @testset "FilteredNGramDocument from string"

@testset "FilteredNGramDocument document distance" begin
    wmd = WordMoversDistance()
    ws1 = ["obama", "speaks", "media", "illinois"]
    d1 = FilteredNGramDocument(join(ws1, " "), wmd)
    ws2 = ["president", "greets", "press", "chicago"]
    d2 = FilteredNGramDocument(join(ws2, " "), wmd)
    @test evaluate(wmd, d1, d2) == evaluate(wmd, ws1, ws2)
end # @testset "FilteredNGramDocument doc distance"

end


text:
using ReinforcementLearningEnvironmentGym, Flux, ReinforcementLearning
# List all envs

listallenvs()

# CartPole example

env = GymEnv("CartPole-v0")
learner = DQN(Chain(Dense(4, 48, relu), Dense(48, 24, relu), Dense(24, 2)),
                  updateevery = 1, updatetargetevery = 100,
                  startlearningat = 50, minibatchsize = 32,
                  doubledqn = false, replaysize = 10^3,
                  opttype = x -> ADAM(x, .0005))
x = RLSetup(learner, env, ConstantNumberEpisodes(10),
            callbacks = [Progress(), EvaluationPerEpisode(TimeSteps()),
                         Visualize(wait = 0)])
@info("Before learning.")
run!(x)
pop!(x.callbacks)
x.stoppingcriterion = ConstantNumberEpisodes(400)
@time learn!(x)
x.stoppingcriterion = ConstantNumberEpisodes(10)
push!(x.callbacks, Visualize(wait = 0))
@info("After learning.")
run!(x)