text:
using DataFrames, DataFramesMeta
import Base.Threads.@threads
"""
    Nucleotide_Essentials.FilterQuality_se

Filters an input .fastq file based upon the encoded Phred+33 or Phred+64 quality scores. The encoding of the reads is automatically deteremined by looking for unique encoding in Phred+33 and Phred+64. Phred+64 encoding is identified by searching for `^`, `a`, `]`, and `f`.

Reads are filtered based upon the number of [expected errors](https://academic.oup.com/bioinformatics/article/31/21/3476/194979?login=true) (``\\mathrm{E}``) based on the error rate based on quality score and the sum of error probabilities, following the equation:

``\\mathrm{E} = \\sum{_ip_i} = \\sum{_i}10^{\\frac{-Q_i}{10}}``

Stringent filtering (```maxEE = 1```) is used by default but can be adjusted by the user.

Reads that pass the filtering parameters are output to a file ending in `_FilteredReads.fastq` in the user-determined directory, as indicated by ```out```.

Supported keyword arguments include:

* ```read1::String```: Path to the reads to undergo quality filtering
* ```out::String```: Path to the directory where reads that pass the quality filtering should be written
* ```maxEE::Int64``` (optional): The max number of expected errors a read can include as the filtering parameter (default: ```maxEE = 1```)
* ```verbose::Bool``` (optional): Whether or not to show some intermediary feedback on the progress of the function (default = false)

# Example:
```julia
FilterQuality_se("forward_R1.fasta", "/outdirectory")
```
"""
function FilterQuality_se(read1::String, out::String, maxEE::Int64 = 1, verbose::Bool = false)

    if endswith(read1, ".gz") === true
        if verbose == true
            println(string("Unzipping ", read1))
        end
        run(`gunzip $read1`); #this works
        unzipped = true
    else
        unzipped = false
    end

    R1 = readFastq(read1)
    quals = R1.quality
    seqs = R1.sequence
    IDs = R1.ID
    Passing = 0
    Failing = 0

    if occursin("^", join(quals)) == true && occursin("a", join(quals)) == true && occursin("]", join(quals)) == true && occursin("f", join(quals)) == true
        encoding = 64;
        println("Detected Phred+64 quality encoding")
    end
    if occursin("^", join(quals)) == false && occursin("a", join(quals)) == false && occursin("]", join(quals)) == false && occursin("f", join(quals)) == false
        println("Detected Phred+33 quality encoding")
        encoding = 33;
    end

    if encoding == 33
        scores_33 = Dict(
            ('!') => "1.0000",
            ('\"') => "0.79433",
            ('#') => "0.63096",
            ('\$') => "0.50119",
            ('%') => "0.39811",
            ('&') => "0.31623",
            ('\'') => "0.25119",
            ('(') => "0.19953",
            (')') => "0.15849",
            ('*') => "0.12589",
            ('+') => "0.10000",
            (',') => "0.07943",
            ('-') => "0.06310",
            ('.') => "0.05012",
            ('/') => "0.03981",
            ('0') => "0.03162",
            ('1') => "0.02512",
            ('2') => "0.01995",
            ('3') => "0.01585",
            ('4') => "0.01259",
            ('5') => "0.01000",
            ('6') => "0.00794",
            ('7') => "0.00631",
            ('8') => "0.00501",
            ('9') => "0.00398",
            (':') => "0.00316",
            (';') => "0.00251",
            ('<') => "0.00200",
            ('=') => "0.00158",
            ('>') => "0.00126",
            ('?') => "0.00100",
            ('@') => "0.00079",
            ('A') => "0.00063",
            ('B') => "0.00050",
            ('C') => "0.00040",
            ('D') => "0.00032",
            ('E') => "0.00025",
            ('F') => "0.00020",
            ('G') => "0.00016",
            ('H') => "0.00013",
            ('I') => "0.00010",
            ('J') => "0.00008",
            ('K') => "0.00006"
        )
    end
    if encoding == 64
        scores_64 = Dict(
            ('@') => "1.0000",
            ('A') => "0.79433",
            ('B') => "0.63096",
            ('C') => "0.50119",
            ('D') => "0.39811",
            ('E') => "0.31623",
            ('F') => "0.25119",
            ('G') => "0.19953",
            ('H') => "0.15849",
            ('I') => "0.12589",
            ('J') => "0.10000",
            ('K') => "0.07943",
            ('L') => "0.06310",
            ('M') => "0.05012",
            ('N') => "0.03981",
            ('O') => "0.03162",
            ('P') => "0.02512",
            ('Q') => "0.01995",
            ('R') => "0.01585",
            ('S') => "0.01259",
            ('T') => "0.01000",
            ('U') => "0.00794",
            ('V') => "0.00631",
            ('W') => "0.00501",
            ('X') => "0.00398",
            ('Y') => "0.00316",
            ('Z') => "0.00251",
            ('[') => "0.00200",
            ('\\') => "0.00158",
            (']') => "0.00126",
            ('^') => "0.00100",
            ('_') => "0.00079",
            ('`') => "0.00063",
            ('a') => "0.00050",
            ('b') => "0.00040",
            ('c') => "0.00032",
            ('d') => "0.00025",
            ('e') => "0.00020",
            ('f') => "0.00016",
            ('g') => "0.00013",
            ('h') => "0.00010",
            ('i') => "0.00008",
            ('j') => "0.00006"
        )
    end

    alreadyprinted = false
    alreadyprinted1 = false
    alreadyprinted2 = false
    alreadyprinted3 = false
    alreadyprinted4 = false
    alreadyprinted5 = false
    alreadyprinted6 = false
    alreadyprinted7 = false
    alreadyprinted8 = false
    alreadyprinted9 = false
    alreadyprinted10 = false
    alreadyprinted11= false
    alreadyprinted12 = false
    alreadyprinted13= false
    alreadyprinted14 = false
    alreadyprinted15 = false
    alreadyprinted16 = false
    alreadyprinted17 = false
    alreadyprinted18 = false

    if verbose == true
        println("Prepping the reads for quality filtering")
    end

    function convert_phred(c)
        if encoding == 33
            scores_33[c]
        else
            scores_64[c]
        end
    end

    for i in 1:length(quals)
        ee1 = sum(parse.(Float64, (map(convert_phred, collect(quals[i])))))
        if ee1 >= maxEE
            Failing = Failing + 1
        else
            Passing = Passing + 1
            if isfile(string(out, "/",split(R1.filename, ".")[1], "_filtered.fastq")) == false
                out_string = string("@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                write(string(out, "/",split(R1.filename, ".")[1], "_filtered.fastq"), out_string)
            else
                out_string2 = string("\n@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                iostream = open(string(out, "/",split(R1.filename, ".")[1], "_filtered.fastq"), "a")
                write(iostream, out_string2)
                close(iostream)
            end
        end
        if verbose == true
            if (round((i/length(seqs))*100)) > 5 && (round((i/length(seqs))*100)) < 10 && alreadyprinted == false
                alreadyprinted = true
                println("Filtering 5% Complete")
            end
            if (round((i/length(seqs))*100)) > 10 && (round((i/length(seqs))*100)) < 15 && alreadyprinted1 == false
                alreadyprinted1 = true
                println("Filtering 10% Complete")
            end
            if (round((i/length(seqs))*100)) > 15 && (round((i/length(seqs))*100)) < 20 && alreadyprinted2 == false
                alreadyprinted2 = true
                println("Filtering 15% Complete")
            end
            if (round((i/length(seqs))*100)) > 20 && (round((i/length(seqs))*100)) < 25 && alreadyprinted3 == false
                alreadyprinted3 = true
                println("Filtering 20% Complete")
            end
            if (round((i/length(seqs))*100)) > 25 && (round((i/length(seqs))*100)) < 30 && alreadyprinted4 == false
                alreadyprinted4 = true
                println("Filtering 25% Complete")
            end
            if (round((i/length(seqs))*100)) > 30 && (round((i/length(seqs))*100)) < 35 && alreadyprinted5 == false
                alreadyprinted5 = true
                println("Filtering 30% Complete")
            end
            if (round((i/length(seqs))*100)) > 35 && (round((i/length(seqs))*100)) < 40 && alreadyprinted6 == false
                alreadyprinted6 = true
                println("Filtering 35% Complete")
            end
            if (round((i/length(seqs))*100)) > 40 && (round((i/length(seqs))*100)) < 45 && alreadyprinted7 == false
                alreadyprinted7 = true
                println("Filtering 40% Complete")
            end
            if (round((i/length(seqs))*100)) > 45 && (round((i/length(seqs))*100)) < 50 && alreadyprinted8 == false
                alreadyprinted8 = true
                println("Filtering 45% Complete")
            end
            if (round((i/length(seqs))*100)) > 50 && (round((i/length(seqs))*100)) < 55 && alreadyprinted9 == false
                alreadyprinted9 = true
                println("Filtering 50% Complete")
            end
            if (round((i/length(seqs))*100)) > 55 && (round((i/length(seqs))*100)) < 60 && alreadyprinted10 == false
                alreadyprinted10 = true
                println("Filtering 55% Complete")
            end
            if (round((i/length(seqs))*100)) > 60 && (round((i/length(seqs))*100)) < 65 && alreadyprinted11 == false
                alreadyprinted11 = true
                println("Filtering 60% Complete")
            end
            if (round((i/length(seqs))*100)) > 65 && (round((i/length(seqs))*100)) < 70 && alreadyprinted12 == false
                alreadyprinted12 = true
                println("Filtering 65% Complete")
            end
            if (round((i/length(seqs))*100)) > 70 && (round((i/length(seqs))*100)) < 75 && alreadyprinted13 == false
                alreadyprinted13 = true
                println("Filtering 70% Complete")
            end
            if (round((i/length(seqs))*100)) > 75 && (round((i/length(seqs))*100)) < 80 && alreadyprinted14 == false
                alreadyprinted14 = true
                println("Filtering 75% Complete")
            end
            if (round((i/length(seqs))*100)) > 80 && (round((i/length(seqs))*100)) < 85 && alreadyprinted15 == false
                alreadyprinted15 = true
                println("Filtering 80% Complete")
            end
            if (round((i/length(seqs))*100)) > 85 && (round((i/length(seqs))*100)) < 90 && alreadyprinted16 == false
                alreadyprinted16 = true
                println("Filtering 85% Complete")
            end
            if (round((i/length(seqs))*100)) > 90 && (round((i/length(seqs))*100)) < 95 && alreadyprinted17 == false
                alreadyprinted17 = true
                println("Filtering 90% Complete")
            end
            if (round((i/length(seqs))*100)) > 95 && (round((i/length(seqs))*100)) < 100 && alreadyprinted18 == false
                alreadyprinted18 = true
                println("Filtering 95% Complete")
            end
        end
    end

    if verbose == true
        println("Cleaning things up ")
    end

    # Zip the outputted files
    outfile1 = string(out, "/",split(R1.filename, ".")[1], "_filtered.fastq")
    run(`gzip $outfile1`); #this works

    if unzipped == true
        if verbose == true
            println(string("Zipping ", read1))
        end
        run(`gzip $read1`); #this works
    end

    println("")
    println("")
    printstyled(string(Passing, " (", (round(Passing/length(seqs)*100, digits = 2)), "%) reads passed quality filtering\n"), color = :green)
    println("")
    println("")
    printstyled(string(Failing, " (", (round(Failing/length(seqs)*100, digits = 2)), "%) reads failed quality filtering\n"), color = :red)
end

"""
    Nucleotide_Essentials.FilterQuality_pe

Filters an input .fastq file based upon the encoded Phred+33 or Phred+64 quality scores. The encoding of the reads is automatically deteremined by looking for unique encoding in Phred+33 and Phred+64. Phred+64 encoding is identified by searching for `^`, `a`, `]`, and `f`.

Reads are filtered based upon the number of [expected errors](https://academic.oup.com/bioinformatics/article/31/21/3476/194979?login=true) (``\\mathrm{E}``) based on the error rate based on quality score and the sum of error probabilities, following the equation:

``\\mathrm{E} = \\sum{_ip_i} = \\sum{_i}10^{\\frac{-Q_i}{10}}``

Stringent filtering (```maxEE = 1```) is used by default but can be adjusted by the user.

Output Files:

* If both the forward and reverse reads pass the filtering parameters:
    * Forward reads are output to a file ending in `R1_Paired_filtered.fastq` in the user-determined directory, as indicated by ```out```
    * Reverse reads are output to a file ending in `R2_Paired_filtered.fastq` in the user-determined directory, as indicated by ```out```
* If only the forward read passes the filtering parameters:
    * Forward reads are output to a file ending in `R1_Unpaired_filtered.fastq` in the user-determined directory, as indicated by ```out```
    * Reverse reads are not written to a file
* If only the reverse read passes the filtering parameters:
    * Reverse reads are output to a file ending in `R2_Unpaired_filtered.fastq` in the user-determined directory, as indicated by ```out```
    * Forward reads are not written to a file

Supported keyword arguments include:

* ```read1::String```: Path to the forward reads to undergo quality filtering
* ```read2::String```: Path to the reverse reads to undergo quality filtering
* ```out::String```: Path to the directory where reads that pass the quality filtering should be written
* ```maxEE::Int64``` (optional): The max number of expected errors a read can include as the filtering parameter (default: ```maxEE = 1```)
* ```verbose::Bool``` (optional): Whether or not to show some intermediary feedback on the progress of the function (default = false)

# Example:
```julia
FilterQuality_pe("forward_R1.fasta", "reverse_R2.fasta", "/outdirectory")

# changing the filtering parameters
FilterQuality_pe("forward_R1.fasta", "reverse_R2.fasta", "/outdirectory", 2, true)
```
"""
function FilterQuality_pe(read1::String, read2::String, out::String, maxEE::Int64 = 1, verbose::Bool = false)

    if endswith(read1, ".gz") === true
        if verbose == true
            println(string("Unzipping ", read1))
        end
        run(`gunzip $read1`); #this works
        unzipped1 = true
    else
        unzipped1 = false
    end
    if endswith(read2, ".gz") === true
        if verbose == true
            println(string("Unzipping ", read2))
        end
        run(`gunzip $read2`); #this works
        unzipped2 = true
    else
        unzipped2 = false
    end

    # Forward reads
    if endswith(read1, ".gz") === true
        R1_fixed = replace(read1, ".fastq.gz"=>".fastq")
        R1 = readFastq(R1_fixed)
    else
        R1 = readFastq(read1)
    end
    quals = R1.quality
    seqs = R1.sequence
    IDs = R1.ID

    # Reverse reads
    if endswith(read2, ".gz") === true
        R2_fixed = replace(read2, ".fastq.gz"=>".fastq")
        R2 = readFastq(R2_fixed)
    else
        R2 = readFastq(read2)
    end
    quals2 = R2.quality
    seqs2 = R2.sequence
    IDs2 = R2.ID

    # Empty objects for tracking passing/failing reads
    Passing = 0
    UnpairedR1Succes = 0
    UnpairedR2Succes = 0
    bothfailed = 0

    alreadyprinted = false
    alreadyprinted1 = false
    alreadyprinted2 = false
    alreadyprinted3 = false
    alreadyprinted4 = false
    alreadyprinted5 = false
    alreadyprinted6 = false
    alreadyprinted7 = false
    alreadyprinted8 = false
    alreadyprinted9 = false
    alreadyprinted10 = false
    alreadyprinted11= false
    alreadyprinted12 = false
    alreadyprinted13= false
    alreadyprinted14 = false
    alreadyprinted15 = false
    alreadyprinted16 = false
    alreadyprinted17 = false
    alreadyprinted18 = false

    if occursin("^", join(quals)) == true && occursin("a", join(quals)) == true && occursin("]", join(quals)) == true && occursin("f", join(quals)) == true
    encoding = 64;
    println("Detected Phred+64 quality encoding")
    end
    if occursin("^", join(quals)) == false && occursin("a", join(quals)) == false && occursin("]", join(quals)) == false && occursin("f", join(quals)) == false
        println("Detected Phred+33 quality encoding")
        encoding = 33;
    end

    if encoding == 33
        scores_33 = Dict(
            ('!') => "1.0000",
            ('\"') => "0.79433",
            ('#') => "0.63096",
            ('\$') => "0.50119",
            ('%') => "0.39811",
            ('&') => "0.31623",
            ('\'') => "0.25119",
            ('(') => "0.19953",
            (')') => "0.15849",
            ('*') => "0.12589",
            ('+') => "0.10000",
            (',') => "0.07943",
            ('-') => "0.06310",
            ('.') => "0.05012",
            ('/') => "0.03981",
            ('0') => "0.03162",
            ('1') => "0.02512",
            ('2') => "0.01995",
            ('3') => "0.01585",
            ('4') => "0.01259",
            ('5') => "0.01000",
            ('6') => "0.00794",
            ('7') => "0.00631",
            ('8') => "0.00501",
            ('9') => "0.00398",
            (':') => "0.00316",
            (';') => "0.00251",
            ('<') => "0.00200",
            ('=') => "0.00158",
            ('>') => "0.00126",
            ('?') => "0.00100",
            ('@') => "0.00079",
            ('A') => "0.00063",
            ('B') => "0.00050",
            ('C') => "0.00040",
            ('D') => "0.00032",
            ('E') => "0.00025",
            ('F') => "0.00020",
            ('G') => "0.00016",
            ('H') => "0.00013",
            ('I') => "0.00010",
            ('J') => "0.00008",
            ('K') => "0.00006"
        )
    end
    if encoding == 64
        scores_64 = Dict(
            ('@') => "1.0000",
            ('A') => "0.79433",
            ('B') => "0.63096",
            ('C') => "0.50119",
            ('D') => "0.39811",
            ('E') => "0.31623",
            ('F') => "0.25119",
            ('G') => "0.19953",
            ('H') => "0.15849",
            ('I') => "0.12589",
            ('J') => "0.10000",
            ('K') => "0.07943",
            ('L') => "0.06310",
            ('M') => "0.05012",
            ('N') => "0.03981",
            ('O') => "0.03162",
            ('P') => "0.02512",
            ('Q') => "0.01995",
            ('R') => "0.01585",
            ('S') => "0.01259",
            ('T') => "0.01000",
            ('U') => "0.00794",
            ('V') => "0.00631",
            ('W') => "0.00501",
            ('X') => "0.00398",
            ('Y') => "0.00316",
            ('Z') => "0.00251",
            ('[') => "0.00200",
            ('\\') => "0.00158",
            (']') => "0.00126",
            ('^') => "0.00100",
            ('_') => "0.00079",
            ('`') => "0.00063",
            ('a') => "0.00050",
            ('b') => "0.00040",
            ('c') => "0.00032",
            ('d') => "0.00025",
            ('e') => "0.00020",
            ('f') => "0.00016",
            ('g') => "0.00013",
            ('h') => "0.00010",
            ('i') => "0.00008",
            ('j') => "0.00006"
        )
    end

    function convert_phred(c)
        if encoding == 33
            scores_33[c]
        else
            scores_64[c]
        end
    end

    if verbose == true
        println("Prepping reads for quality filtering")
    end

    for i in 1:length(quals)
        ee1 = sum(parse.(Float64, (map(convert_phred, collect(quals[i])))))
        ee2 = sum(parse.(Float64, (map(convert_phred, collect(quals2[i])))))
        if ee1 < maxEE && ee2 < maxEE
            Passing = Passing + 1
            if isfile(string(out, "/",split(R1.filename, ".")[1], "_Paired_filtered.fastq")) == false
                out_string1 = string("@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                write(string(out, "/",split(R1.filename, ".")[1], "_Paired_filtered.fastq"), out_string1)
            else
                out_string2 = string("\n@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                iostream = open(string(out, "/",split(R1.filename, ".")[1], "_Paired_filtered.fastq"), "a")
                write(iostream, out_string2)
                close(iostream)
            end
            if isfile(string(out, "/",split(R2.filename, ".")[1], "_Paired_filtered.fastq")) == false
                out_string3 = string("@", IDs2[i], "\n", seqs2[i], "\n+\n", quals2[i])
                write(string(out, "/",split(R2.filename, ".")[1], "_Paired_filtered.fastq"), out_string3)
            else
                out_string4 = string("\n@", IDs2[i], "\n", seqs2[i], "\n+\n", quals2[i])
                iostream = open(string(out, "/",split(R2.filename, ".")[1], "_Paired_filtered.fastq"), "a")
                write(iostream, out_string4)
                close(iostream)
            end
        elseif ee1 < maxEE && ee2 > maxEE
            UnpairedR1Succes = UnpairedR1Succes + 1
            if isfile(string(out, "/",split(R1.filename, ".")[1], "_Unpaired_filtered.fastq")) == false
                out_string5 = string("@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                write(string(out, "/",split(R1.filename, ".")[1], "_Unpaired_filtered.fastq"), out_string5)
            else
                out_string6 = string("\n@", IDs[i], "\n", seqs[i], "\n+\n", quals[i])
                iostream = open(string(out, "/",split(R1.filename, ".")[1], "_Unpaired_filtered.fastq"), "a")
                write(iostream, out_string6)
                close(iostream)
            end
        elseif ee1 >= maxEE && ee2 < maxEE
            UnpairedR2Succes = UnpairedR2Succes + 1
            if isfile(string(out, "/",split(R2.filename, ".")[1], "_Unpaired_filtered.fastq")) == false
                out_string7 = string("@", IDs2[i], "\n", seqs2[i], "\n+\n", quals2[i])
                write(string(out, "/",split(R2.filename, ".")[1], "_Unpaired_filtered.fastq"), out_string7)
            else
                out_string8 = string("\n@", IDs2[i], "\n", seqs2[i], "\n+\n", quals2[i])
                iostream = open(string(out, "/",split(R2.filename, ".")[1], "_Unpaired_filtered.fastq"), "a")
                write(iostream, out_string8)
                close(iostream)
            end
        elseif ee1 >= maxEE && ee2 >= maxEE
            bothfailed = bothfailed + 1
        end
        if verbose == true
            if (round((i/length(seqs))*100)) > 5 && (round((i/length(seqs))*100)) < 10 && alreadyprinted == false
                alreadyprinted = true
                println("Filtering 5% Complete")
            end
            if (round((i/length(seqs))*100)) > 10 && (round((i/length(seqs))*100)) < 15 && alreadyprinted1 == false
                alreadyprinted1 = true
                println("Filtering 10% Complete")
            end
            if (round((i/length(seqs))*100)) > 15 && (round((i/length(seqs))*100)) < 20 && alreadyprinted2 == false
                alreadyprinted2 = true
                println("Filtering 15% Complete")
            end
            if (round((i/length(seqs))*100)) > 20 && (round((i/length(seqs))*100)) < 25 && alreadyprinted3 == false
                alreadyprinted3 = true
                println("Filtering 20% Complete")
            end
            if (round((i/length(seqs))*100)) > 25 && (round((i/length(seqs))*100)) < 30 && alreadyprinted4 == false
                alreadyprinted4 = true
                println("Filtering 25% Complete")
            end
            if (round((i/length(seqs))*100)) > 30 && (round((i/length(seqs))*100)) < 35 && alreadyprinted5 == false
                alreadyprinted5 = true
                println("Filtering 30% Complete")
            end
            if (round((i/length(seqs))*100)) > 35 && (round((i/length(seqs))*100)) < 40 && alreadyprinted6 == false
                alreadyprinted6 = true
                println("Filtering 35% Complete")
            end
            if (round((i/length(seqs))*100)) > 40 && (round((i/length(seqs))*100)) < 45 && alreadyprinted7 == false
                alreadyprinted7 = true
                println("Filtering 40% Complete")
            end
            if (round((i/length(seqs))*100)) > 45 && (round((i/length(seqs))*100)) < 50 && alreadyprinted8 == false
                alreadyprinted8 = true
                println("Filtering 45% Complete")
            end
            if (round((i/length(seqs))*100)) > 50 && (round((i/length(seqs))*100)) < 55 && alreadyprinted9 == false
                alreadyprinted9 = true
                println("Filtering 50% Complete")
            end
            if (round((i/length(seqs))*100)) > 55 && (round((i/length(seqs))*100)) < 60 && alreadyprinted10 == false
                alreadyprinted10 = true
                println("Filtering 55% Complete")
            end
            if (round((i/length(seqs))*100)) > 60 && (round((i/length(seqs))*100)) < 65 && alreadyprinted11 == false
                alreadyprinted11 = true
                println("Filtering 60% Complete")
            end
            if (round((i/length(seqs))*100)) > 65 && (round((i/length(seqs))*100)) < 70 && alreadyprinted12 == false
                alreadyprinted12 = true
                println("Filtering 65% Complete")
            end
            if (round((i/length(seqs))*100)) > 70 && (round((i/length(seqs))*100)) < 75 && alreadyprinted13 == false
                alreadyprinted13 = true
                println("Filtering 70% Complete")
            end
            if (round((i/length(seqs))*100)) > 75 && (round((i/length(seqs))*100)) < 80 && alreadyprinted14 == false
                alreadyprinted14 = true
                println("Filtering 75% Complete")
            end
            if (round((i/length(seqs))*100)) > 80 && (round((i/length(seqs))*100)) < 85 && alreadyprinted15 == false
                alreadyprinted15 = true
                println("Filtering 80% Complete")
            end
            if (round((i/length(seqs))*100)) > 85 && (round((i/length(seqs))*100)) < 90 && alreadyprinted16 == false
                alreadyprinted16 = true
                println("Filtering 85% Complete")
            end
            if (round((i/length(seqs))*100)) > 90 && (round((i/length(seqs))*100)) < 95 && alreadyprinted17 == false
                alreadyprinted17 = true
                println("Filtering 90% Complete")
            end
            if (round((i/length(seqs))*100)) > 95 && (round((i/length(seqs))*100)) < 100 && alreadyprinted18 == false
                alreadyprinted18 = true
                println("Filtering 95% Complete")
            end
        end
    end

    # Zip the outputted files
    if verbose == true
        println("Cleaning things up ")
    end
    outfile1 = string(out, "/",split(R1.filename, ".")[1], "_Paired_filtered.fastq")
    outfile2 = string(out, "/",split(R2.filename, ".")[1], "_Paired_filtered.fastq")
    outfile3 = string(out, "/",split(R1.filename, ".")[1], "_Unpaired_filtered.fastq")
    outfile4 = string(out, "/",split(R2.filename, ".")[1], "_Unpaired_filtered.fastq")
    if isfile(outfile1) == true
        if verbose == true
            println(string("""|
                        |----Zipping """, outfile1))
        end
        run(`gzip $outfile1`);
    end
    if isfile(outfile2) == true
        if verbose == true
            println(string("""|
                        |----Zipping """, outfile2))
        end
        run(`gzip $outfile2`);
    end
    if isfile(outfile3) == true
        if verbose == true
            println(string("""|
                        |----Zipping """, outfile3))
        end
        run(`gzip $outfile3`);
    end
    if isfile(outfile4) == true
        if verbose == true
            println(string("""|
                        |----Zipping """, outfile4))
        end
        run(`gzip $outfile4`);
    end

    if unzipped1 == true
        if verbose == true
            println(string("""|
                        |----Zipping """, R1_fixed))
        end
        run(`gzip $R1_fixed`); #this works
    end
    if unzipped2 == true
        if verbose == true
            println(string("""|
                        |----Zipping """, R2_fixed))
        end
        run(`gzip $R2_fixed`); #this works
    end

    println("")
    printstyled(string(Passing, " (", (round(Passing/length(seqs)*100, digits = 2)), "%) read pairs passed quality filtering\n"), color = :green)
    println("")
    printstyled(string(bothfailed, " (", (round(bothfailed/length(seqs)*100, digits = 2)), "%) read pairs failed quality filtering\n"), color = :red)
    println("")
    printstyled(string(UnpairedR1Succes, " (", (round(UnpairedR1Succes/length(seqs)*100, digits = 2)), "%) forward reads passed without a matching reverse read\n"), color = :yellow)
    println("")
    printstyled(string(UnpairedR2Succes, " (", (round(UnpairedR2Succes/length(seqs)*100, digits = 2)), "%) reverse reads passed without a matching forward read\n"), color = :yellow)
end


text:
# Genetic Algorithms
# ==================
#         objfun: Objective fitness function
#              N: Search space dimensionality
# initPopulation: Search space dimension ranges as a vector, or initial population values as matrix,
#                 or generation function which produce individual population entities.
# populationSize: Size of the population
#  crossoverRate: The fraction of the population at the next generation, not including elite children,
#                 that is created by the crossover function.
#   mutationRate: Probability of chromosome to be mutated
#              ɛ: Positive integer specifies how many individuals in the current generation
#                 are guaranteed to survive to the next generation.
#                 Floating number specifies fraction of population.
#
function ga(objfun::Function, N::Int;
            initPopulation::Individual = ones(N),
            lowerBounds::Union(Nothing, Vector) = nothing,
            upperBounds::Union(Nothing, Vector) = nothing,
            populationSize::Int = 50,
            crossoverRate::Float64 = 0.8,
            mutationRate::Float64 = 0.1,
            ɛ::Real = 0,
            selection::Function = ((x,n)->1:n),
            crossover::Function = ((x,y)->(y,x)),
            mutation::Function = (x->x),
            iterations::Integer = 100*N,
            tol = 0.0,
            tolIter = 10,
            verbose = false,
            debug = false,
            interim = false)

    store = Dict{Symbol,Any}()

    # Setup parameters
    elite = isa(ɛ, Int) ? ɛ : @compat round(Int, ɛ * populationSize)
    fitFunc = inverseFunc(objfun)

    # Initialize population
    individual = getIndividual(initPopulation, N)
    fitness = zeros(populationSize)
    population = Array(typeof(individual), populationSize)
    offspring = similar(population)

    # Generate population
    for i in 1:populationSize
        if isa(initPopulation, Vector)
            population[i] = initPopulation.*rand(N)
        elseif isa(initPopulation, Matrix)
            population[i] = initPopulation[:, i]
        elseif isa(initPopulation, Function)
            population[i] = initPopulation(N) # Creation function
        else
            error("Cannot generate population")
        end
        fitness[i] = fitFunc(population[i])
        debug && println("INIT $(i): $(population[i]) : $(fitness[i])")
    end
    fitidx = sortperm(fitness, rev = true)
    keep(interim, :fitness, copy(fitness), store)

    # Generate and evaluate offspring
    itr = 1
    bestFitness = 0.0
    bestIndividual = 0
    fittol = 0.0
    fittolitr = 1
    while true
        debug && println("BEST: $(fitidx)")

        # Select offspring
        selected = selection(fitness, populationSize)

        # Perform mating
        offidx = randperm(populationSize)
        for i in 1:2:populationSize
            j = (i == populationSize) ? i-1 : i+1
            if rand() < crossoverRate
                debug && println("MATE $(offidx[i])+$(offidx[j])>: $(population[selected[offidx[i]]]) : $(population[selected[offidx[j]]])")
                offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]])
                debug && println("MATE >$(offidx[i])+$(offidx[j]): $(offspring[i]) : $(offspring[j])")
            else
                offspring[i], offspring[j] = population[selected[i]], population[selected[j]]
            end
        end

        # Perform mutation
        for i in 1:populationSize
            if rand() < mutationRate
                debug && println("MUTATED $(i)>: $(offspring[i])")
                mutation(offspring[i])
                debug && println("MUTATED >$(i): $(offspring[i])")
            end
        end

        # Elitism
        if elite > 0
            for i in 1:elite
                subs = rand(1:populationSize)
                debug && println("ELITE $(fitidx[i])=>$(subs): $(population[fitidx[i]]) => $(offspring[subs])")
                offspring[subs] = population[fitidx[i]]
            end
        end

        # New generation
        for i in 1:populationSize
            population[i] = offspring[i]
            fitness[i] = fitFunc(offspring[i])
            debug && println("FIT $(i): $(fitness[i])")
        end
        fitidx = sortperm(fitness, rev = true)
        bestIndividual = fitidx[1]
        curGenFitness = @compat Float64(objfun(population[bestIndividual]))
        fittol = abs(bestFitness - curGenFitness)
        bestFitness = curGenFitness

        keep(interim, :fitness, copy(fitness), store)
        keep(interim, :bestFitness, bestFitness, store)

        # Verbose step
        verbose &&  println("BEST: $(bestFitness): $(population[bestIndividual]), G: $(itr)")

        # Terminate:
        #  if fitness tolerance is met for specified number of steps
        if fittol < tol
            if fittolitr > tolIter
                break
            else
                fittolitr += 1
            end
        else
            fittolitr = 1
        end
        # if number of iterations more then specified
        if itr >= iterations
            break
        end
        itr += 1
    end

    return population[bestIndividual], bestFitness, itr, fittol, store
end

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Rotations from an Earth-Fixed, Earth-Centered (ECEF) reference frame to an
#   Earth-Fixed Inertial (ECI) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_ecef_to_eci

"""
    r_ecef_to_eci([T,] ECEF, ECI, jd_utc::Number [, eop_data])

Compute the rotation from an Earth-Centered, Earth-Fixed (`ECEF`) reference
frame to an Earth-Centered Inertial (`ECI`) reference frame at the Julian Day
[UTC] `jd_utc`. The rotation description that will be used is given by `T`,
which can be `DCM` or `Quaternion`. The ECEF frame is selected by the input
`ECEF` and the `ECI` frame is selected by the input `ECI`. The possible values
are listed below. The model used to compute the rotation is specified by the
selection of the origin and destination frames. Currently, there are two models
supported: IAU-76/FK5 and IAU-2006 with 2010 conventions (CIO and equinox
approaches).

# Rotation description

The rotations that aligns the ECEF with ECI can be described by Direction Cosine
Matrices or Quaternions. This is selected by the parameter `T`. The possible
values are:

- `DCM`: The rotation will be described by a Direction Cosine Matrix.
- `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported.**

# ECEF Frame

The ECEF frame is selected by the parameter `ECEF`. The possible values are:

- `ITRF()`: ECEF will be selected as the International Terrestrial Reference
    Frame (ITRF).
- `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
    frame.
- `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
    System (TIRS).

# ECI Frame

The ECI frame is selected by the parameter `ECI`. The possible values are:

- `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
    reference frame.
- `TOD()`: ECI will be selected as the True of Date (TOD).
- `MOD()`: ECI will be selected as the Mean of Date (MOD).
- `J2000()`: ECI will be selected as the J2000 reference frame.
- `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
    (GCRF).
- `CIRS()`: ECI will be selected as the Celestial Intermediate Reference System
    (CIRS).
- `ERS()`: ECI will be selected as the Earth Reference System (ERS).
- `MOD06()`: ECI will be selected as the Mean of Date (MOD) according to the
    definition in IAU-2006/2010 theory.
- `MJ2000()`: ECI will be selected as the J2000 mean equatorial frame (MJ2000).

!!! note
    The frames `MOD()` and `MOD06()` are virtually the same. However, we
    selected different names to make clear which theory are being used since
    mixing transformation between frames from IAU-76/FK5 and IAU-2006/2010 must
    be performed with caution.

# EOP Data

The conversion between the frames depends on EOP Data (see
[`get_iers_eop`](@ref) and [`read_iers_eop`](@ref)). If IAU-76/FK5 model is
used, then the type of `eop_data` must be [`EOPData_IAU1980`](@ref). Otherwise,
if IAU-2006/2010 model is used, then the type of `eop_data` must be
[`EOPData_IAU2000A`](@ref). The following table shows the requirements for EOP
data given the selected frames.

|   Model                     |  ECEF  |   ECI    |    EOP Data     |
|:----------------------------|:-------|:---------|:----------------|
| IAU-76/FK5                  | `ITRF` | `GCRF`   | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `J2000`  | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `MOD`    | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `TOD`    | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `TEME`   | EOP IAU1980     |
| IAU-76/FK5                  | `PEF`  | `GCRF`   | EOP IAU1980     |
| IAU-76/FK5                  | `PEF`  | `J2000`  | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `MOD`    | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `TOD`    | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `TEME`   | Not required¹   |
| IAU-2006/2010 CIO-based     | `ITRF` | `CIRS`   | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `ITRF` | `GCRF`   | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `TIRS` | `CIRS`   | Not required¹   |
| IAU-2006/2010 CIO-based     | `TIRS` | `GCRF`   | Not required¹ ² |
| IAU-2006/2010 Equinox-based | `ITRF` | `ERS`    | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `ITRF` | `MOD06`  | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `ITRF` | `MJ2000` | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `TIRS` | `ERS`    | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `TIRS` | `MOD06`  | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `TIRS` | `MJ2000` | Not required¹ ³ |

`¹`: In this case, UTC will be assumed equal to UT1 to compute the Greenwich
Mean Sidereal Time. This is an approximation, but should be sufficiently
accurate for some applications. Notice that, if EOP Data is provided, UT1 will
be accurately computed.

`²`: In this case, the terms that account for the free core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

`³`: In this case, the terms that corrects the nutation in obliquity and in
longitude due to the free core nutation will not be available, reducing the
precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the ECEF reference
frame into alignment with the ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(Val(:IAU1980));

julia> r_ecef_to_eci(DCM, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33509e-5
 -0.000797312   0.00106478   0.999999

julia> r_ecef_to_eci(ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33509e-5
 -0.000797312   0.00106478   0.999999

julia> r_ecef_to_eci(PEF(), J2000(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45854e-5
 -0.000796885   0.00106622   0.999999

julia> r_ecef_to_eci(PEF(), J2000(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00133066
 -0.78518      -0.619267     3.45854e-5
 -0.000796879   0.00106623   0.999999

julia> r_ecef_to_eci(Quaternion, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.43631 - 0.000590997⋅i + 0.000305106⋅j + 0.000305106⋅k

julia> eop_IAU2000A = get_iers_eop(Val(:IAU2000A));

julia> r_ecef_to_eci(ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33516e-5
 -0.000797311   0.00106478   0.999999

julia> r_ecef_to_eci(TIRS(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45884e-5
 -0.000796885   0.00106623   0.999999

julia> r_ecef_to_eci(Quaternion, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.43631 - 0.000590997⋅i + 0.000305106⋅j + 0.000305106⋅k
```
"""
@inline function r_ecef_to_eci(
    T_ECEF::T_ECEFs,
    T_ECI::T_ECIs,
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, jd_utc, eop_data)
end

@inline function r_ecef_to_eci(
    T_ECEF::T_ECEFs_IAU_2006,
    T_ECI::T_ECIs_IAU_2006,
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, jd_utc, eop_data)
end

# Specializations for those cases that EOP Data is not needed.
@inline function r_ecef_to_eci(
    T_ECEF::Val{:PEF},
    T_ECI::Union{Val{:J2000}, Val{:MOD}, Val{:TOD}, Val{:TEME}},
    jd_utc::Number
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, jd_utc)
end

@inline function r_ecef_to_eci(
    T_ECEF::Val{:TIRS},
    T_ECI::T_ECIs_IAU_2006,
    jd_utc::Number
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, jd_utc)
end

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                 ITRF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:GCRF},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(jd_utc) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    return r_itrf_to_gcrf_fk5(T, jd_ut1, jd_tt, x_p, y_p, δΔϵ_1980, δΔψ_1980)
end

#                                ITRF => J2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:J2000},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Compute the rotation.
    return r_itrf_to_gcrf_fk5(T, jd_ut1, jd_tt, x_p, y_p, 0, 0)
end

#                                 ITRF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MOD},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(jd_utc) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    r_pef_itrf = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_mod_pef  = r_pef_to_mod_fk5(T, jd_ut1, jd_tt, δΔϵ_1980, δΔψ_1980)

    return compose_rotation(r_pef_itrf, r_mod_pef)
end

#                                 ITRF => TOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:TOD},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    r_pef_itrf = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_tod_pef  = r_pef_to_tod_fk5(T, jd_ut1, jd_tt, δΔψ_1980)

    return compose_rotation(r_pef_itrf, r_tod_pef)
end

#                                 ITRF => TEME
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:TEME},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Compute the rotation.
    r_pef_itrf = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_teme_pef = r_pef_to_teme(T, jd_ut1)

    return compose_rotation(r_pef_itrf, r_teme_pef)
end

#                                 PEF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:GCRF},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(jd_utc) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    r_mod_pef  = r_pef_to_mod_fk5(T, jd_ut1, jd_tt, δΔϵ_1980, δΔψ_1980)
    r_gcrf_mod = r_mod_to_gcrf_fk5(T, jd_tt)

    return compose_rotation(r_mod_pef, r_gcrf_mod)
end

#                                PEF => J2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:J2000},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    r_mod_pef  = r_pef_to_mod_fk5(T, jd_ut1, jd_tt, 0, 0)
    r_gcrf_mod = r_mod_to_gcrf_fk5(T, jd_tt)

    return compose_rotation(r_mod_pef, r_gcrf_mod)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:J2000}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    r_mod_pef  = r_pef_to_mod_fk5(T, jd_ut1, jd_tt, 0, 0)
    r_gcrf_mod = r_mod_to_gcrf_fk5(T, jd_tt)

    return compose_rotation(r_mod_pef, r_gcrf_mod)
end

#                                 PEF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:MOD},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(jd_utc) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    return r_pef_to_mod_fk5(T, jd_ut1, jd_tt, δΔϵ_1980, δΔψ_1980)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:MOD}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    return r_pef_to_mod_fk5(T, jd_ut1, jd_tt, 0, 0)
end

#                                 PEF => TOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:TOD},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(jd_utc) * arcsec_to_rad
    y_p      = eop_data.y(jd_utc) * arcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    return r_pef_to_tod_fk5(T, jd_ut1, jd_tt, δΔψ_1980)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:TOD}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    return r_pef_to_tod_fk5(T, jd_ut1, jd_tt, 0)
end

#                                 PEF => TEME
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:TEME},
    jd_utc::Number,
    eop_data::EOPData_IAU1980
)
    # Get the time in UT1.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)

    # Compute the rotation.
    return r_pef_to_teme(T, jd_ut1)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:TEME}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc

    # Compute the rotation.
    return r_pef_to_teme(T, jd_ut1)
end

################################################################################
#                           IAU-2006/2010 CIO-based
################################################################################

#                                 ITRF => CIRS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:CIRS},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Compute the rotation.
    r_tirs_itrf = r_itrf_to_tirs_iau2006(T, jd_tt, x_p, y_p)
    r_cirs_tirs = r_tirs_to_cirs_iau2006(T, jd_ut1)

    return compose_rotation(r_tirs_itrf, r_cirs_tirs)
end

#                                 ITRF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:GCRF},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad
    dx  = eop_data.dX(jd_utc) * milliarcsec_to_rad
    dy  = eop_data.dY(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    r_tirs_itrf = r_itrf_to_tirs_iau2006(T, jd_tt, x_p, y_p)
    r_cirs_tirs = r_tirs_to_cirs_iau2006(T, jd_ut1)
    r_gcrf_cirs = r_cirs_to_gcrf_iau2006(T, jd_tt, dx, dy)

    return compose_rotation(r_tirs_itrf, r_cirs_tirs, r_gcrf_cirs)
end

#                                 TIRS => CIRS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:CIRS},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)

    # Compute the rotation.
    return r_tirs_to_cirs_iau2006(T, jd_ut1)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:CIRS}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc

    # Compute the rotation.
    return r_tirs_to_cirs_iau2006(T, jd_ut1)
end

#                                 TIRS => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:GCRF},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    dx = eop_data.dX(jd_utc) * milliarcsec_to_rad
    dy = eop_data.dY(jd_utc) * milliarcsec_to_rad

    # Compute the rotation.
    r_cirs_tirs = r_tirs_to_cirs_iau2006(T, jd_ut1)
    r_gcrf_cirs = r_cirs_to_gcrf_iau2006(T, jd_tt, dx, dy)

    return compose_rotation(r_cirs_tirs, r_gcrf_cirs)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:GCRF}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc

    # Get the time in TT.
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    r_cirs_tirs = r_tirs_to_cirs_iau2006(T, jd_ut1)
    r_gcrf_cirs = r_cirs_to_gcrf_iau2006(T, jd_tt)

    return compose_rotation(r_cirs_tirs, r_gcrf_cirs)
end

################################################################################
#                         IAU-2006/2010 equinox-based
################################################################################

#                                 ITRF => ERS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:ERS},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Obtain the correction of the nutation in longitude.
    ~, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_tirs_itrf = r_itrf_to_tirs_iau2006(T, jd_tt, x_p, y_p)
    r_ERS_TIRS = r_tirs_to_ers_iau2006(T, jd_ut1, jd_tt, δΔΨ_2000)

    return compose_rotation(r_tirs_itrf, r_ERS_TIRS)
end

#                                 ITRF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MOD06},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_tirs_itrf = r_itrf_to_tirs_iau2006(T, jd_tt, x_p, y_p)
    r_mod_tirs = r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt, δΔϵ_2000, δΔΨ_2000)

    return compose_rotation(r_tirs_itrf, r_mod_tirs)
end

#                                ITRF => MJ2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MJ2000},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(jd_utc) * arcsec_to_rad
    y_p = eop_data.y(jd_utc) * arcsec_to_rad

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_tirs_itrf  = r_itrf_to_tirs_iau2006(T, jd_tt, x_p, y_p)
    r_mod_tirs   = r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt, δΔϵ_2000, δΔΨ_2000)
    r_mj2000_mod = r_mod_to_mj2000_iau2006(T, jd_tt)

    return compose_rotation(r_tirs_itrf, r_mod_tirs, r_mj2000_mod)
end

# NOTE: We do not implement the conversion ITRF => GCRF using equinox-based
# IAU-2006/2010 theory because the CIO-based approach is faster and more
# precise.

#                                  TIRS => ERS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:ERS},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Obtain the correction of the nutation in longitude.
    ~, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    return r_tirs_to_ers_iau2006(T, jd_ut1, jd_tt, δΔΨ_2000)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:ERS}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    return r_tirs_to_ers_iau2006(T, jd_ut1, jd_tt)
end

#                                 TIRS => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:MOD06},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    return r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt, δΔϵ_2000, δΔΨ_2000)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:MOD06}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    return r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt)
end

#                                TIRS => MJ2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:MJ2000},
    jd_utc::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    jd_ut1 = jd_utc_to_ut1(jd_utc, eop_data)
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, jd_utc)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_mod_tirs   = r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt, δΔϵ_2000, δΔΨ_2000)
    r_mj2000_mod = r_mod_to_mj2000_iau2006(T, jd_tt)

    return compose_rotation(r_mod_tirs, r_mj2000_mod)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:MJ2000}, jd_utc::Number)
    # Since we do not have EOP Data, assume that jd_utc is equal to jd_ut1.
    jd_ut1 = jd_utc
    jd_tt  = jd_utc_to_tt(jd_utc)

    # Compute the rotation.
    r_mod_tirs   = r_tirs_to_mod_iau2006(T, jd_ut1, jd_tt)
    r_mj2000_mod = r_mod_to_mj2000_iau2006(T, jd_tt)

    return compose_rotation(r_mod_tirs, r_mj2000_mod)
end

# NOTE: We do not implement the conversion TIRS => GCRF using equinox-based
# IAU-2006/2010 theory because the CIO-based approach is faster and more
# precise.


text:
#Provide utility functions for interoperability with Compose.jl
#Define conversions from Shapefile.jl shape primitives to their Compose.jl counterparts

#Convert Shapefile rectangle to Compose rectangle
import Compose: rectangle
rectangle{T<:Real}(R::Shapefile.Rect{T}) = rectangle(R.left,R.top,R.right-R.left,R.bottom-R.top)

#Compose polygons cannot be disjoint but Shapefile.Polygons can
#Need to convert Shapefile.Polygon to list of Compose polygons
function convert(::Type{Vector{Compose.Form{Compose.PolygonPrimitive}}},
        shape::Shapefile.Polygon)
    points = {}
    polygons={}
    currentpart=2
    for (i,p) in enumerate(shape.points)
        push!(points, p)
        if i==length(shape.points) || (currentpart≤length(shape.parts) && i==shape.parts[currentpart])
            push!(polygons, polygon([(p.x,p.y) for p in points]))
            currentpart += 1
            points = {}
        end
    end
    polygons
end

polygons(shape::Shapefile.Polygon) = convert(Vector{Compose.Form{Compose.PolygonPrimitive}},
    shape)


text:
module Slicer

using Luxor

include("types.jl")
include("util.jl")
include("stl.jl")
include("intersect.jl")
include("aabb.jl")
include("endmatching.jl")
include("shoelace.jl")

function debug(polygons)
    @svg begin
        for polygon in polygons
            foreach(((a, b),) -> arrow(Point(a[1], a[2]), Point(b[1], b[2]), arrowheadlength=.4, linewidth=.1),
            edges(polygon))
        end
    end
end

function slice(path)
    triangles = load(path)
    polylines = intersect(triangles, Plane((0, 0, 1), (0, 0, 1)))
    polygons = polylines |> match |> simplify |> correct
    debug(polygons)
end

end # module


text:
using Revise
using DifferentialEquations
using JLD2
using Dates
using Infiltrator
using PyPlot
using Bem2d

"""
    derivsconst!()

Time derivatives of stresses, velocities, and state variables for
the constant displacement element case.  Calculation is in place for speed.
"""
function derivsconst!(dudt, u, p, t)
    intidx, nintidx, bemsliptotractotal, els, eta, thetalaw, dc, blockvxglobal, blockvyglobal, dthetadt, dvxdt, dvydt, vx, vy, dtracxglobaldt, dtracyglobaldt, dtracglobaldt = p
    @views Bem2d.multmatvec!(vx, vy, els.rotmatinv[intidx, :, :], u[1:3:end], u[2:3:end])
    @views dtracglobaldt = bemsliptotractotal * [blockvxglobal .- u[1:3:end] blockvyglobal .- u[2:3:end]]'[:]
    @views Bem2d.multmatvec!(dtracxglobaldt, dtracyglobaldt, els.rotmatinv[intidx, :, :], dtracglobaldt[1:2:end], dtracglobaldt[2:2:end])
    for i in 1:nintidx
        @views dthetadt[i] = thetalaw(abs(vx[i]), abs(u[3:3:end][i]), dc)
        @views dvxdt[i] = 1 / (eta / els.normalstress[intidx[i]] + els.a[intidx[i]] / abs(vx[i])) * (dtracxglobaldt[i] / els.normalstress[intidx[i]] - els.b[intidx[i]] * dthetadt[i] / u[3:3:end][i])
        dvydt[i] = 0
    end
    @views Bem2d.multmatvec!(dudt[1:3:end], dudt[2:3:end], els.rotmat[intidx, :, :], dvxdt, dvydt)
    dudt[3:3:end] = dthetadt
    return nothing
end

"""
    derivsquad!()

Time derivatives of stresses, velocities, and state variables for
the quadratic displacement element case.  Calculation is in place for speed.
"""
function derivsquad!(dudt, u, p, t)
    intidx, nintidx, partials, els, eta, thetalaw, dc, blockvxglobal, blockvyglobal, dthetadt, dvxdt, dvydt, vx, vy, dtracxglobaldt, dtracyglobaldt, dtracglobaldt = p
    @views multmatvecquad!(vx, vy, els.rotmatinv[intidx, :, :], u[1:3:end], u[2:3:end])
    @views dtracglobaldt = partials["trac"]["fault"]["fault"] * [blockvxglobal .- u[1:3:end] blockvyglobal .- u[2:3:end]]'[:]
    @views multmatvecquad!(dtracxglobaldt, dtracyglobaldt, els.rotmatinv[intidx, :, :], dtracglobaldt[1:2:end], dtracglobaldt[2:2:end])
    for i in 1:3*nintidx
        elidx = Int64(floor((i - 1) / 3) + 1) # Change w/ every 3rd node
        @views dthetadt[i] = thetalaw(abs(vx[i]), abs(u[3:3:end][i]), dc)
        @views dvxdt[i] = 1 / (eta / els.normalstress[intidx[elidx]] + els.a[intidx[elidx]] / abs(vx[i])) * (dtracxglobaldt[i] / els.normalstress[intidx[elidx]] - els.b[intidx[elidx]] * dthetadt[i] / u[3:3:end][i])
        dvydt[i] = 0
    end
    @views multmatvecquad!(dudt[1:3:end], dudt[2:3:end], els.rotmat[intidx, :, :], dvxdt, dvydt)
    dudt[3:3:end] = dthetadt
    return nothing
end

"""
    qdtopography()

Quasidynamic model of earthquake cycle on a non-planar fault
beneath non-planar topography.
"""
function qdtopography()
    # Constants
    nsteps = 10000
    nfreesurf = 100
    nfault = 200
    printstep = 100
    amplitude = 1.0
    outfilename = string(now()) * ".jld2"
    siay = 365.25 * 24 * 60 * 60
    tspan = (0, siay * 1000)
    abstol = 1e-4
    reltol = 1e-4
    mu = 3e10
    nu = 0.25
    rho = 2700.0
    eta = mu / (2.0 * sqrt(mu / rho))
    dc = 0.05
    blockvelx = 1e-9
    blockvely = 0.0

    # Create fault elements
    els = Bem2d.Elements(Int(1e5))
    faultwidth = 10000

    # Curved fault
    x1, y1, x2, y2 = Bem2d.discretizedline(-7e3, 0e3, 0, 0, nfault)
    y1 = 3e3 * atan.(x1 / 1e3)
    y2 = 3e3 * atan.(x2 / 1e3)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.a[els.endidx + i] = 0.015
        els.b[els.endidx + i] = 0.020
        els.normalstress[els.endidx + i] = 50e7
        els.name[els.endidx + i] = "fault"
    end
    Bem2d.standardize_elements!(els)

    # Topographic free surface
    x1, y1, x2, y2 = Bem2d.discretizedline(-20e3, 0, 20e3, 0, nfreesurf)
    y1 = -1e3 * atan.(x1 / 1e3)
    y2 = -1e3 * atan.(x2 / 1e3)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurftopo"
        els.a[els.endidx + i] = 0.000
        els.b[els.endidx + i] = 0.000
        els.normalstress[els.endidx + i] = 0e6
    end
    Bem2d.standardize_elements!(els)

    # Create convience tools
    idx = Bem2d.getidxdict(els)
    partialsconst = Bem2d.initpartials(els)

    # Calculate slip to traction partials on the fault
    println("Calculating velocity to traction matrix")
    @time _, _, partialsconst["trac"]["fault"]["fault"] = Bem2d.partialsconstdispstress(Bem2d.slip2dispstress, els, idx["fault"], idx["fault"], mu, nu)
    @time _, _, partialsconst["trac"]["fault"]["freesurftopo"] = Bem2d.partialsconstdispstress(Bem2d.slip2dispstress, els, idx["fault"], idx["freesurftopo"], mu, nu)
    @time _, _, partialsconst["trac"]["freesurftopo"]["freesurftopo"] = Bem2d.partialsconstdispstress(Bem2d.slip2dispstress, els, idx["freesurftopo"], idx["freesurftopo"], mu, nu)
    @time _, _, partialsconst["trac"]["freesurftopo"]["fault"] = Bem2d.partialsconstdispstress(Bem2d.slip2dispstress, els, idx["freesurftopo"], idx["fault"], mu, nu)

    #
    # Just trying out some new notation...need to think about it.
    # parC["T"]["fault"]["fault"]
    # parQ["S"]["fault"]["fault"]
    #

    # This matrix allows us to go from fault slip to surface displacements
    # We then need to go from surface displacements to tractions on the fault
    # Do we need tractions on the surface too to go to fault tractions?  No because it's a free surface.
    @time faultsliptosurfacedisp = inv(partialsconst["trac"]["freesurftopo"]["freesurftopo"]) * partialsconst["trac"]["fault"]["freesurftopo"]
    faultsliptosurfacedisptofaulttraction = partialsconst["trac"]["freesurftopo"]["fault"] * faultsliptosurfacedisp
    bemsliptotractotal = partialsconst["trac"]["fault"]["fault"] - faultsliptosurfacedisptofaulttraction

    # CS elements - Euler style stress integration
    nnodes = 1 * nfault
    ics = zeros(3 * nnodes)
    ics[1:3:end] = 1e-3 * blockvelx * ones(nnodes)
    ics[2:3:end] = 0.0 * blockvely * ones(nnodes)
    ics[3:3:end] = 1e8 * ones(nnodes)
    intidx = idx["fault"] # indices of elements to integrate
    nintidx = length(intidx)
    dthetadt = zeros(nintidx)
    dvxdt = zeros(nintidx)
    dvydt = zeros(nintidx)
    vx = zeros(nintidx)
    vy = zeros(nintidx)
    dtracxglobaldt = zeros(nintidx)
    dtracyglobaldt = zeros(nintidx)
    dtracglobaldt = zeros(2 * nintidx)
    p = (intidx, nintidx, bemsliptotractotal, els, eta, Bem2d.thetaaginglaw, dc, blockvelx, blockvely, dthetadt, dvxdt, dvydt, vx, vy, dtracxglobaldt, dtracyglobaldt, dtracglobaldt)
    prob = DifferentialEquations.ODEProblem(derivsconst!, ics, tspan, p)
    integrator = DifferentialEquations.init(prob, Vern7(), abstol = abstol, reltol = reltol)
    @time for i in 1:nsteps
        DifferentialEquations.step!(integrator)
        if mod(i, printstep) == 0
            println("step: " * string(i) * " of " * string(nsteps) * ", time: " * string(integrator.sol.t[end] / siay))
        end
    end
    PyPlot.close("all")
    Bem2d.plotqdtimeseries(integrator.sol, 3, nfault)

    # 3QN elements - Euler style stress integration
    # nnodes = 3 * nfault
    # ics = zeros(3 * nnodes)
    # ics[1:3:end] = 1e-3 * blockvelx * ones(nnodes)
    # ics[2:3:end] = 0.0 * blockvely * ones(nnodes)
    # ics[3:3:end] = 1e8 * ones(nnodes)
    # intidx = collect(1:1:els.endidx) # indices of elements to integrate
    # nintidx = length(intidx)
    # dthetadt = zeros(3 * nintidx)
    # dvxdt = zeros(3 * nintidx)
    # dvydt = zeros(3 * nintidx)
    # vx = zeros(3 * nintidx)
    # vy = zeros(3 * nintidx)
    # dtracxglobaldt = zeros(3 * nintidx)
    # dtracyglobaldt = zeros(3 * nintidx)
    # dtracglobaldt = zeros(2 * 3 * nintidx)
    # p = (intidx, nintidx, partialsquad, els, eta, thetasliplaw, dc, blockvelx, blockvely, dthetadt, dvxdt, dvydt, vx, vy, dtracxglobaldt, dtracyglobaldt, dtracglobaldt)
    # prob = ODEProblem(derivsquad!, ics, tspan, p)
    # integrator = init(prob, Vern7(), abstol = abstol, reltol = reltol)
    # @time for i in 1:nsteps
    #     step!(integrator)
    #     if mod(i, printstep) == 0
    #         println("step: " * string(i) * " of " * string(nsteps) * ", time: " * string(integrator.sol.t[end] / siay))
    #     end
    # end
    # plotqdtimeseriesquad(integrator.sol, 3, nfault)

    # Create new folder and save timestamped .jld2 file to that folder
    basename = string(now())
    outfoldername = "/home/meade/Desktop/data/qdvis/" * basename
    basefilename = outfoldername * "/" * basename
    outfilename = basefilename * ".jld2"
    mkpath(outfoldername)
    @time @save outfilename integrator.sol els mu nu

    # Plot and dave basic time series information
    close("all")
    plotqdtimeseries(integrator.sol, 3, nfault)
    figure(1); savefig(basefilename * "_timeseries.pdf")
    figure(2); savefig(basefilename * "_contours.pdf")

    println("Wrote integration results to:")
    println(outfilename)
end
qdtopography()


text:
# Mutable

"G_M_1"
G_M_1 = 1

"""G_M_2"""
G_M_2 = 1

# Constant

"G_C_1"
const G_C_1 = 1

"""G_C_2"""
const G_C_2 = 1

dummy_1() = ()


text:
# Test residual block
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using Test, LinearAlgebra, InvertibleNetworks

# Input
nx1 = 32
nx2 = 32
nx_in = 8
n_hidden = 16   # same for x and y
batchsize = 2

ny1 = 64
ny2 = 24
ny_in = 1

# Residual blocks
RB = ConditionalResidualBlock(nx1, nx2, nx_in, ny1, ny2, ny_in, n_hidden, batchsize)

function loss(RB, X, Y)
    Zx, Zy = RB.forward(X, Y)
    f = -log_likelihood(Zx) - log_likelihood(Zy)
    ΔZx = -∇log_likelihood(Zx)
    ΔZy = -∇log_likelihood(Zy)
    ΔX, ΔY = RB.backward(ΔZx, ΔZy, X, Y)[1:2]
    return f, ΔX, ΔY, RB.W0.grad, RB.W1.grad, RB.W2.grad, RB.W3.grad,
        RB.b0.grad, RB.b1.grad, RB.b2.grad
end

###################################################################################################
# Gradient tests

# Input image and data
X0 = glorot_uniform(nx1, nx2, nx_in, batchsize)
X = glorot_uniform(nx1, nx2, nx_in, batchsize)
Y0 = glorot_uniform(ny1, ny2, ny_in, batchsize)
Y = glorot_uniform(ny1, ny2, ny_in, batchsize)
dX = X - X0
dY = Y - Y0

# Gradient test w.r.t. input
f0, ΔX, ΔY = loss(RB, X0, Y0)[1:3]
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = loss(RB, X0 + h*dX, Y0 + h*dY)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX) - h*dot(dY, ΔY))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
RB0 = ConditionalResidualBlock(nx1, nx2, nx_in, ny1, ny2, ny_in, n_hidden, batchsize)
RBini = deepcopy(RB0)
dW0 = RB.W0.data - RB0.W0.data
dW1 = RB.W1.data - RB0.W1.data
dW2 = RB.W2.data - RB0.W2.data
dW3 = RB.W3.data - RB0.W3.data
f0, ΔX, ΔY, ΔW0, ΔW1, ΔW2, ΔW3 = loss(RB0, X, Y)[1:7]
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    RB0.W0.data = RBini.W0.data + h*dW0
    RB0.W1.data = RBini.W1.data + h*dW1
    RB0.W2.data = RBini.W2.data + h*dW2
    RB0.W3.data = RBini.W3.data + h*dW3
    f = loss(RB0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW0, ΔW0) - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2) - h*dot(dW3, ΔW3))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
using Arpack, SparseArrays, LinearAlgebra

# define functions used here
include("ades.jl")
include("acre.jl")
include("pascaltriangle.jl")
include("in2b.jl")
include("b2in.jl")
include("epsilon.jl")
include("Hintfunc.jl")
include("cutMsizeEne.jl")

function Hsocfunc(Msize0::Int64, Np::Int64, ksoc::Float64, Omega::Float64)

    Msize = Msize0*2

    matp = zeros(Int,Msize+1,Np+1);
    pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    maxmatp = matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_ns

    # defines vectors and matrices
    vecmbnn = spzeros(Int64,Msize+1);
    vecmbnnj = spzeros(Int64,Msize+1);
    vecmbnnij = spzeros(Int64,Msize+1);
    Hsoc = spzeros(ComplexF64,maxmatp,maxmatp);

    # define a matrix for the Hamiltonian
    for nn = 1:maxmatp

        vecmbnn = in2b(nn,Msize,Np)
        energyij = 0.

        # free terms
        for jj = 1:Msize

            vecmbnnj = ades(jj,vecmbnn)
            if vecmbnnj[Msize+1] == 0
               continue
            end

            for ii = 1:Msize

                vecmbnnij = acre(ii,vecmbnnj)

                energyij = epsilon(ii,jj,ksoc,Omega)
                if isapprox(energyij,0) #energyij == 0 # energy is Int and zero if ii==jj
                   continue #vecmbnnij[1+Msize] = 0
                end

                for mm = 1:nn

                    vecmbmm = in2b(mm,Msize,Np)
                    if vecmbnnij[1:Msize] == vecmbmm[1:Msize]
                       Hsoc[mm,nn] = Hsoc[mm,nn] + energyij*sqrt(vecmbnnij[Msize+1])
                    end

                end

            end

        end

        # # Interactions
        # for ll = 1:Msize
        #
        # end

    end

    # since the Hamiltonian is helmitian
    Hsoc = Hsoc + Hsoc' - spdiagm(diag(Hsoc))

    return Hsoc

end

function diagonaliseHsoc(Hsoc::SparseMatrixCSC{ComplexF64})

    lambda, phi = eigs(Hsoc,nev=6,which=:SR)
    # lambda = real(lambda)

    return lambda, phi

end

function main(gdown::Float64, gup::Float64, gdu::Float64, ksoc::Float64, Omega::Float64, Msize0::Int64, Np::Int64)

    # Msize = Msize0*2
    # indvec = cutMsizeEne(Msize0,Np,Ene0minumhalf)
    # Msizecut = length(indvec)

    mat0 = Hsocfunc(Msize0,Np,ksoc,Omega)
    mat1, mat2, mat3 = Hintfunc(Msize0,Np)
    lambda, phi = diagonaliseHsoc(mat0+gdown*mat1+gup*mat2+gdu*mat3)

    return lambda

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 1, 100]
# operator
ntimesperdecade = 10
nfreqsperdecade = 5
nkᵣeval = 60
zTx = -120
zRx = -80
x_rx = -115.
y_rx = 0.
rx_roll = 10. #ACHTUNG!
rx_pitch = 10. #ACHTUNG!
rx_yaw = 10. #ACHTUNG!
tx_roll = 0.
tx_pitch = 0.
tx_yaw = 0.
addprimary = true
tempest = transD_GP.TEMPEST1DInversion.Bfield(
  					ntimesperdecade = ntimesperdecade,
            nfreqsperdecade = nfreqsperdecade,
					  nkᵣeval = nkᵣeval,
            times  = times,
            ramp   = ramp,
					  zTx    = zTx,
					  zRx    = zRx,
					  x_rx   = x_rx,
					  y_rx   = y_rx,
					  rx_roll = rx_roll,
					  rx_pitch = rx_pitch,
					  rx_yaw = rx_yaw,
					  tx_roll = tx_roll,
					  tx_pitch = tx_pitch,
					  tx_yaw = tx_yaw,
					  addprimary=addprimary)
# X and Z components in fT from GA-AEM
ross_primary = [28.5962      -13.689]
ross_secondary = [
   8.03371     -6.620856
  5.849011     -5.572887
  4.957888     -5.037114
  4.362005     -4.626088
  3.907595     -4.288714
  3.510647     -3.988984
  3.048846      -3.63051
  2.496073     -3.168608
  1.889611     -2.605418
  1.286286     -1.967748
0.7633423     -1.324871
0.3876874    -0.7776925
0.1682725    -0.3943537
0.06388109    -0.1749779
0.02070349   -0.06618546]

if addprimary
	ross_secondary .+= ross_primary
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using Replay
using Test

const repl_script = """
x = 1 + 1
println("Hello World")
@show x
"""

@testset "replay: color=yes" begin
    color = "yes"
    buf = replay(repl_script, IOBuffer(); cmd="-q --color=yes")
    out = buf |> take! |> String
    #=
    open(joinpath(@__DIR__, "references", "replay_color_$(color)_julia.txt"), "w") do f
        write(f, out)
    end
    =#
    reftxt = joinpath(@__DIR__, "references", "replay_color_yes_julia.txt")
    ref = join(readlines(reftxt), "\r\n")
    @test out == ref
end

@testset "replay: color=no" begin
    color = "no"
    buf = replay(repl_script, IOBuffer(); cmd="-q --color=no")
    out = buf |> take! |> String
    #=
    open(joinpath(@__DIR__, "references", "replay_color_$(color)_julia.txt"), "w") do f
        write(f, out)
    end
    =#
    reftxt = joinpath(@__DIR__, "references", "replay_color_no_julia.txt")
    ref = join(readlines(reftxt), "\r\n")
    @test out == ref
end

#=
using Replay, Test
endexamples_dir = joinpath(pkgdir(Replay), "examples")
examples_dir = joinpath(pkgdir(Replay), "examples")
for example in readdir(examples_dir)
    example == "disable_color" && continue
    @testset "$example" begin
        apppath = joinpath(examples_dir, example, "app.jl")
        julia_exepath = joinpath(Sys.BINDIR::String, Base.julia_exename())
        @info "running $(example)"
        r = run(`$(julia_exepath) $(apppath)`)
        @test r.exitcode == 0
    end
end
=#


text:
using BinDeps

@BinDeps.setup

liblsoda = library_dependency("liblsoda")
BinDeps.provides(Sources,URI("https://github.com/sdwfrost/liblsoda/archive/v0.1.1.tar.gz"),liblsoda,unpacked_dir="liblsoda-0.1.1")

lsodadir = BinDeps.depsdir(liblsoda)

libdir=joinpath(lsodadir,"usr","lib")

srcdir = joinpath(lsodadir,"src","liblsoda-0.1.1")
builddir = joinpath(srcdir,"src")

provides(Binaries,
    URI("https://dl.bintray.com/sdwfrost/generic/liblsoda.7z"),
    [liblsoda], unpacked_dir="bin$(Sys.WORD_SIZE)",
    os = :Windows)

if Sys.KERNEL == :Darwin
    suffix="dylib"
elseif Sys.KERNEL == :NT
    suffix="dll"
else
    suffix="so"
end

BinDeps.provides(SimpleBuild,
    (@build_steps begin
        BinDeps.GetSources(liblsoda)
	BinDeps.CreateDirectory(libdir)
        @build_steps begin
            BinDeps.ChangeDirectory(srcdir)
            BinDeps.MakeTargets()
        `cp $builddir/liblsoda.$suffix $libdir`
        end
end), [liblsoda], os = :Unix)
@BinDeps.install Dict(:liblsoda => :liblsoda)


text:
## IOStream

const sizeof_ios_t = int(ccall(:jl_sizeof_ios_t, Int32, ()))

type IOStream <: IO
    handle::Ptr{Void}
    ios::Array{Uint8,1}
    name::String
    mark::Int64

    IOStream(name::String, buf::Array{Uint8,1}) = new(pointer(buf), buf, name, -1)
end
# TODO: delay adding finalizer, e.g. for memio with a small buffer, or
# in the case where we takebuf it.
function IOStream(name::String, finalize::Bool)
    buf = zeros(Uint8,sizeof_ios_t)
    x = IOStream(name, buf)
    if finalize
        finalizer(x, close)
    end
    return x
end
IOStream(name::String) = IOStream(name, true)

convert(T::Type{Ptr{Void}}, s::IOStream) = convert(T, s.ios)
show(io::IO, s::IOStream) = print(io, "IOStream(", s.name, ")")
fd(s::IOStream) = int(ccall(:jl_ios_fd, Clong, (Ptr{Void},), s.ios))
stat(s::IOStream) = stat(fd(s))
close(s::IOStream) = ccall(:ios_close, Void, (Ptr{Void},), s.ios)
isopen(s::IOStream) = bool(ccall(:ios_isopen, Cint, (Ptr{Void},), s.ios))
function flush(s::IOStream)
    sigatomic_begin()
    systemerror("flush", ccall(:ios_flush, Cint, (Ptr{Void},), s.ios) != 0)
    sigatomic_end()
    s
end
iswritable(s::IOStream) = bool(ccall(:ios_get_writable, Cint, (Ptr{Void},), s.ios))
isreadable(s::IOStream) = bool(ccall(:ios_get_readable, Cint, (Ptr{Void},), s.ios))
modestr(s::IO) = modestr(isreadable(s), iswritable(s))
modestr(r::Bool, w::Bool) = r ? (w ? "r+" : "r") : (w ? "w" : error("Neither readable nor writable"))

function truncate(s::IOStream, n::Integer)
    systemerror("truncate", ccall(:ios_trunc, Int32, (Ptr{Void}, Uint), s.ios, n) != 0)
    return s
end

function seek(s::IOStream, n::Integer)
    ret = ccall(:ios_seek, FileOffset, (Ptr{Void}, FileOffset), s.ios, n)
    systemerror("seek", ret == -1)
    ret < -1 && error("seek failed")
    return s
end

seekstart(s::IO) = seek(s,0)

function seekend(s::IOStream)
    systemerror("seekend", ccall(:ios_seek_end, FileOffset, (Ptr{Void},), s.ios) != 0)
    return s
end

function skip(s::IOStream, delta::Integer)
    ret = ccall(:ios_skip, FileOffset, (Ptr{Void}, FileOffset), s.ios, delta)
    systemerror("skip", ret == -1)
    ret < -1 && error("skip failed")
    return s
end

function position(s::IOStream)
    pos = ccall(:ios_pos, FileOffset, (Ptr{Void},), s.ios)
    systemerror("position", pos == -1)
    return pos
end

eof(s::IOStream) = bool(ccall(:jl_ios_eof, Int32, (Ptr{Void},), s.ios))

# For interfacing with C FILE* functions

immutable CFILE
    ptr::Ptr{Void}
end

function CFILE(s::IO)
    @unix_only FILEp = ccall(:fdopen, Ptr{Void}, (Cint, Ptr{Uint8}), convert(Cint, fd(s)), modestr(s))
    @windows_only FILEp = ccall(:_fdopen, Ptr{Void}, (Cint, Ptr{Uint8}), convert(Cint, fd(s)), modestr(s))
    systemerror("fdopen", FILEp == C_NULL)
    seek(CFILE(FILEp), position(s))
end

convert(::Type{CFILE}, s::IO) = CFILE(s)

function seek(h::CFILE, offset::Integer)
    systemerror("fseek", ccall(:fseek, Cint, (Ptr{Void}, Clong, Cint),
                               h.ptr, convert(Clong, offset), int32(0)) != 0)
    h
end

position(h::CFILE) = ccall(:ftell, Clong, (Ptr{Void},), h.ptr)

## constructing and opening streams ##

# "own" means the descriptor will be closed with the IOStream
function fdio(name::String, fd::Integer, own::Bool=false)
    s = IOStream(name)
    ccall(:ios_fd, Ptr{Void}, (Ptr{Void}, Clong, Int32, Int32),
          s.ios, fd, 0, own);
    return s
end
fdio(fd::Integer, own::Bool=false) = fdio(string("<fd ",fd,">"), fd, own)

function open(fname::String, rd::Bool, wr::Bool, cr::Bool, tr::Bool, ff::Bool)
    s = IOStream(string("<file ",fname,">"))
    systemerror("opening file $fname",
                ccall(:ios_file, Ptr{Void},
                      (Ptr{Uint8}, Ptr{Uint8}, Int32, Int32, Int32, Int32),
                      s.ios, fname, rd, wr, cr, tr) == C_NULL)
    if ff
        systemerror("seeking to end of file $fname", ccall(:ios_seek_end, FileOffset, (Ptr{Void},), s.ios) != 0)
    end
    return s
end
open(fname::String) = open(fname, true, false, false, false, false)

function open(fname::String, mode::String)
    mode == "r"  ? open(fname, true , false, false, false, false) :
    mode == "r+" ? open(fname, true , true , false, false, false) :
    mode == "w"  ? open(fname, false, true , true , true , false) :
    mode == "w+" ? open(fname, true , true , true , true , false) :
    mode == "a"  ? open(fname, false, true , true , false, true ) :
    mode == "a+" ? open(fname, true , true , true , false, true ) :
    error("invalid open mode: ", mode)
end

function open(f::Function, args...)
    io = open(args...)
    try
        f(io)
    finally
        close(io)
    end
end

## low-level calls ##

write(s::IOStream, b::Uint8) = int(ccall(:ios_putc, Int32, (Uint8, Ptr{Void}), b, s.ios))

function write{T}(s::IOStream, a::Array{T})
    if isbits(T)
        if !iswritable(s)
            error("attempt to write to a read-only IOStream")
        end
        int(ccall(:ios_write, Uint, (Ptr{Void}, Ptr{Void}, Uint),
                  s.ios, a, length(a)*sizeof(T)))
    else
        invoke(write, (IO, Array), s, a)
    end
end

function write(s::IOStream, p::Ptr, nb::Integer)
    if !iswritable(s)
        error("attempt to write to a read-only IOStream")
    end
    int(ccall(:ios_write, Uint, (Ptr{Void}, Ptr{Void}, Uint), s.ios, p, nb))
end

function write{T,N,A<:Array}(s::IOStream, a::SubArray{T,N,A})
    if !isbits(T) || stride(a,1)!=1
        return invoke(write, (Any, AbstractArray), s, a)
    end
    colsz = size(a,1)*sizeof(T)
    if N<=1
        return write(s, pointer(a, 1), colsz)
    else
        cartesianmap((idxs...)->write(s, pointer(a, idxs), colsz),
                     tuple(1, size(a)[2:end]...))
        return colsz*trailingsize(a,2)
    end
end

# num bytes available without blocking
nb_available(s::IOStream) = ccall(:jl_nb_available, Int32, (Ptr{Void},), s.ios)

function read(s::IOStream, ::Type{Uint8})
    b = ccall(:ios_getc, Int32, (Ptr{Void},), s.ios)
    if b == -1
        throw(EOFError())
    end
    uint8(b)
end

function read!{T}(s::IOStream, a::Array{T})
    if isbits(T)
        nb = length(a)*sizeof(T)
        if ccall(:ios_readall, Uint,
                 (Ptr{Void}, Ptr{Void}, Uint), s.ios, a, nb) < nb
            throw(EOFError())
        end
    else
        invoke(read!, (IO, Array), s, a)
    end
    a
end

## text I/O ##

function write(s::IOStream, c::Char)
    if !iswritable(s)
        error("attempt to write to a read-only IOStream")
    end
    int(ccall(:ios_pututf8, Int32, (Ptr{Void}, Char), s.ios, c))
end
read(s::IOStream, ::Type{Char}) = ccall(:jl_getutf8, Char, (Ptr{Void},), s.ios)

takebuf_string(s::IOStream) =
    ccall(:jl_takebuf_string, ByteString, (Ptr{Void},), s.ios)

takebuf_array(s::IOStream) =
    ccall(:jl_takebuf_array, Vector{Uint8}, (Ptr{Void},), s.ios)

function takebuf_raw(s::IOStream)
    sz = position(s)
    buf = ccall(:jl_takebuf_raw, Ptr{Uint8}, (Ptr{Void},), s.ios)
    return buf, sz
end

function sprint(size::Integer, f::Function, args...)
    s = IOBuffer(Array(Uint8,size), true, true)
    truncate(s,0)
    f(s, args...)
    takebuf_string(s)
end

sprint(f::Function, args...) = sprint(0, f, args...)

write(x) = write(STDOUT::IO, x)

function readuntil(s::IOStream, delim::Uint8)
    ccall(:jl_readuntil, Array{Uint8,1}, (Ptr{Void}, Uint8), s.ios, delim)
end

function readbytes!(s::IOStream, b::Array{Uint8}, nb=length(b))
    olb = lb = length(b)
    nr = 0
    while !eof(s) && nr < nb
        if lb < nr+1
            lb = max(65536, (nr+1) * 2)
            resize!(b, lb)
        end
        nr += int(ccall(:ios_readall, Uint,
                        (Ptr{Void}, Ptr{Void}, Uint),
                        s.ios, pointer(b, nr+1), min(lb-nr, nb-nr)))
    end
    if lb > olb
        resize!(b, nr) # shrink to just contain input data if was resized
    end
    return nr
end

readbytes(s::IOStream) = readbytes(s, filesize(s))

## Character streams ##
const _chtmp = Array(Char, 1)
function peekchar(s::IOStream)
    if ccall(:ios_peekutf8, Int32, (Ptr{Void}, Ptr{Char}), s, _chtmp) < 0
        return char(-1)
    end
    return _chtmp[1]
end

function peek(s::IOStream)
    ccall(:ios_peekc, Int32, (Ptr{Void},), s)
end

function skipchars(s::IOStream, pred; linecomment::Char=char(0xffffffff))
    ch = peekchar(s); status = int(ch)
    while status >= 0 && (pred(ch) || ch == linecomment)
        if ch == linecomment
            readline(s)
        else
            read(s, Char)  # advance one character
        end
        ch = peekchar(s); status = int(ch)
    end
    return s
end

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
function partition_interval(nbins::Integer, minval::Real, maxval::Real)
    edges = (range(0,step=(maxval - minval) / nbins,length=nbins)) .+ minval
end

function partition_interval(nbins::Integer, minval::AbstractGray, maxval::AbstractGray)
    partition_interval(nbins, minval.val, maxval.val)
end

"""
```
edges, count = build_histogram(img, nbins)
edges, count = build_histogram(img, nbins; minval, maxval)
edges, count = build_histogram(img, edges)
```

Generates a histogram for the image over `nbins` spread between `[minval, maxval]`.
Color images are automatically converted to grayscale.

# Output
Returns `edges` which is a `AbstractRange` type that specifies how the  interval
`[minval, maxval]` is divided into bins, and an array `count` which records the
concomitant bin frequencies. In particular, `count` has the following
properties:

* `count[0]` is the number satisfying `x < edges[1]`
* `count[i]` is the number of values `x` that satisfy `edges[i] <= x < edges[i+1]`
* `count[end]` is the number satisfying `x >= edges[end]`.
* `length(count) == length(edges)+1`.

# Details

One can consider a histogram as a piecewise-constant model of a probability
density function ``f`` [1]. Suppose that ``f`` has support on some interval ``I =
[a,b]``.  Let ``m`` be an integer and ``a = a_1 < a_2 < \\ldots < a_m < a_{m+1} =
b`` a sequence of real numbers. Construct a sequence of intervals

```math
I_1 = [a_1,a_2], I_2 = (a_2, a_3], \\ldots, I_{m} = (a_m,a_{m+1}]
```

which partition ``I`` into subsets ``I_j`` ``(j = 1, \\ldots, m)`` on which
``f`` is constant. These subsets satisfy ``I_i \\cap I_j = \\emptyset, \\forall
i \\neq j``, and are commonly referred to as *bins*. Together they encompass the
entire range of data values such that ``\\sum_j |I_j | = | I |``. Each bin has
width ``w_j = |I_j| = a_{j+1} - a_j`` and height ``h_j`` which is the constant
probability density over the region of the bin. Integrating the constant
probability density over the width of the bin ``w_j`` yields a probability mass
of ``\\pi_j = h_j w_j`` for the bin.

For a sample ``x_1, x_2, \\ldots, x_N``, let

```math
n_j = \\sum_{n = 1}^{N}\\mathbf{1}_{(I_j)}(x_n),
\\quad \\text{where} \\quad
\\mathbf{1}_{(I_j)}(x) =
\\begin{cases}
 1 & \\text{if} \\; x \\in I_j,\\\\
 0 & \\text{otherwise},
\\end{cases},
```
represent the number of samples falling into the interval ``I_j``. An estimate
for the probability mass of the ``j``th bin is given by the relative frequency
``\\hat{\\pi} = \\frac{n_j}{N}``, and the histogram estimator of the probability
density function is defined as
```math
\\begin{aligned}
\\hat{f}_n(x)  & = \\sum_{j = 1}^{m}\\frac{n_j}{Nw_j} \\mathbf{1}_{(I_j)}(x) \\\\
& = \\sum_{j = 1}^{m}\\frac{\\hat{\\pi}_j}{w_j} \\mathbf{1}_{(I_j)}(x) \\\\
& = \\sum_{j = 1}^{m}\\hat{h}_j \\mathbf{1}_{(I_j)}(x).
\\end{aligned}
```

The function ``\\hat{f}_n(x)`` is a genuine density estimator because ``\\hat{f}_n(x)  \\ge 0`` and
```math
\\begin{aligned}
\\int_{-\\infty}^{\\infty}\\hat{f}_n(x) \\operatorname{d}x & = \\sum_{j=1}^{m} \\frac{n_j}{Nw_j} w_j \\\\
& = 1.
\\end{aligned}
```

# Options
Various options for the parameters of this function are described in more detail
below.

## Choices for `nbins`
You can specify the number of discrete bins for the histogram. When specifying
the number of bins consider the maximum number of graylevels that your image
type supports. For example, with an image of type `N0f8` there is a maximum
of 256 possible graylevels. Hence, if you request more than 256 bins for
that type of image you should expect to obtain zero counts for numerous bins.

## Choices for `minval`
You have the option to specify the lower bound of the interval over which the
histogram will be computed.  If `minval` is not specified then the minimum
value present in the image is taken as the lower bound.

## Choices for `maxval`
You have the option to specify the upper bound of the interval over which the
histogram will be computed.  If `maxval` is not specified then the maximum
value present in the image is taken as the upper bound.

## Choices for `edges`
If you do not designate the number of bins, nor the lower or upper bound of the
interval, then you have the option to directly stipulate how the intervals will
be divided by specifying a `AbstractRange` type.

# Example

Compute the histogram of a grayscale image.
```julia

using TestImages, FileIO, ImageView

img =  testimage("mandril_gray");
edges, counts  = build_histogram(img, 256, minval = 0, maxval = 1)
```

Given a color image, compute the histogram of the red channel.
```julia
img = testimage("mandrill")
r = red.(img)
edges, counts  = build_histogram(r, 256, minval = 0, maxval = 1)
```

# References
[1] E. Herrholz, "Parsimonious Histograms," Ph.D. dissertation, Inst. of Math. and Comp. Sci., University of Greifswald, Greifswald, Germany, 2011.
"""
function build_histogram(img::AbstractArray, nbins::Integer = 256; minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray})
    edges = partition_interval(nbins, minval, maxval)
    build_histogram(img, edges)
end

# Warning is due to: https://github.com/JuliaLang/julia/issues/9498
# Might resolve this based on this idea: https://discourse.julialang.org/t/fun-with-kwargs-methods/15711/8
function build_histogram(img::AbstractArray, nbins::Integer = 256)
    build_histogram(img, nbins, minval =  minfinite(img), maxval = maxfinite(img))
end

function build_histogram(img::AbstractArray{T}, edges::AbstractRange) where {T<:Color3}
    build_histogram(Gray.(img),edges)
end

function build_histogram(img::AbstractArray, edges::AbstractRange)
    Base.has_offset_axes(edges) && throw(ArgumentError("edges must be indexed starting with 1"))
    lb = first(axes(edges,1))-1
    ub = last(axes(edges,1))
    first_edge = first(edges)
    step_size = step(edges)
    counts = fill(0, lb:ub)
    for val in img
         if isnan(val)
             continue
         else
            if val >= edges[end]
                counts[end] += 1
                continue
            elseif val < first(edges)
                counts[lb] += 1
            else
                index = Int(Base.div(val-first_edge,step_size)) + 1
                counts[index] += 1
            end
        end
    end
    edges, counts
end


text:
using Boot

using TestSetExtensions

using Base.Test

using Suppressor

@testset ExtendedTestSet "All Tests" begin

  include("core/index.jl")

end


text:
using RedPitayaDAQServer, Documenter

makedocs(
    format = Documenter.HTML(prettyurls = false),
    modules = [RedPitayaDAQServer],
    sitename = "RP DAQ Server",
    authors = "Tobias Knopp, Jonas Beuke, Matthias Gräser",
    pages = [
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Architecture" => "architecture.md",
        "Connections" => "connections.md",
        "Cluster" => "cluster.md",
        "Client" => "client.md",
        "Examples" => Any["Simple" => "examples/simple.md",
                          "SlowADC" => "examples/slowADC.md"],
        "SCPI Interface" => "scpi.md",
        "FPGA Development" => "fpga.md",
        #"Getting Started" => "overview.md",
    ]
#    html_prettyurls = false, #!("local" in ARGS),
)

deploydocs(repo   = "github.com/tknopp/RedPitayaDAQServer.jl.git",
           target = "build")

max_stars_repo_name: fredrikekre/RedPitayaDAQServer

max_issues_repo_name: fredrikekre/RedPitayaDAQServer

max_forks_repo_name: fredrikekre/RedPitayaDAQServer