text:
#import Base.==
using AutoHashEquals
@auto_hash_equals mutable struct swahiliTag
    word::String
    pos::String
    msd::String
    behind::String
    ahead::String
    cls::String
    tam::String
    subj::String
    obj::String
    rel::String
    gloss::String
    syntax::String
    lemma::String
end

"""
Method for measuring equality between swahili tag objects
"""
#function ==(tag1::swahiliTag, tag2::swahiliTag)
#    if tag1.word == tag2.word &&
#        tag1.pos == tag2.pos &&
#        tag1.msd == tag2.msd &&
#        tag1.behind == tag2.behind &&
#        tag1.ahead == tag2.ahead &&
#        tag1.cls == tag2.ahead &&
#        tag1.tam == tag2.tam &&
#        tag1.subj == tag2.subj &&
#        tag1.obj == tag2.obj &&
#        tag1.rel == tag2.rel &&
#        tag1.gloss == tag2.gloss &&
#        tag1.syntax == tag2.syntax &&
#        tag1.lemma == tag2.lemma
#        return true
#    else
#        return false
#    end
#end
#function ==(tagList1::Array{swahiliTag, 1}, tagList2::Array{swahiliTag, 1})
#    if length(tagList1) != length(tagList2)
#        return false
#    else
#        for tagI = 1:length(tagList1)
#            if tagList1[tagI] != tagList2[tagI]
#                return false
#            end
#        end
#        return true
#    end
#end
#function isequal(tag1::swahiliTag, tag2::swahiliTag)
#   return tag1 == tag2
#end

#Outer constructor for swahiliTag to provide default arguments
swahiliTag(word,
            pos,
            msd;
            behind="-",
            ahead="-",
            cls="-",
            tam="-",
            subj="-",
            obj="-",
            rel="-",
            gloss="-",
            syntax="-",
            lemma="-") =
        swahiliTag(word, pos, msd,
                        behind, ahead,
                        cls, tam, subj,
                        obj, rel, gloss,
                        syntax,lemma)
mutable struct swahiliQuery
    jsonFileName::String
    words::Array{swahiliTag,1}
    data::Dict{String, Any}
    uniques::Array{String, 1}
end

#Outer constructor for swahiliQuery to provide default arguments
swahiliQuery(jsonFileName,
            words=Array{swahiliTag,1}(),
            data=Dict{String, Any}(),
            uniques=Array{String,1}()) =
        swahiliQuery(jsonFileName, words, data, uniques)

#type for holding sentence information taken from the helsinki corpus

"""
This class holds the information for a single sentence in an interlinear text
This is used to extract the tag sequences for nltk's hmm tagger

There are no default arguments for this because both attributes are things that cannot be assigned defaults
"""
mutable struct sentence
    words::Array{swahiliTag, 1}
    gloss::String
end

#Additional swahliQuery type that

struct dependency
    head::Int
    dep::Int
    label::String
end
# outer constructor for dependency
dependency(head, dep, label="-") = dependency(head, dep, label)

mutable struct dependencySent
    tokenList::Array{swahiliTag,1}
    depList::Array{dependency, 1}
end

function getDepParse(sent::dependencySent)
    tokenList = sent.tokenList
    deps = sent.dependency
    for tokenI = 1:length(tokenList)
        token = tokenList[tokenI]
        if token.syntax == "@<DN"
            unifyTag = swahiliTag("-", "N", "-")
            nounID = seekTagIndex(sent, tokenI, unifyingTag=unifyTag, direction="Aft")
            push!(deps, dependency(nounID, tokenI))
        end
    end
end
function seekTagIndex(sent::Array{swahiliTag,1}, startIndex; unifyTag::swahiliTag=swahiliTag(), direction::String="None")
    i = startIndex
    found = false
    if direction != "None"
        while i >= 1 && i <= length(sent)
            if direction == "Fore"
                i += 1
            else
                i -= 1
            end
        end
    else
        while i + startIndex <= length(sent) && i - startIndex >= 1

        end
    end
end
#function isUnifiable()
#
#end
function getCoNLLU(sent::dependencySent)
    conll = ""
    for tokI = 1:length(sent.tokenList)
        token = sent.tokenList[depI]
        head, label = [(arc.head, arc.label) for arc in sent.depList if arc.dep == tokI]
        form = token.word
        lemma = token.lemma
        upos = getUPosTag(token.pos)
        xpos = token.pos

        # The structure of this is
        # ID FORM LEMMA UPOSTAG XPOSTAG FEATS HEAD DEPREL DEPS MISC
        line = string(tokI) * token.word * token.lemma * token.getUPosTag(token.pos)
    end
    return conll
end

"""
This function does translation of tags from the Helsinki Corpus of Swahili
to Universal Part of Speech tags
"""
function getUPosTag(tag::swahiliTag)

    tagCorr = Dict(
                    "ABBR"      => "SYM",   #Probably a good idea to use the syntactic tag for these
                                            #All of the examples I looked at were ADVs though
                    "ADJ"       => "ADJ",
                    "ADV"       => "ADV",
                    "AG-PART"   => "ADP",   #This refers to the  use of 'na' to represent oblique
                                            #subjects in passives
                    "CC"        => "CCONJ",
                    "CONJ"      => "SCONJ", #These are all @CS in the data I looked at, it appears
                                            # that CC is reserved for coordinating conjunctions and
                                            #Everything else falls under here
                    "CONJ/CC"   => "CCONJ",
                    "DEM"       => "DET",
                    "EXCLAM"    => "INTJ",  #need an exception here because ! is also EXLCLAM
                    "GEN-CON"   => "ADP",   #perhaps a feature is needed to differentiate these?
                    "GEN-CON"   => "ADP",
                    "INTERROG"  => "PRON",
                    "N"         => "NOUN",
                    "NUM"       => "NUM",
                    "NUM-ROM"   => "NUM",
                    "POSS-PRON" => "SCONJ", #These are the mwenye, penye, zenye style relative markers
                                            #meaning something like 'which has'
                    "PREP"      => "ADP",
                    "PREP/ADV"  => "ADV",   #These are actually a false compund
                                            # The hcs counts these as a single word but they have a space intervening
                    "PRON"      => "PRON",
                    "PROPNAME"  => "PROPN",
                    "REL-LI"    => "PRON",  #UD treats inflected relative clauses as pronouns
                    "REL-LI-VYO"=> "PRON",  # Not subordinating conjunctions
                    "REL-SI"    => "PRON",
                    "REL-SI-VYO"=> "PRON",
                    "TITLE"     => "PROPN",
                    "V"         => "VERB",
                    "V-BE"      => "AUX",
                    "V-DEF"     => "VERB",  # These so-called 'defective verbs' are almost all
                                            # verbs of arabic origin that have frozen, uninflected forms
                    )
    puncts = [  "COLON",
                "COMMA",
                "DOUBLE-QUOTE",
                "DOUBLE-QUOTE-CLOSING",
                "DOUBLE-QUOTE-OPENING",
                "HYPHEN",
                "LEFT-PARENTHESIS",
                "PERCENT-MARK",
                "QUESTION-MARK",
                "RIGHT-PARENTHESIS",
                "SEMI-COLON",
                "SINGLE-QUOTE",
                "SINGLE-QUOTE-CLOSING",
                "SINGLE-QUOTE-OPENING",
                "SLASH",
                "STOP"]
    #Catch !
    if tag.word == "!" || tag.pos in puncts
        return "PUNCT"
    elseif contains(tag.pos, "REL-LI")
        if contains(tag.syntax, "V")    #Some of these are functioning as verbs and some are functioning as nouns in sentences
                                        #the only clue is in the syntax piece of the word
            return "VERB"
        else
            return "PRON"
        end
    else
        return tagCorr[tag.pos]
    end
end
#function getUDFeats(tag::swahiliTag)
#    """
#    This function uses the information in the tag to determine what the appropriate
#    features using the UD feature inventory are.
#    """
#    break
#end

"""
This function removes the msd tags that specify individual words, these are too easy to guess
"""
function stripMSD!(tag::swahiliTag)
    tag.msd = replace(tag.msd, r"\[.*\]", "")
    tag.msd = strip(replace(tag.msd, r"\:.*[ ]?", " "))
end

"""
This function sets the cls attribute of a tag if it is a noun, if not it is left in the default
"""
function setCLS!(tag::swahiliTag)
    if tag.pos == "N"
        classSearch = eachmatch(r"\s*([0-9/]+)-((?:PLSG)|(?:SG)|(?:PL))\s*", tag.msd)
        if !isempty(classSearch)
            classRes = start(classSearch).match
            number = classRes
            if classRes[2] == "SG"
                if contains(number, "/")
                    tag.cls = split(number, "/")[1]
                else
                    tag.cls = number
                end
            elseif classRes[2] == "PL"
                if contains(number, "/")
                    tag.cls = split(number, "/")[2]
                else
                    tag.cls = number
                end
            else
                tag.cls = number
            end
        end
    end
end
"""
This function sets the tam, subj, obj and rel attributes of a tag if it is a verb, if not these are left alone
(For non-verbs, these tags are not relevant)
"""
function setVerbTags!(tag::swahiliTag)
    if tag.pos == "V"
        classSearch = eachmatch(r"TAM=([A-Z]+)(?:-([A-Z]+)){0,1}(?:\:([a-z]+)){0,1}", tag.msd)
        if !isempty(classSearch)
            tag.tam = start(classSearch).match
        end
        haystackPieces = split(tag.msd)
        for piece in haystackPieces
            index = search(piece, r"SUB-PREF=.*")
            if index != 0:-1
                tag.subj = piece[index]
            end
            index = search(piece, r"OBJ-PREF=.*")
            if index != 0:-1
                tag.obj = piece[index]
            end
            index = search(piece, r"REL-PREF=.*")
            if index != 0:-1
                tag.rel = piece[index]
            end
        end
    end
end

function setWordContext!(tag::swahiliTag; behind=nothing, ahead=nothing, prefix=true)
    if prefix == true
        if behind != nothing
            if prefix==true && length(behind) > 2
                tag.behind = lowercase(behind)[1:2]
            else
                tag.behind = lowercase(behind)
            end
        else
            tag.behind = "nothingBOS"
        end
        if ahead != nothing
            if prefix == true && length(ahead) > 2
                tag.ahead = lowercase(ahead)[1:2]
            else
                tag.ahead = lowercase(ahead)
            end
        else
            tag.ahead = "nothingEOS"
        end
    end
end

function getCharNGrams(tag::swahiliTag, n::Int)
    grams = []
    i = 1
    #Add buffer
    word = " " * tag.word * " "^n
    for i in eachindex(word[1:end-n])
        push!(grams, word[i:i+n-1])
    end
    return grams
end

"""
    getFeatureValues(swahiliTag, uniques)

uniques must be an array of strings. Fill in a dictionary with the count or simple boolean values

"""
function getFeatureValues(tag::swahiliTag, uniques::Array{String}, counts=true)
    #initialize the dictionary
    initData = zeros(Int, length(uniques))
    featureDict = Dict(zip(uniques, initData))
    for n in [1,2,3,4]
        for gram in getCharNGrams(tag, n)
            try
                featureDict[gram] = counts ? featureDict[gram] + 1 : 1
            end
        end
    end
    if word.behind != "nothingBOS"
        try
            featureDict[word.behind] = counts ? featureDict[word.behind] + 1 : 1
        end
    end
    if word.ahead != "nothingEOS"
        try
            featureDict[word.ahead] = counts ? featureDict[word.ahead] + 1 : 1
        end
    end
    return featureDict
end
"""
    getFeatureDict(swahiliTag)

Fill in a dictionary with the count or simple boolean values

"""
function getFeatureValues(word::swahiliTag, counts=true)
    #initialize the dictionary
    featureDict = Dict{String, Int}()
    for n in [1,2,3,4]
        for gram in getCharNGrams(word, n)
            try
                featureDict[gram] = counts ? featureDict[gram] + 1 : 1
            catch err
                if isa(err, KeyError)
                    featureDict[gram] = 1
                end
            end
        end
    end
    if word.behind != "nothingBOS"
        try
            featureDict[word.behind] = counts ? featureDict[word.behind] + 1 : 1
        catch err
            if isa(err, KeyError)
                featureDict[word.behind] = 1
            end
        end
    end
    if word.ahead != "nothingEOS"
        try
            featureDict[word.ahead] = counts ? featureDict[word.ahead] + 1 : 1
        catch err
            if isa(err, KeyError)
                featureDict[word.ahead] = 1
            end
        end
    end
    return featureDict
end

"""
    fillQuery!(query)

Fill the query's words and unique attributes using the json file in
`swahiliQuery.jsonFileName`

cases where the filename is actually a directory are handled appropriately.
"""
function fillQuery!(query::swahiliQuery)
    uniques = []
    sents = []
    if isdir(query.jsonFileName)
        files = readdir(query.jsonFileName)
        res = @sync @parallel vcat for file in files
            tempUniques = []
            println("Processing file " * file)
            words = getWordsFromJSON(query.jsonFileName * file)
            sents = vcat(sents, words)
            for word in words
                append!(tempUniques, keys(getFeatureValues(word)))
            end
            (tempUniques, sents)
        end
        tempUniques, tempWords = collect(zip(res...))
        #query.uniques = unique(uniques)
        query.uniques = unique(tempUniques[1])#[result[1] for result in res]
        query.words = tempWords[1]#[result[2] for result in res]
        println(length(query.uniques))
        println(length(query.words))
    else
        query.words = getWordsFromJSON(query.jsonFileName)
        uniques = []
        for word in query.words
            append!(uniques,keys(getFeatureValues(word)))
        end
        query.uniques = unique(uniques)
    end
end

"""
    getWordsFromJSON(filepath)

Generate a list of swahiliTag objects for the words that occur in the specified json filepath

"""
function getWordsFromJSON(filepath::String, pullUniques=true)
    data = JSON.parsefile(filepath)["kwic"]
    words = Array{swahiliTag, 1}()
    for hit in data
        hitSize = length(hit["tokens"])
        for tokenI = 1:hitSize
            token = hit["tokens"][tokenI]
            if true #token["pos"] in ["V","N"] Uncomment this to filter pos tags
                #try
                    token = hit["tokens"][tokenI]
                    for key in keys(token)
                        if token[key] == nothing
                            token[key] = "-"
                        end
                    end

                    currentTag = swahiliTag(token["word"],
                                                token["pos"],
                                                token["msd"],
                                                gloss=token["gloss"],
                                                syntax=token["syntax"]
                                                )
                    previousToken = tokenI != 1 ? hit["tokens"][tokenI - 1]["word"] : nothing
                    nextToken = tokenI != hitSize ? hit["tokens"][tokenI + 1]["word"] : nothing
                    stripMSD!(currentTag)
                    setCLS!(currentTag)
                    setVerbTags!(currentTag)
                    setWordContext!(currentTag, behind=previousToken, ahead=nextToken)
                    push!(words, currentTag)
                #catch a
                #    print(a)
                #    println("token")
                #    println(token["word"])
                #    println(token["pos"])
                #    println(token["msd"])
                #    println(token["gloss"])
                #    println(token["syntax"])
                #end
            end
        end
    end
    return words
end
"""
    getSentsFromJSON(filepath)

Generate a list of sentences.
Each sentence is itself a list of swahiliTag objects for the words that occur in the specified json filepath

"""
function getSentsFromJSON(filepath::String)
    data = JSON.parsefile(filepath)["kwic"]
    sents = Array{Array{swahiliTag,1},1}()
    for hit in data
        words = Array{swahiliTag, 1}()
        hitSize = length(hit["tokens"])
        for tokenI = 1:hitSize
            token = hit["tokens"][tokenI]
            try
                token = hit["tokens"][tokenI]
                for key in keys(token)
                    if token[key] == nothing
                        token[key] = "-"
                    end
                end
                currentTag = swahiliTag(token["word"],
                                            token["pos"],
                                            token["msd"],
                                            gloss=token["gloss"],
                                            syntax=token["syntax"]
                                            )

                previousToken = tokenI != 1 ? hit["tokens"][tokenI - 1]["word"] : nothing
                nextToken = tokenI != hitSize ? hit["tokens"][tokenI + 1]["word"] : nothing
                stripMSD!(currentTag)
                setCLS!(currentTag)
                setVerbTags!(currentTag)
                push!(words, currentTag)
            catch err
                println("ERROR" * err.msg)
                println("token")
                println(token["word"])
                println(token["msd"])
                println(token["gloss"])
                println(token["syntax"])
                exit()
            end
        end
        push!(sents, words)
    end
    return sents
end

"""
This function converts sents into a flattened swahiliQuery objects
"""
function flattenSents(sents)
    #Define a new swahili query object
    output = swahiliQuery("-")
    words = []
    feats = []
    print()
    for sentI = 1:length(sents)
        sent = sents[sentI]
        if sentI % 100 == 0
            print("\r" * string(sentI))
        end
        for word in sent
            push!(words, word)
            append!(feats, keys(getFeatureValues(word)))
        end
    end
    output.words = words
    output.uniques = unique(feats)
    print()
    return output
end

function extractSents(filepath)
    sents = []
    if isdir(filepath)
        files = readdir(filepath)
        sents = @sync @parallel append! for file in files
            print("\r" * " "^50)
            print("\rProcessing " * filepath * file)
            getSentsFromJSON(filepath * file)
        end
    else
        sents = getSentsFromJSON(filepath)
    end
    println()
    return sents
end

#::Array{Array{swahiliTag, 1},1}
function convertSentsToHmm(sents)
    mutatedSents = @sync @parallel vcat for sent in sents
        toHmmTuple(sent)
    end
    return mutatedSents
end

function toHmmTuple(sent)
    hmmTuple = []#Array{Tuple{String, String}, 1}
    for word in sent
        push!(hmmTuple, (word.word, word.pos))
    end
    return [hmmTuple]
end
"""
This function changes the array of sparse vectors I was using before into a sparse matrix
"""
function sparseMatrix(data)
    rowIndexes::Array{Int64, 1} = []
    colIndexes::Array{Int64, 1} = []
    dataVals = []
    for rowI = 1:length(data)
        vecIndexes = data[rowI].nzind
        vecVals = data[rowI].nzval
        for (index, val) in zip(vecIndexes, vecVals)
            push!(rowIndexes, rowI)
            push!(colIndexes, index)
            push!(dataVals, val)
        end
    end
    res = SparseMatrixCSC{Int64, Int64}(sparse(rowIndexes, colIndexes, dataVals))
    return res
end
#function writeCONLL(query:: )
#=function createMatrix(query::swahiliQuery, cols::Array)
    println("There are " * string(length(cols)) * " features being extracted")
    println("There are " * string(length(query.words)) * " words to gather features for")
    words = query.words
    @inbounds classes = [word.cls for word in words]
    @inbounds pos = [word.pos for word in words]
    @inbounds msd = [word.msd for word in words]
    @inbounds obj = [word.obj for word in words]
    @inbounds subj = [word.subj for word in words]
    @inbounds tam = [word.tam for word in words]
    @inbounds rel = [word.rel for word in words]

    println("Current Progress")
    dataArray = []
    println()
    for wordI = 1:length(query.words)
        if wordI % 100 == 0
            print("\r" * string(wordI))
        end
        featsForExample = getFeatureValues(words[wordI])
        @inbounds tempRow = SparseArrays.sparsevec([feat in keys(featsForExample) ? Int64(featsForExample[feat]) : Int64(0) for feat in cols])
        #[tempRow]
        push!(dataArray, tempRow)
    end
    return Dict("classes"=>classes,
                "features"=>cols,
                "data"=>dataArray,
                "pos"=>pos,
                "msd"=>msd,
                "obj"=>obj,
                "subj"=>subj,
                "tam"=>tam,
                "rel"=>rel)
end=#
function createMatrix(query::swahiliQuery, cols::Array)
    println("There are " * string(length(cols)) * " features being extracted")
    println("There are " * string(length(query.words)) * " words to gather features for")
    words = query.words
    @inbounds classes = [word.cls for word in words]
    @inbounds pos = [word.pos for word in words]
    @inbounds msd = [word.msd for word in words]
    @inbounds obj = [word.obj for word in words]
    @inbounds subj = [word.subj for word in words]
    @inbounds tam = [word.tam for word in words]
    @inbounds rel = [word.rel for word in words]

    println("Current Progress")
    rowIndexes = []
    colIndexes = []
    vals = []
    println()
    for wordI = 1:length(query.words)
        if wordI % 100 == 0
            print("\r" * string(wordI))
        end
        featsForExample = getFeatureValues(words[wordI])
        tempRow = [(featI, Int64(featsForExample[cols[featI]])) for featI in eachindex(cols) if cols[featI] in keys(featsForExample)]
        tempColIndexes = [first(x) for x in tempRow]
        tempVals = [last(x) for x in tempRow]
        tempRowIndexes = fill(wordI, length(tempVals))
        append!(rowIndexes, tempRowIndexes)
        append!(colIndexes, tempColIndexes)
        append!(vals, tempVals)
    end
    dataArray = SparseMatrixCSC{Int64, Int64}(sparse(rowIndexes, colIndexes, vals))
    return Dict("classes"=>classes,
                "features"=>cols,
                "data"=>dataArray,
                "pos"=>pos,
                "msd"=>msd,
                "obj"=>obj,
                "subj"=>subj,
                "tam"=>tam,
                "rel"=>rel)
end
function sampleMatrix(n_samples, matrix)
    sampledTrain = Dict(key => sample(n_samples, matrix[key]) for key in keys(trainData))
    return sampledTrain
end

function sample(n_samples, iterables...)
    length_iter = length(iterables[1])
    #Check that iterables all match in length
    for iterable in iterables
        if length(iterable) != length_iter
            jl_error("Iterables differ in length")
            jl_error("Target Length: " * str(length_iter))
            jl_error("Current Iterable Length: " * str(length(iterable)))
        end
    end
    srand(32109)
    indexes = rand(1:length_iter, n_samples)
    new_iterables = []
    for iterable in iterables
        sampled_iter = [iterable[ind] for ind in indexes]
        push!(new_iterables, sampled_iter)
    end
    return new_iterables
end

#=function writeData(data::Dict, path::String, key::String="pos")
    println("Writing labels")
    outLabels = open(path * ".pkl", "w")
    pickle.dump(Dict(key => data[key] for key in keys(data)), outLabels) #if key != "data"), outLabels)
    close(outLabels)
    println("Writing data")
    outData = open(path, "w")
    for row in data["data"]
        fullRow = full(row)
        for cell in fullRow[1:end-1]
            write(outData, string(cell) * "\t")
        end
        write(outData, string(fullRow[end]) * "\n")
    end
    close(outData)
end

function writeData(data, path::String)
    println("Beginning to write out data")
    newData = py_sparse.csc_matrix(data["data"])
    data["data"] = newData
    println("Data has been changed into python sparse matrix")
    writePkl(data, path)
end
=#
function writeData(data, path::String)
    println("Writing data")
    mmwrite(path * ".mtx", data["data"])
    println("Writing pickle")
    writePkl(Dict(key=>data[key] for key in keys(data) if key != "data"), path * ".pkl")
end
function writePkl(data, path::String)
    out = open(path, "w")
    pickle.dump(data, out)
    close(out)
end

function split_datasets(ratio::T, iterables...) where {T<:AbstractFloat}
    length_iter = length(iterables[1])
    n_train_samples = Integer(ceil(ratio * length_iter))
    #Check that iterables all match in length
    for iterable in iterables
        if length(iterable) != length_iter
            jl_error("Iterables differ in length")
            jl_error("Target Length: " * str(length_iter))
            jl_error("Current Iterable Length: " * str(length(iterable)))
        end
    end
    srand(32109)
    indexes = rand(1:length_iter, n_train_samples)
    skipped_indexes = setdiff(1:length_iter, indexes)
    train_iterables = []
    test_iterables = []
    for iterable in iterables
        sampled_iter = [iterable[ind] for ind in indexes]
        push!(train_iterables, sampled_iter)
        sampled_iter = [iterable[ind] for ind in skipped_indexes]
        push!(test_iterables, sampled_iter)
    end
    return train_iterables, test_iterables
end

max_stars_repo_name: ksteimel/SwahiliMorphemeTagging

max_issues_repo_name: ksteimel/SwahiliMorphemeTagging

max_forks_repo_name: ksteimel/SwahiliMorphemeTagging


text:
"""

    heattransfer(myARGS)

Extract command-line arguments for computing the radiative heat transfer among a collection
of molecules, then compute the radiative heat transfer appropriately.

!!! note

    Any frequencies will be converted to their absolute values. Frequencies are not allowed
    to be exactly zero.

!!! note

    This computes the integrand
    ``\\operatorname{Tr}(\\operatorname{asym}(\\mathbb{V}_{m}^{-1}) \\mathbb{P}_{m} \\mathbb{T}^{\\dagger} \\operatorname{asym}(\\mathbb{P}_{n} \\mathbb{G}^{\\mathrm{env}}) \\mathbb{T}\\mathbb{P}_{m})``
    at specified real frequencies, without doing any integration. This includes neither
    numerical prefactors outside of the trace, nor the Planck function.

!!! tip

    This function can be embarrassingly parallelized over frequency and wavevector. However,
    the integrand at each frequency and wavevector is calculated for every transformation,
    as this maximizes the efficiency of the program by virtue of reusing quantities
    independent of geometric transformations.

"""
function heattransfer(myARGS)

    mollistfilename = readRequiredFilenameARG(myARGS, "mollistfilename", "AARMBEM_heat.jl");
    outfilename = readRequiredFilenameARG(myARGS, "outfilename", "AARMBEM_heat.jl");

    translistfilename = readOptionalFilenameARG(myARGS, "translistfilename");
    periodicfilename = readOptionalFilenameARG(myARGS, "periodicfilename");

    Genvstr = readGenvstr(myARGS);

    nophonons, outfilename = readOptionalBoolcondARG(myARGS, "nophonons", outfilename);
    nonretarded, outfilename = readOptionalBoolcondARG(myARGS, "nonretarded", outfilename);

    FT = readOptionalFloatType(myARGS);
    println("FT is ", FT);

    myAllMolData = readMolSystem(mollistfilename, nophonons, FT);
    myTransData = readTransData(myAllMolData, translistfilename);
    myPeriodicData = readPeriodicData(FT, periodicfilename);

    freqlist = Array{FT, 1}(undef, 0);
    readfreqlist!(freqlist, myARGS, "AARMBEM_heat.jl");

    klist = readklist(myARGS, "AARMBEM_heat.jl", myPeriodicData.numdims, FT);

    freqklist = hcat(vec(repeat(transpose(freqlist), size(klist, 1), 1)),
                     repeat(klist, length(freqlist), 1));

    heattransfer_general(outfilename, myAllMolData, myTransData, myPeriodicData, Genvstr,
                      nonretarded, freqklist)

end

"""

    heattransfer_general(outfilename, myAllMolData, myTransData, myPeriodicData, Genvstr,
                         nonretarded, freqklist)

Compute the actual radiative heat transfer integrand and print it to `outfilename`.

"""
function heattransfer_general(outfilename::AbstractString,
                              myAllMolData::MolSystem{<:OneMol{FT}},
                              myTransData::TransData{FT},
                              myPeriodicData::PeriodicData{FT},
                              Genvstr::AbstractString, nonretarded::Bool,
                              freqklist::Array{FT, 2}) where FT<:AbstractFloat

    GFSCAGG! = GFVACGG_sca!;
    if (Genvstr == "PEC")
        GFSCAGG! = GFPECGG!;
    end

    totalTinvmat = SharedArray{Complex{FT}}(3*myAllMolData.cumnumatomslist[end],
                                            3*myAllMolData.cumnumatomslist[end]);
    totalGenvmat = SharedArray{Complex{FT}}(3*myAllMolData.cumnumatomslist[end],
                                            3*myAllMolData.cumnumatomslist[end]);
    totalTmat = SharedArray{Complex{FT}}(3*myAllMolData.cumnumatomslist[end],
                                         3*myAllMolData.cumnumatomslist[end]);
    Phimn = zeros(FT, myAllMolData.nummols, myAllMolData.nummols);

    for ll=1:size(freqklist, 1)

        currfreq = abs(freqklist[ll, 1]);

        currfreqG = currfreq;
        if (nonretarded)
            currfreqG = zero(currfreqG);
        end
        currk = vec(freqklist[ll, 2:4]);
        fill!(totalTinvmat, zero(eltype(totalTinvmat)));
        fill!(totalGenvmat, zero(eltype(totalGenvmat)));
        fill!(totalTmat, zero(eltype(totalTmat)));

        for nn=1:myAllMolData.nummols

            currumolind = myAllMolData.uniquemolinds[nn];

            if (currumolind > -1)

                constructMolAlpha!(myAllMolData.uniqueMolDataArray[currumolind],
                                   myPeriodicData, currfreq, currk);

                assemble1MolGvacinf!(myAllMolData.uniqueMolDataArray[currumolind],
                                     myPeriodicData, currfreqG, currk);

                assemble1MolGvacinfdiag!(myAllMolData.uniqueMolDataArray[currumolind],
                                         currfreqG);

            end

        end

        for tt=1:myTransData.numTrans

            transformAtomPos!(myTransData, myAllMolData, tt);

            for nn2=1:myAllMolData.nummols

                startidx2 = 1;
                if (nn2 > 1)
                    startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
                end
                currmolind2 = myAllMolData.uniquemolinds[nn2];
                if (currmolind2 < 0)
                    currmolind2 = myAllMolData.dupemolinds[nn2];
                end

                if (myTransData.changedFromBefore[nn2, tt])
                    assembleTinvGenvDiagBlock!(totalTinvmat, totalGenvmat, myAllMolData,
                                               myTransData, nn2, currmolind2, startidx2,
                                               tt, myPeriodicData, GFSCAGG!, currfreqG,
                                               currk);
                end

            end

            for nn2=1:myAllMolData.nummols

                startidx2 = 1;
                if (nn2 > 1)
                    startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
                end

                for nn1=(nn2+1):myAllMolData.nummols

                    startidx1 = 3*myAllMolData.cumnumatomslist[nn1-1] + 1;

                    if (myTransData.changedFromBefore[nn2, tt] ||
                        myTransData.changedFromBefore[nn1, tt])

                        assembleTinvGenvOffDiagBlock!(totalTinvmat, totalGenvmat,
                                                      myAllMolData, myTransData, nn1, nn2,
                                                      startidx1, startidx2, myPeriodicData,
                                                      GFSCAGG!, currfreqG, currk);
                    end
                end

            end

            if (myPeriodicData.numdims > 0)

                for nn2=1:myAllMolData.nummols

                    startidx2 = 1;
                    if (nn2 > 1)
                        startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
                    end

                    for nn1=(nn2+1):myAllMolData.nummols

                        startidx1 = 3*myAllMolData.cumnumatomslist[nn1-1] + 1;

                        if (myTransData.changedFromBefore[nn2, tt] ||
                            myTransData.changedFromBefore[nn1, tt])

                            assembleTinvGenvOffDiagBlock!(totalTinvmat, totalGenvmat,
                                                          myAllMolData, myTransData, nn2, nn1,
                                                          startidx2, startidx1, myPeriodicData,
                                                          GFSCAGG!, currfreqG, currk);
                        end
                    end

                end
            end

            totalTmat .= inv(totalTinvmat);

            assemblePhimn!(Phimn, myAllMolData, myPeriodicData, totalTmat, totalGenvmat);

            outfile = open(outfilename, "a+");
            if (myPeriodicData.numdims == 0)

                @printf(outfile, "%s %.16e", myTransData.transLabels[tt], abs(currfreq));
                for nn1=1:myAllMolData.nummols
                    for nn2=1:myAllMolData.nummols
                        @printf(outfile, " %.16e", Phimn[nn1, nn2]);;
                    end
                end
                @printf(outfile, "\n");
            else
                @printf(outfile, "%s %.16e %.16e %.16e %.16e", myTransData.transLabels[tt],
                        abs(currfreq), currk[1], currk[2], currk[3]);
                for nn1=1:myAllMolData.nummols
                    for nn2=1:myAllMolData.nummols
                        @printf(outfile, " %.16e", Phimn[nn1, nn2]);
                    end
                end
                @printf(outfile, "\n%s %.16e %.16e %.16e %.16e",
                        myTransData.transLabels[tt], abs(currfreq), -1*currk[1],
                        -1*currk[2], -1*currk[3]);
                for nn1=1:myAllMolData.nummols
                    for nn2=1:myAllMolData.nummols
                        @printf(outfile, " %.16e", Phimn[nn2, nn1]);
                    end
                end
                @printf(outfile, "\n");
            end
            close(outfile);

            untransformAtomPos!(myTransData, myAllMolData, tt);

        end

    end

end

"""

    assemblePhimn!(Phimn::Array{FT, 2},
                   myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                   myPeriodicData::PeriodicData{FT},
                   totalTmat::SharedArray{Complex{FT}, 2},
                   totalGenvmat::SharedArray{Complex{FT}, 2})

Assemble the matrix `Phimn`, whose elements `Phimn[nn1, nn2]` represent the energy transfer
from molecule `nn1` to molecule `nn2`, computed from trace expressions involving the
matrices `totalTmat` and `totalGenvmat`, in-place. If there is no spatial periodicity,
`Phimn` should be symmetric. If there is spatial periodicity, the transpose of `Phimn` at
wavevector ``\\vec{k}`` is the equivalent of computing `Phimn` at ``-\\vec{k}``. This assumes
the presence of phonons.

"""
function assemblePhimn!(Phimn::Array{FT, 2},
                        myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                        myPeriodicData::PeriodicData{FT},
                        totalTmat::SharedArray{Complex{FT}, 2},
                        totalGenvmat::SharedArray{Complex{FT}, 2}) where FT<:AbstractFloat

    for nn2=1:myAllMolData.nummols

        startidx2 = 1;
        if (nn2 > 1)
            startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
        end
        Genvview = view(totalGenvmat,
                        (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                        1:3*myAllMolData.cumnumatomslist[end]);

        startmol1 = 1;
        if (myPeriodicData.numdims == 0)
            startmol1 = nn2;
        end

        for nn1=startmol1:myAllMolData.nummols

            startidx1 = 1;
            if (nn1 > 1)
                startidx1 = 3*myAllMolData.cumnumatomslist[nn1-1] + 1;
            end
            Tview = view(totalTmat,
                         1:3*myAllMolData.cumnumatomslist[end],
                         (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]));

            currmolind1 = 1;
            if (myAllMolData.dupemolinds[nn1] == -1)
                currmolind1 = myAllMolData.uniquemolinds[nn1];
            else
                currmolind1 = myAllMolData.dupemolinds[nn1];
            end
            Phimn[nn1, nn2] =
            imag(tr(ahermpart(myAllMolData.uniqueMolDataArray[currmolind1].alphainv) *
                    adjoint(Tview)[:, (startidx2-1) .+
                                   (1:3*myAllMolData.numatomslist[nn2])] *
                    Genvview * Tview));
            if (myPeriodicData.numdims == 0)
                Phimn[nn2, nn1] = Phimn[nn1, nn2];
            end
        end

    end

end

"""

    assemblePhimn!(Phimn::Array{FT, 2},
                   myAllMolData::MolSystem{OneMol_NoPhonons{FT}},
                   myPeriodicData::PeriodicData{FT},
                   totalTmat::SharedArray{Complex{FT}, 2},
                   totalGenvmat::SharedArray{Complex{FT}, 2})

Assemble the matrix `Phimn`, whose elements `Phimn[nn1, nn2]` represent the energy transfer
from molecule `nn1` to molecule `nn2`, computed from trace expressions involving the
matrices `totalTmat` and `totalGenvmat`, in-place. If there is no spatial periodicity,
`Phimn` should be symmetric. If there is spatial periodicity, the transpose of `Phimn` at
wavevector ``\\vec{k}`` is the equivalent of computing `Phimn` at ``-\\vec{k}``. This assumes
the absence of phonons.

"""
function assemblePhimn!(Phimn::Array{FT, 2},
                        myAllMolData::MolSystem{OneMol_NoPhonons{FT}},
                        myPeriodicData::PeriodicData{FT},
                        totalTinvmat::SharedArray{Complex{FT}, 2},
                        totalGenvmat::SharedArray{Complex{FT}, 2}) where FT<:AbstractFloat

    for nn2=1:myAllMolData.nummols

        startidx2 = 1;
        if (nn2 > 1)
            startidx2 = 3*myAllMolData.cumnumatomslist[nn2-1] + 1;
        end
        Genvview = view(totalGenvmat,
                        (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                        1:3*myAllMolData.cumnumatomslist[end]);

        startmol1 = 1;
        if (myPeriodicData.numdims == 0)
            startmol1 = nn2;
        end

        for nn1=startmol1:myAllMolData.nummols

            startidx1 = 1;
            if (nn1 > 1)
                startidx1 = 3*myAllMolData.cumnumatomslist[nn1-1] + 1;
            end
            Tview = view(totalTmat,
                         1:3*myAllMolData.cumnumatomslist[end],
                         (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]));

            currmolind1 = 1;
            if (myAllMolData.dupemolinds[nn1] == -1)
                currmolind1 = myAllMolData.uniquemolinds[nn1];
            else
                currmolind1 = myAllMolData.dupemolinds[nn1];
            end
            Phimn[nn1, nn2] =
            imag(tr(mulMatDiagleft(adjoint(Tview)[:, (startidx2-1) .+
                                                  (1:3*myAllMolData.numatomslist[nn2])] *
                                   Genvview * Tview,
                                   vecNto3N(imag.(myAllMolData.uniqueMolDataArray[currmolind1].alphainv)))));
            if (myPeriodicData.numdims == 0)
                Phimn[nn2, nn1] = Phimn[nn1, nn2];
            end
        end

    end

end

"""

    assembleTinvGenvDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                               totalGenvmat::SharedArray{Complex{FT}, 2},
                               myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                               myTransData::TransData{FT}, nn2::Integer,
                               currmolind2::Integer, startidx2::Integer,
                               tt::Integer, myPeriodicData::PeriodicData{FT},
                               GFSCAGG!::Function,
                               currfreqG::Union{FT, Complex{FT}},
                               currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the molecular diagonal block of the scattering Green's function interaction matrix
of molecule `nn2` (whose unique molecular index `currmolind2` is precomputed and passed as
an argument) at frequency `currfreqG` and wavevector `currk` using `GFSCAGG!`, then add the
vacuum contribution of molecule `nn2`, performing any required matrix rotations
corresponding to transformation `tt` for molecule `nn2`. Stamp that added to coincident
contributions into
`totalGenvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`,
or that (without coincident contributions) subtracted from the inverse molecular
susceptibility into
`totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`.
Finally, undo any rotation corresponding to transformation `tt` in the inverse T-operator of
molecule `nn2` for reuse with future transformations. This assumes the presence of phonons.

"""
function assembleTinvGenvDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                                    totalGenvmat::SharedArray{Complex{FT}, 2},
                                    myAllMolData::MolSystem{OneMol_WithPhonons{FT}},
                                    myTransData::TransData{FT}, nn2::Integer,
                                    currmolind2::Integer, startidx2::Integer,
                                    tt::Integer, myPeriodicData::PeriodicData{FT},
                                    GFSCAGG!::Function,
                                    currfreqG::Union{FT, Complex{FT}},
                                    currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    tempview = view(totalGenvmat, (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                    (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]));

    fill!(tempview, zero(eltype(tempview)));

    assembleGscaMolBlock!(tempview, myAllMolData, myTransData, nn2, nn2, 1, 1,
                          myPeriodicData, GFSCAGG!, currfreqG, currk);

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].alphainv,
                 myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                 myTransData.baselineRot[3*nn2-2:3*nn2, :]);

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                 myTransData.baselineRot[3*nn2-2:3*nn2, :]);

    tempview .+= myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf[:, :];

    totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                 (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])] .=
    myAllMolData.uniqueMolDataArray[currmolind2].alphainv .- tempview;

    addtodiag!(tempview,
               vecNto3N(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinfdiag));

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].alphainv,
                 transpose(myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                           myTransData.baselineRot[3*nn2-2:3*nn2, :]));

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 transpose(myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                           myTransData.baselineRot[3*nn2-2:3*nn2, :]));

end

"""

    assembleTinvGenvDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                               totalGenvmat::SharedArray{Complex{FT}, 2},
                               myAllMolData::MolSystem{OneMol_NoPhonons{FT}},
                               myTransData::TransData{FT}, nn2::Integer,
                               currmolind2::Integer,
                               startidx2::Integer,
                               tt::Integer, myPeriodicData::PeriodicData{FT},
                               GFSCAGG!::Function,
                               currfreqG::Union{FT, Complex{FT}},
                               currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the molecular diagonal block of the scattering Green's function interaction matrix
of molecule `nn2` (whose unique molecular index `currmolind2` is precomputed and passed as
an argument) at frequency `currfreqG` and wavevector `currk` using `GFSCAGG!`, then add the
vacuum contribution of molecule `nn2`, performing any required matrix rotations
corresponding to transformation `tt` for molecule `nn2`. Stamp that added to coincident
contributions into
`totalGenvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`,
or that (without coincident contributions) subtracted from the inverse molecular
susceptibility into
`totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`.
Finally, undo any rotation corresponding to transformation `tt` in the inverse T-operator of
molecule `nn2` for reuse with future transformations. This assumes the absence of phonons.

"""
function assembleTinvGenvDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                                    totalGenvmat::SharedArray{Complex{FT}, 2},
                                    myAllMolData::MolSystem{OneMol_NoPhonons{FT}},
                                    myTransData::TransData{FT}, nn2::Integer,
                                    currmolind2::Integer,
                                    startidx2::Integer,
                                    tt::Integer, myPeriodicData::PeriodicData{FT},
                                    GFSCAGG!::Function,
                                    currfreqG::Union{FT, Complex{FT}},
                                    currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    tempview = view(totalGenvmat, (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                    (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]));

    fill!(tempview, zero(eltype(tempview)));

    assembleGscaMolBlock!(tempview, myAllMolData, myTransData, nn2, nn2, 1, 1,
                          myPeriodicData, GFSCAGG!, currfreqG, currk);

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                 myTransData.baselineRot[3*nn2-2:3*nn2, :]);

    tempview .+= myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf[:, :];

    totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                 (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])] .=
    subtractfromdiag(tempview,
                     vecNto3N(myAllMolData.uniqueMolDataArray[currmolind2].alphainv));

    addtodiag!(tempview,
               vecNto3N(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinfdiag));

    rotate3Nmat!(myAllMolData.uniqueMolDataArray[currmolind2].GvacGGinf,
                 transpose(myTransData.rotArray[3*nn2-2:3*nn2, 3*tt-2:3*tt] *
                           myTransData.baselineRot[3*nn2-2:3*nn2, :]));

end

"""

    assembleTinvGenvOffDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                                  totalGenvTinvmat::SharedArray{Complex{FT}, 2},
                                  myAllMolData::MolSystem{<:OneMol{FT}},
                                  myTransData::TransData{FT},
                                  nn1::Integer, nn2::Integer,
                                  startidx1::Integer,
                                  startidx2::Integer,
                                  myPeriodicData::PeriodicData{FT},
                                  GFSCAGG!::Function,
                                  currfreqG::Union{FT, Complex{FT}},
                                  currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

Assemble the molecular off-diagonal block of the scattering Green's function interaction
matrix of molecules `nn2` as the source and `nn1` as the field (whose respective unique
molecular indices `currmolind1` and `currmolind2` are precomputed and passed as arguments)
at frequency `currfreqG` and wavevector `currk` using `GFSCAGG!`, multiply on
the left by the electronic polarizabilities, add the vacuum contribution, and stamp that
into
`totalGenvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`
and its negative into
`totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]), (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])]`.
If there is no spatial periodicity, do the same for the (unconjugated) transpose, swapping
`nn2` and `nn1` without performing further Green's function computations.

"""
function assembleTinvGenvOffDiagBlock!(totalTinvmat::SharedArray{Complex{FT}, 2},
                                       totalGenvmat::SharedArray{Complex{FT}, 2},
                                       myAllMolData::MolSystem{<:OneMol{FT}},
                                       myTransData::TransData{FT},
                                       nn1::Integer, nn2::Integer,
                                       startidx1::Integer,
                                       startidx2::Integer,
                                       myPeriodicData::PeriodicData{FT},
                                       GFSCAGG!::Function,
                                       currfreqG::Union{FT, Complex{FT}},
                                       currk::Array{FT, 1}=zeros(FT, 3)) where FT<:AbstractFloat

    tempview = view(totalGenvmat, (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]),
                    (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]));
    fill!(tempview, zero(eltype(tempview)));

    assembleGvacMolBlock!(tempview, myAllMolData, myTransData, nn1, nn2, 1, 1,
                          myPeriodicData, currfreqG, currk);
    assembleGscaMolBlock!(tempview, myAllMolData, myTransData, nn1, nn2, 1, 1,
                          myPeriodicData, GFSCAGG!, currfreqG, currk);

    totalTinvmat[(startidx1-1).+(1:3*myAllMolData.numatomslist[nn1]),
                 (startidx2-1).+(1:3*myAllMolData.numatomslist[nn2])] .=
    -1 .* tempview;

    if (myPeriodicData.numdims == 0)
        totalGenvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                     (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1])] .=
        transpose(deepcopy(tempview));
        totalTinvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                     (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1])] .= -1 .*
        totalGenvmat[(startidx2-1).+(1:3*myAllMolData.numatomslist[nn2]),
                     (startidx1-1).+(1:3*myAllMolData.numatomslist[nn1])];
    end

end

ext: jl


text:
using Cxx
using Libdl

#const path_to_lib = "/home/jarv1s/pixhawk_sensor_data"
#const path_to_header = joinpath(@__DIR__, "..", "deps", "usr", "include")

#addHeaderDir(path_to_lib, kind=C_System)
#Libdl.dlopen(path_to_lib * "/libPixhawkData.so", Libdl.RTLD_GLOBAL)

cxx"#include <cstdlib>"
cxx"#include <stdio.h>"
cxx"#include <unistd.h>"
cxx"#include <fcntl.h>"
cxx"#include <termios.h>"
cxx"#include <pthread.h>"
cxx"#include <signal.h>"

cxx"#include <mavlink/include/mavlink/v2.0/common/mavlink.h>"

cxxinclude("serial_port.h")

#=cxx"""
class Serial_Port
{
public:

    Serial_Port();
    Serial_Port(const char *uart_name_, int baudrate_);
    void initialize_defaults();
    ~Serial_Port();

    bool debug;
    const char *uart_name;
    int  baudrate;
    int  status;

    int read_message(mavlink_message_t &message);
    int write_message(const mavlink_message_t &message);

    void open_serial();
    void close_serial();

    void start();
    void stop();

    void handle_quit( int sig );

private:

    int  fd;
    mavlink_status_t lastStatus;
    pthread_mutex_t  lock;

    int  _open_port(const char* port);
    bool _setup_port(int baud, int data_bits, int stop_bits, bool parity, bool hardware_control);
    int  _read_port(uint8_t &cp);
    int _write_port(char *buf, unsigned len);
};
"""=#

Serial_Port = @cxxnew Serial_Port()

cxx"""
Serial_Port::
Serial_Port(const char *uart_name_ , int baudrate_)
{
	initialize_defaults();
	uart_name = uart_name_;
	baudrate  = baudrate_;
}

Serial_Port::
Serial_Port()
{
	initialize_defaults();
}

Serial_Port::
~Serial_Port()
{
	// destroy mutex
	pthread_mutex_destroy(&lock);
}
"""

Serial_Port1() = @cxx ~Serial_Port()

cxx"""
void
Serial_Port::
initialize_defaults()
{
	// Initialize attributes
	debug  = false;
	fd     = -1;
	status = SERIAL_PORT_CLOSED;

	uart_name = (char*)"/dev/ttyUSB0";
	baudrate  = 57600;

	// Start mutex
	int result = pthread_mutex_init(&lock, NULL);
	if ( result != 0 )
	{
		printf("\n mutex init failed\n");
		throw 1;
	}
}
"""

initialize_defaults() = @cxx initialize_defaults()

cxx"""
int
Serial_Port::
read_message(mavlink_message_t &message)
{
	uint8_t          cp;
	mavlink_status_t status;
	uint8_t          msgReceived = false;

	// this function locks the port during read
	int result = _read_port(cp);

	if (result > 0)
	{
		// the parsing
		msgReceived = mavlink_parse_char(MAVLINK_COMM_1, cp, &message, &status);

		// check for dropped packets
		if ( (lastStatus.packet_rx_drop_count != status.packet_rx_drop_count) && debug )
		{
			printf("ERROR: DROPPED %d PACKETS\n", status.packet_rx_drop_count);
			unsigned char v=cp;
			fprintf(stderr,"%02x ", v);
		}
		lastStatus = status;
	}

	// Couldn't read from port
	else
	{
		fprintf(stderr, "ERROR: Could not read from fd %d\n", fd);
	}

	if(msgReceived && debug)
	{
		// Report info
		printf("Received message from serial with ID #%d (sys:%d|comp:%d):\n", message.msgid, message.sysid, message.compid);

		fprintf(stderr,"Received serial data: ");
		unsigned int i;
		uint8_t buffer[MAVLINK_MAX_PACKET_LEN];

		// check message is write length
		unsigned int messageLength = mavlink_msg_to_send_buffer(buffer, &message);

		// message length error
		if (messageLength > MAVLINK_MAX_PACKET_LEN)
		{
			fprintf(stderr, "\nFATAL ERROR: MESSAGE LENGTH IS LARGER THAN BUFFER SIZE\n");
		}

		// print out the buffer
		else
		{
			for (i=0; i<messageLength; i++)
			{
				unsigned char v=buffer[i];
				fprintf(stderr,"%02x ", v);
			}
			fprintf(stderr,"\n");
		}
	}

	// Done!
	return msgReceived;
}
"""

read_message() = @cxx read_message()

cxx"""
int
Serial_Port::
write_message(const mavlink_message_t &message)
{
	char buf[300];

	// Translate message to buffer
	unsigned len = mavlink_msg_to_send_buffer((uint8_t*)buf, &message);

	// Write buffer to serial port, locks port while writing
	int bytesWritten = _write_port(buf,len);

	return bytesWritten;
}
"""

write_message() = @cxx write_message()

cxx"""
void
Serial_Port::
open_serial()
{

	// --------------------------------------------------------------------------
	//   OPEN PORT
	// --------------------------------------------------------------------------
	//printf("OPEN PORT\n");

	fd = _open_port(uart_name);

	// Check success
	if (fd == -1)
	{
		printf("failure, could not open port.\n");
		throw EXIT_FAILURE;
	}

	// --------------------------------------------------------------------------
	//   SETUP PORT
	// --------------------------------------------------------------------------
	bool success = _setup_port(baudrate, 8, 1, false, false);

	// --------------------------------------------------------------------------
	//   CHECK STATUS
	// --------------------------------------------------------------------------
	if (!success)
	{
		printf("failure, could not configure port.\n");
		throw EXIT_FAILURE;
	}
	if (fd <= 0)
	{
		printf("Connection attempt to port %s with %d baud, 8N1 failed, exiting.\n", uart_name, baudrate);
		throw EXIT_FAILURE;
	}

	// --------------------------------------------------------------------------
	//   CONNECTED!
	// --------------------------------------------------------------------------
	//printf("Connected to %s with %d baud, 8 data bits, no parity, 1 stop bit (8N1)\n", uart_name, baudrate);
	lastStatus.packet_rx_drop_count = 0;

	status = true;

	printf("\n");

	return;

}
"""

open_serial() = @cxx open_serial()

cxx"""
void
Serial_Port::
close_serial()
{
	//printf("CLOSE PORT\n");

	int result = close(fd);

	if ( result )
	{
		fprintf(stderr,"WARNING: Error on port close (%i)\n", result );
	}

	status = false;

	printf("\n");

}
"""

close_serial() = @cxx close_serial()

cxx"""
void
Serial_Port::
start()
{
	open_serial();
}
"""

start() = @cxx start()

cxx"""
void
Serial_Port::
stop()
{
	close_serial();
}
"""

stop() = @cxx stop()

cxx"""
void
Serial_Port::
handle_quit( int sig )
{
	try {
		stop();
	}
	catch (int error) {
		fprintf(stderr,"Warning, could not stop serial port\n");
	}
}
"""

handle_quit() = @cxx handle_quit()

cxx"""
int
Serial_Port::
_open_port(const char* port)
{
	// Open serial port
	// O_RDWR - Read and write
	// O_NOCTTY - Ignore special chars like CTRL-C
	fd = open(port, O_RDWR | O_NOCTTY | O_NDELAY);

	// Check for Errors
	if (fd == -1)
	{
		/* Could not open the port. */
		return(-1);
	}

	// Finalize
	else
	{
		fcntl(fd, F_SETFL, 0);
	}

	// Done!
	return fd;
}
"""

_open_port() = @cxx _open_port()

cxx"""
bool
Serial_Port::
_setup_port(int baud, int data_bits, int stop_bits, bool parity, bool hardware_control)
{
	// Check file descriptor
	if(!isatty(fd))
	{
		fprintf(stderr, "\nERROR: file descriptor %d is NOT a serial port\n", fd);
		return false;
	}

	// Read file descritor configuration
	struct termios  config;
	if(tcgetattr(fd, &config) < 0)
	{
		fprintf(stderr, "\nERROR: could not read configuration of fd %d\n", fd);
		return false;
	}

	config.c_iflag &= ~(IGNBRK | BRKINT | ICRNL |
						INLCR | PARMRK | INPCK | ISTRIP | IXON);

	config.c_oflag &= ~(OCRNL | ONLCR | ONLRET |
						 ONOCR | OFILL | OPOST);

	#ifdef OLCUC
		config.c_oflag &= ~OLCUC;
	#endif

	#ifdef ONOEOT
		config.c_oflag &= ~ONOEOT;
	#endif

	config.c_lflag &= ~(ECHO | ECHONL | ICANON | IEXTEN | ISIG);

	config.c_cflag &= ~(CSIZE | PARENB);
	config.c_cflag |= CS8;

	// One input byte is enough to return from read()
	// Inter-character timer off
	config.c_cc[VMIN]  = 1;
	config.c_cc[VTIME] = 10; // was 0

	// Apply baudrate
	switch (baud)
	{
		case 1200:
			if (cfsetispeed(&config, B1200) < 0 || cfsetospeed(&config, B1200) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;
		case 1800:
			cfsetispeed(&config, B1800);
			cfsetospeed(&config, B1800);
			break;
		case 9600:
			cfsetispeed(&config, B9600);
			cfsetospeed(&config, B9600);
			break;
		case 19200:
			cfsetispeed(&config, B19200);
			cfsetospeed(&config, B19200);
			break;
		case 38400:
			if (cfsetispeed(&config, B38400) < 0 || cfsetospeed(&config, B38400) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;
		case 57600:
			if (cfsetispeed(&config, B57600) < 0 || cfsetospeed(&config, B57600) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;
		case 115200:
			if (cfsetispeed(&config, B115200) < 0 || cfsetospeed(&config, B115200) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;

		case 460800:
			if (cfsetispeed(&config, B460800) < 0 || cfsetospeed(&config, B460800) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;
		case 921600:
			if (cfsetispeed(&config, B921600) < 0 || cfsetospeed(&config, B921600) < 0)
			{
				fprintf(stderr, "\nERROR: Could not set desired baud rate of %d Baud\n", baud);
				return false;
			}
			break;
		default:
			fprintf(stderr, "ERROR: Desired baud rate %d could not be set, aborting.\n", baud);
			return false;

			break;
	}

	// Finally, apply the configuration
	if(tcsetattr(fd, TCSAFLUSH, &config) < 0)
	{
		fprintf(stderr, "\nERROR: could not set configuration of fd %d\n", fd);
		return false;
	}

	// Done!
	return true;
}
"""

_setup_port() = @cxx _setup_port()

cxx"""
int
Serial_Port::
_read_port(uint8_t &cp)
{

	// Lock
	pthread_mutex_lock(&lock);

	int result = read(fd, &cp, 1);

	// Unlock
	pthread_mutex_unlock(&lock);

	return result;
}
"""

_read_port() = @cxx _read_port()

cxx"""
int
Serial_Port::
_write_port(char *buf, unsigned len)
{

	// Lock
	pthread_mutex_lock(&lock);

	// Write packet via serial link
	const int bytesWritten = static_cast<int>(write(fd, buf, len));

	// Wait until all data has been written
	tcdrain(fd);

	// Unlock
	pthread_mutex_unlock(&lock);

	return bytesWritten;
}
"""

_write_port() = @cxx _write_port()


text:
using HorizonSideRobots

function putpointcorners(r)
    start = []
    start = gocorner(r)
    pointcorners(r)
    goback(r,start)
end

function counter(r,side)
    c = 0
    while isborder(r,side) == 0
        move!(r,side)
        c += 1
    end
    return c
end

function gocorner(r)
start = []
    while isborder(r,Nord)==0 || isborder(r,West)==0
        x=counter(r,Nord)
        y=counter(r,West)
        push!(start,x)
        push!(start,y)
    end
    return start
end

function pointcorners(r)
    for side in [Sud,Ost,Nord,West]
        while isborder(r,side)==0
            move!(r,side)
        end
    putmarker!(r)
    end
end

function goback(r,start)
    reverse!(start)
    for i = 1:length(start)
        if i%2==0
            cmove!(r,Sud,start[i])
        else
            cmove!(r,Ost,start[i])
        end
    end
end

function cmove!(r,side,c)
    for i in 1:c
        move!(r,side)
    end
end

putpointcorners(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
using KitBase, Plots
using KitBase.ProgressMeter: @showprogress

set = Setup(case = "sod", space = "1d1f1v", maxTime = 0.12)
ps = PSpace1D(0.0, 1.0, 100, 1)
vs = VSpace1D(-5.0, 5.0, 72)
gas = Gas(Kn = 1e-4, K = 0.0, γ = 3.0)
fw, ff, bc, p = config_ib(set, ps, vs, gas)
ib = IB1F(fw, ff, bc, p)

ks = SolverSet(set, ps, vs, gas, ib)
ctr, face = init_fvm(ks, :static_array; structarray = true)

t = 0.0
dt = KitBase.timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int
@showprogress for iter = 1:nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, zeros(3))
end

plot(ks, ctr)


text:
"""
Function to obtain the output current time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
the dynamic device and bus voltage. It is dispatched for device type to compute the specific current.

"""
function compute_output_current(
    res::SimulationResults,
    dynamic_device::G,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    dt::Union{Nothing, Float64},
) where {G <: PSY.DynamicGenerator}

    #Obtain Data
    sys = get_system(res)

    #Get machine
    machine = PSY.get_machine(dynamic_device)
    Sbase = PSY.get_base_power(sys)
    basepower = PSY.get_base_power(dynamic_device)
    base_power_ratio = basepower / Sbase
    return _machine_current(
        machine,
        PSY.get_name(dynamic_device),
        V_R,
        V_I,
        base_power_ratio,
        res,
        dt,
    )
end

"""
Function to obtain the field current time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
the dynamic device and bus voltage. It is dispatched for device type to compute the specific current.

"""
function compute_field_current(
    res::SimulationResults,
    dynamic_device::G,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    dt::Union{Nothing, Float64},
) where {G <: PSY.DynamicGenerator}

    #Obtain Data
    sys = get_system(res)

    #Get machine
    machine = PSY.get_machine(dynamic_device)
    return _field_current(machine, PSY.get_name(dynamic_device), V_R, V_I, res, dt)
end

"""
Function to obtain the field voltage time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
the dynamic device and bus voltage. It is dispatched for device type to compute the specific voltage.

"""
function compute_field_voltage(
    res::SimulationResults,
    dynamic_device::G,
    dt::Union{Nothing, Float64},
) where {G <: PSY.DynamicGenerator}

    #Get AVR
    avr = PSY.get_avr(dynamic_device)
    return _field_voltage(avr, PSY.get_name(dynamic_device), res, dt)
end

"""
Function to obtain the output current time series of a Classic Machine model out of the DAE Solution. It is dispatched via the machine type.

"""
function _machine_current(
    machine::PSY.BaseMachine,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)

    R = PSY.get_R(machine)
    Xd_p = PSY.get_Xd_p(machine)
    eq_p = PSY.get_eq_p(machine)

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ)
    I_I = similar(δ)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]
        #Obtain electric current
        i_dq[1] = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * (eq_p .- V_q) - R * V_d)  #15.36
        i_dq[2] = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * V_d + R * (eq_p .- V_q)) #15.36

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the output current time series of a One-D-One-Q model out of the DAE Solution. It is dispatched via the machine type.
"""
function _machine_current(
    machine::PSY.OneDOneQMachine,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ed_p = post_proc_state_series(res, (name, :ed_p), dt)

    R = PSY.get_R(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ, Float64)
    I_I = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]
        #Obtain electric current
        i_dq[1] =
            (1.0 / (R^2 + Xd_p * Xq_p)) * (Xq_p * (eq_p[ix] - V_q) + R * (ed_p[ix] - V_d))  #15.32
        i_dq[2] =
            (1.0 / (R^2 + Xd_p * Xq_p)) * (-Xd_p * (ed_p[ix] - V_d) + R * (eq_p[ix] - V_q))  #15.32

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the output current time series of a SimpleMarconato or SimpleAndersonFouad model out of the DAE Solution. It is dispatched via the machine type.
"""
function _machine_current(
    machine::PSY.Union{PSY.SimpleMarconatoMachine, PSY.SimpleAFMachine},
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_pp = post_proc_state_series(res, (name, :eq_pp), dt)
    _, ed_pp = post_proc_state_series(res, (name, :ed_pp), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ, Float64)
    I_I = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]
        #Obtain electric current
        i_dq[1] =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (Xq_pp * (eq_pp[ix] - V_q) + R * (ed_pp[ix] - V_d))      #15.25
        i_dq[2] =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (-Xd_pp * (ed_pp[ix] - V_d) + R * (eq_pp[ix] - V_q))      #15.25

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the output current time series of a Marconato or AndersonFouad model out of the DAE Solution. It is dispatched via the machine type.
"""
function _machine_current(
    machine::Union{PSY.MarconatoMachine, PSY.AndersonFouadMachine},
    name::String,
    ::Vector{Float64},
    ::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_pp = post_proc_state_series(res, (name, :eq_pp), dt)
    _, ed_pp = post_proc_state_series(res, (name, :ed_pp), dt)
    _, ψd = post_proc_state_series(res, (name, :ψd), dt)
    _, ψq = post_proc_state_series(res, (name, :ψq), dt)

    #Get parameters
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ, Float64)
    I_I = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]

        #Obtain electric current
        i_dq[1] = (1.0 / Xd_pp) * (eq_pp[ix] - ψd[ix])      #15.18
        i_dq[2] = (1.0 / Xq_pp) * (-ed_pp[ix] - ψq[ix])     #15.18

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the output current time series of a GENROU/GENROE model out of the DAE Solution. It is dispatched via the machine type.
"""
function _machine_current(
    machine::Union{PSY.RoundRotorQuadratic, PSY.RoundRotorExponential},
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ed_p = post_proc_state_series(res, (name, :ed_p), dt)
    _, ψ_kd = post_proc_state_series(res, (name, :ψ_kd), dt)
    _, ψ_kq = post_proc_state_series(res, (name, :ψ_kq), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    ψq_pp = γ_q1 * ed_p + ψ_kq * (1 - γ_q1)
    ψd_pp = γ_d1 * eq_p + γ_d2 * (Xd_p - Xl) * ψ_kd

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ, Float64)
    I_I = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]

        #Obtain electric current
        i_dq[1] =
            (1.0 / (R^2 + Xq_pp * Xd_pp)) *
            (-R * (V_d - ψq_pp[ix]) + Xq_pp * (-V_q + ψd_pp[ix]))
        i_dq[2] =
            (1.0 / (R^2 + Xq_pp * Xd_pp)) *
            (Xd_pp * (V_d - ψq_pp[ix]) + R * (-V_q + ψd_pp[ix]))

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the output current time series of a GENSAL/GENSAE model out of the DAE Solution. It is dispatched via the machine type.
"""
function _machine_current(
    machine::Union{PSY.SalientPoleQuadratic, PSY.SalientPoleExponential},
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    base_power_ratio::Float64,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ψ_kd = post_proc_state_series(res, (name, :ψ_kd), dt)
    _, ψq_pp = post_proc_state_series(res, (name, :ψq_pp), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)

    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd

    i_dq = Vector{Float64}(undef, 2)
    I_R = similar(δ, Float64)
    I_I = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]

        #Obtain electric current
        i_dq[1] =
            (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp[ix]) + Xd_pp * (ψd_pp[ix] - V_q))
        i_dq[2] =
            (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d + ψq_pp[ix]) + R * (ψd_pp[ix] - V_q))

        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
    end
    return ts, I_R, I_I
end

"""
Function to obtain the field current time series of a Dynamic Generator. It is dispatched via the machine type.
By default, machine does not have support for field current

"""
function _field_current(
    machine::M,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    res::SimulationResults,
    dt::Union{Nothing, Float64},
) where {M <: PSY.Machine}
    @warn("Field current is not supported in the machine type $(M). Returning zeros.")
    ts, _ = post_proc_state_series(res, (name, :δ), dt)
    I_fd = zeros(length(ts))
    return ts, I_fd
end

"""
Function to obtain the field current time series of a Dynamic Generator with machine type GENROU/GENROE.

"""
function _field_current(
    machine::M,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    res::SimulationResults,
    dt::Union{Nothing, Float64},
) where {M <: Union{PSY.RoundRotorQuadratic, PSY.RoundRotorExponential}}
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ed_p = post_proc_state_series(res, (name, :ed_p), dt)
    _, ψ_kd = post_proc_state_series(res, (name, :ψ_kd), dt)
    _, ψ_kq = post_proc_state_series(res, (name, :ψ_kd), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #Additional Fluxes
    ψq_pp = γ_q1 * ed_p + ψ_kq * (1 - γ_q1)
    ψd_pp = γ_d1 * eq_p + γ_d2 * (Xd_p - Xl) * ψ_kd
    ψ_pp = sqrt.(ψd_pp .^ 2 + ψq_pp .^ 2)

    I_fd = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]

        #Obtain electric current
        I_d =
            (1.0 / (R^2 + Xq_pp * Xd_pp)) *
            (-R * (V_d - ψq_pp[ix]) + Xq_pp * (-V_q + ψd_pp[ix]))

        Se = saturation_function(machine, ψ_pp[ix])
        I_fd[ix] =
            eq_p[ix] +
            (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd[ix] + γ_d2 * eq_p[ix]) +
            Se * ψd_pp[ix]
    end
    return ts, I_fd
end

"""
Function to obtain the field current time series of a Dynamic Generator with machine type GENSAL.

"""
function _field_current(
    machine::PSY.SalientPoleQuadratic,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ψ_kd = post_proc_state_series(res, (name, :ψ_kd), dt)
    _, ψq_pp = post_proc_state_series(res, (name, :ψq_pp), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd

    I_fd = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]

        #Obtain electric current
        I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp[ix]) + Xd_pp * (ψd_pp[ix] - V_q))
        Se = saturation_function(machine, eq_p[ix])
        I_fd[ix] =
            eq_p[ix] +
            Se * eq_p[ix] +
            (Xd - Xd_p) * (I_d + γ_d2 * (eq_p[ix] - ψ_kd[ix] - (Xd_p - Xl) * I_d))
    end
    return ts, I_fd
end

"""
Function to obtain the field current time series of a Dynamic Generator with machine type GENSAE.

"""
function _field_current(
    machine::PSY.SalientPoleExponential,
    name::String,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, δ = post_proc_state_series(res, (name, :δ), dt)
    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
    _, ψ_kd = post_proc_state_series(res, (name, :ψ_kd), dt)
    _, ψq_pp = post_proc_state_series(res, (name, :ψq_pp), dt)

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd
    ψ_pp = sqrt.(ψd_pp .^ 2 + ψq_pp .^ 2)

    I_fd = similar(δ, Float64)

    for ix in 1:length(δ)
        v = δ[ix]
        V_d, V_q = ri_dq(v) * [V_R[ix]; V_I[ix]]

        #Obtain electric current
        I_d =
            (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp[ix]) + Xq_pp * (-V_q + ψd_pp[ix]))
        Se = saturation_function(machine, ψ_pp[ix])
        I_fd[ix] =
            eq_p[ix] +
            Se * ψd_pp[ix] +
            (Xd - Xd_p) * (I_d + γ_d2 * (eq_p[ix] - ψ_kd[ix] - (Xd_p - Xl) * I_d))
    end
    return ts, I_fd
end

"""
Function to obtain the field voltage time series of a Dynamic Generator with avrs that have
the field voltage as a state. By default it is assumed that the models have that state.

"""
function _field_voltage(
    ::A,
    name::String,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
) where {A <: PSY.AVR}
    return post_proc_state_series(res, (name, :Vf), dt)
end

"""
Function to obtain the field voltage time series of a Dynamic Generator with avr AVRFixed.

"""
function _field_voltage(
    avr::PSY.AVRFixed,
    name::String,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    Vf0 = PSY.get_Vf(avr)
    ts, _ = post_proc_state_series(res, (name, :δ), dt)
    Vf = Vf0 * ones(length(ts))
    return ts, Vf
end

"""
Function to obtain the field voltage time series of a Dynamic Generator with avr ESAC1A and EXAC1.

"""
function _field_voltage(
    avr::Union{PSY.ESAC1A, PSY.EXAC1},
    name::String,
    res::SimulationResults,
    dt::Union{Nothing, Float64},
)
    ts, Ve = post_proc_state_series(res, (name, :Ve), dt)
    _, Xad_Ifd = post_proc_field_current_series(res, name, dt)
    Kc = PSY.get_Kc(avr)
    I_N = Kc * Xad_Ifd ./ Ve
    Vf = Ve .* rectifier_function.(I_N)
    return ts, Vf
end


text:
using BioSequences, DataFrames, CairoMakie

import ..enzyme_list
using ..wgregseq: design.import_primer

dir = @__DIR__
home_dir = joinpath(split(dir, "/")[1:end-1])

function check_dataframe(df; print_results=true, site_start=1, site_end=nothing)
    gdf = groupby(df, :promoter)
    max_list, min_list = Float64[], Float64[]
    output=""
    for _df in gdf
        promoter = _df.promoter |> unique

        if print_results
            println("Promoter: $promoter")
            println("-------------------------")
        else
            output *= "Promoter: $promoter\n"
            output *= "-------------------------\n"
        end

        mut_rate = mutation_coverage(_df, site_start, site_end)
        min_rate, max_rate = minimum(mut_rate), maximum(mut_rate)
        if print_results
            println("Minimum mutation rate: $(min_rate)")
            println("Maximum mutation rate: $(max_rate)")
            println()

        else
            output *= "Minimum mutation rate: $(min_rate)\n"
            output *= "Maximum mutation rate: $(max_rate)\n"
            output *= "\n"

        end
        push!(min_list, min_rate)
        push!(max_list, max_rate)
        if check_cut_sites(_df)
            if print_results
                println("Cut sites are correct.")
                println()
            else
                output *= "Cut sites are correct.\n"
                output *= "\n"
            end
        else
            if print_results
                println("Cut sites NOT are correct.")
                println()
            else
                output *= "Cut sites NOT are correct.\n"
                output *= "\n"
            end

        end

        if check_primers(_df)
            if print_results
                println("Primers are correct.")
                println()
            else
                output *= "Primers are correct.\n"
                output *= "\n"
            end
        else
            if print_results
                println("Primes NOT are correct.")
                println()
            else
                output *= "Primes NOT are correct.\n"
                output *= "\n"
            end

        end

        if print_results
            println()
            println()
        else
            output *= "\n"
            output *= "\n"
        end
    end

    if ~print_results
        open("/$home_dir/qc_output.txt", "w") do file
            write(file, output)
        end
    end
    check_primers_off_target(df)
    promoter = df.promoter |> unique
    fig = Figure(resolution=(15*length(promoter), 800))
    ax = Axis(fig[1, 1])

    scatter!(ax, 1:length(promoter), max_list, label="Maximum")
    scatter!(ax, 1:length(promoter), min_list, label="Minimum")
    lines!(ax, [1, length(promoter)], [0.1, 0.1], color="grey", linestyle=:dash)
    ax.ylabel = "Mutation Rate"
    ax.xticklabelrotation = pi/4
    ax.xticks = 1:length(promoter)
    ax.xtickformat = x -> string.(promoter)
    axislegend()
    save("/$home_dir/figures/min_max_mutation_rates.pdf", fig)
    return fig
end

function mutation_coverage(df, site_start=1, site_end=nothing)

    if isnothing(site_end)
        site_end = length(df.sequence[1])
    end

    seq_list = [seq[site_start:site_end] for seq in df.sequence]
    mat = PFM(seq_list)
    mat = mat ./ length(seq_list)
    mut_rate = 1 .- [maximum(A) for A in eachcol(mat)]
    return mut_rate
end

function check_cut_sites(df)
    gdf = groupby(df, ["upstream_re_site", "downstream_re_site"])
    sites_correct = true
    for _df in gdf
        enz1 = unique(_df.upstream_re_site)[1]
        enz2 = unique(_df.downstream_re_site)[1]

        site1 = enzyme_list[enzyme_list.enzyme .== enz1, "site"][1]
        site2 = enzyme_list[enzyme_list.enzyme .== enz2, "site"][1]
        sites_correct *= ~any(map(seq -> seq[21:26] != LongDNA{4}(site1),  _df.sequence))
        sites_correct *= ~any(map(seq -> seq[187:192] != LongDNA{4}(site2),  _df.sequence))
    end
    return sites_correct
end

function check_primers(df)
    primer_cols = filter(x -> occursin("primer", x), names(df))
    gdf = groupby(df, primer_cols)
    primers_correct = true
    for _df in gdf
        for col in primer_cols
            direction = split(col, "_")[1] |> string
            ind, positions  = unique(_df[!, col])[1]

            primer_seq = import_primer(ind, direction)[1:(positions[2]-positions[1]+1)]
            primers_correct *= ~any(map(seq -> seq[positions[1]:positions[2]] != LongDNA{4}(primer_seq),  _df.sequence))
        end
    end
    return primers_correct
end

function check_primers_off_target(df)
    primer_cols = filter(x -> occursin("rev", x), names(df))
    for col in primer_cols
        primers = unique(df[!, col])
        for primer in primers
            primer_len = primer[2][2] - primer[2][1] + 1
            println(primer)
            primer_seq = import_primer(primer[1], "rev")[1:primer_len]
            query = ApproximateSearchQuery(primer_seq)
            inds0 = [~isnothing(findfirst(query, 0, seq)) for seq in df.sequence]
            x0 = (inds0 .!= map(x -> x == primer, df[!, col]))

            x1 = [findfirst(query, 1, seq) for seq in df.sequence]
            x1 = x1[.~ isnothing.(x1)]
            x1 = x1[length.(collect.(x1)) .== primer_len]

            x2 = [findfirst(query, 2, seq) for seq in df.sequence]
            x2 = x2[.~ isnothing.(x2)]
            x2 = x2[length.(collect.(x2)) .== primer_len]

            if sum(x0) > 0
                println("$(sum(x0)) off target binding for primer $(primer[1]) with no mismatches.")
            elseif length(x1) > 0
                println("$(length(x1)) off target binding for primer $(primer[1]) with one mismatch.")
            elseif length(x2) > 0
                println("$(length(x2)) off target binding for primer $(primer[1]) with two mismatches.")
            else
                println("No off target binding for primer $(primer[1]).")
            end
            println()
        end
    end
end


text:
#=
test_coverage:
- Julia version:
- Author: Berk
- Date: 2022-01-26
=#

# Run julia17 --project=test/. --code-coverage test/all.jl
# Then restart julia from within OptimalConstraintTree and run this file.
using Pkg
Pkg.activate("test")
global PROJECT_ROOT = @__DIR__
ENV["CODECOV_TOKEN"]= "dfef6217-c401-4f55-91b7-efbf1afc6ca0"
using Coverage
# process '*.cov' files
coverage = process_folder() # defaults to src/; alternatively, supply the folder name as argument
# process '*.info' files
coverage = merge_coverage_counts(coverage, filter!(
    let prefixes = (joinpath(pwd(), "src", ""))
        c -> any(p -> startswith(c.filename, p), prefixes)
    end,
    LCOV.readfolder("test")))
# Get total coverage for all Julia files
covered_lines, total_lines = get_summary(coverage)
# Or process a single file
@show get_summary(process_file(joinpath("src", "OptimalTrees.jl")))
println("Covered lines: ", covered_lines)
println("Total lines: ", total_lines)
println("Ratio: ", covered_lines/total_lines)

# Submit to coverage
Codecov.submit_local(coverage)


text:
"""
    upperbound(x)

The largest value that _x_ may actually hold.

See also: [`lowerbound`](@Ref), [`interval`](@ref), [`setinterval`](@ref), [`ball`](@ref), [`setball`](@ref)
""" upperbound

"""
    lowerbound(x)

The smallest value that _x_ may actually hold.

See also: [`upperbound`](@Ref), [`interval`](@ref), [`setinterval`](@ref), [`ball`](@ref), [`setball`](@ref)
""" lowerbound

"""
    interval(x)

The smallest value and the largest value that _x_ may actually hold.

See also: [`lowerbound`](@Ref), [`uppebround`](@ref), [`setinterval`](@ref), [`ball`](@ref), [`setball`](@ref)
""" interval

"""
    setinterval(lowerbound, upperbound)

Returns the bounded interval as an `ArbReal` or an `ArbComplex`.

See also: [`lowerbound`](@Ref), [`uppebround`](@ref), [`interval`](@ref), [`ball`](@ref), [`setball`](@ref)
""" setinterval

"""
    midpoint(x)

The midpoint of the ball represented by _x_.

See also: [`radius`](@ref), [`ball`](@ref), [`setball`](@ref), [`interval`](@ref), [`setinterval`](@ref)
""" midpoint

"""
    radius(x)

The radius (half-width) of the ball represented by _x_.

See also: [`midpoint`](@ref), [`ball`](@ref) [`setball`](@ref), [`interval`](@ref), [`setinterval`](@ref)
""" radius

"""
    ball(x)

The midpoint and radius of the ball represented by _x_.

See also:  [`midpoint`](@ref), [`radius`](@ref), [`setball`](@ref), [`interval`](@ref), [`setinterval`](@ref)
""" ball

"""
    setball(midpoint, radius)

Returns the ball given as _midpoint_ and _radius_.

See also:  [`midpoint`](@ref), [`radius`](@ref), [`ball`](@ref), [`interval`](@ref), [`setinterval`](@ref)
""" setball

function setinterval(lo::ArbFloat{P}, hi::ArbFloat{P}) where {P}
    lo > hi && return setinterval(hi, lo)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_interval_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}, Ref{ArbFloat}, Clong), z, lo, hi, P)
    return z
end

function setinterval(lo::ArbReal{P}, hi::ArbReal{P}) where {P}
    lo > hi && return setinterval(hi, lo)
    setinterval(lowerbound(lo, ArbFloat), upperbound(hi, ArbFloat))
end

function interval(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    ArbFloat{P}(lowerbound(x)), ArbFloat{P}(upperbound(x))
end

function interval(x::ArbReal{P}) where {P}
    lowerbound(x), upperbound(x)
end

function setball(mid::ArbFloat{P}, rad::ArbFloat{P}) where {P}
    signbit(rad) && throw(ErrorException("nonnegative radius required ($rad)"))
    lbound = mid - rad
    ubound = mid + rad
    setinterval(lbound, ubound)
end

function setball(mid::ArbReal{P}, rad::ArbReal{P}) where {P}
    signbit(rad) && throw(ErrorException("nonnegative radius required ($rad)"))
    setball(ArbFloat{P}(mid), ArbFloat{P}(rad))
end

@inline ArbReal{P}(mid::ArbFloat{P}, rad::ArbFloat{P}) where {P} = setball(mid, rad)
@inline ArbReal{P}(mid::ArbReal{P}, rad::ArbReal{P}) where {P} = setball(midpoint(mid), midpoint(rad))

function ball(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    ArbFloat{P}(midpoint(x)), ArbFloat{P}(radius(x))
end

function ball(x::ArbReal{P}) where {P}
    midpoint(x), radius(x)
end

function ball(x::ArbFloat{P}) where {P}
    x, zero(typeof(x))
end

@inline function upperbound(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_ubound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function upperbound(x::ArbReal{P}) where {P}
    w = upperbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function lowerbound(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_lbound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function lowerbound(x::ArbReal{P}) where {P}
    w = lowerbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function upperbound_abs(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_abs_ubound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function upperbound_abs(x::ArbReal{P}) where {P}
    w = upperbound_abs(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function lowerbound_abs(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_abs_lbound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function lowerbound_abs(x::ArbReal{P}) where {P}
    w = lowerbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

function interval_abs(x::ArbReal{P}) where {P}
    lo = lowerbound_abs(x)
    hi = upperbound_abs(x)
    return lo, hi
end

function increase_radius(x::ArbReal{P}, err::ArbFloat{P}) where {P}
    err >= 0 || throw(ErrorException("nonnegative err required ($err)"))
    ccall(@libarb(arb_add_error_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, err)
    return x
end

function increase_radius(x::ArbReal{P}, err::ArbReal{P}) where {P}
    err >= 0 || throw(ErrorException("nonnegative err required ($err)"))
    ccall(@libarb(arb_add_error_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, err)
    return x
end

increase_radius(x::ArbReal{P}) where {P} = increase_radius(x, ulp(x))

function decrease_radius(x::ArbReal{P}, err::ArbFloat{P}) where {P}
    err = -abs(err)
    ccall(@libarb(arb_add_error_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, err)
    return x
end

function decrease_radius(x::ArbReal{P}, err::ArbReal{P}) where {P}
    err = -abs(err)
    ccall(@libarb(arb_add_error_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, err)
    return x
end

decrease_radius(x::ArbReal{P}) where {P} = decrease_radius(x, ulp(x))


text:
# ------------ Factor ------------ #

@inline function (ctx::GenerateContext)(fn::typeof(factor), arg)
    ctx.score += arg
    ctx.weight += arg
    return arg
end


text:
# function Mvx(vx_pml::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dx::Float64, dt::Float64, ext::Int64, iflag::Int64)
#     rho = modExpand(rho, ext, iflag)
#     (m,n) = size(rho)
#     a1 = 9/8;   a2 = -1/24;
#     c1 = a1/dx; c2 = a2/dx;
#     C3 = zeros(m*n)
#     denum = zeros(m*n)
#     for ix = 1 : n
#         for iz = 1 : m
#             if ix < n
#                a = (rho[iz,ix+1]+rho[iz,ix]) / (2*dt)
#             else
#                a = rho[iz,ix] / dt
#             end
#             b = vx_pml[iz,ix] / 2
#             denum[(ix-1)*m+iz] = 1 / (a+b)
#             C3[(ix-1)*m+iz]= (a-b) / (a+b)
#         end
#     end
#     tmp = spzeros(n, n)
#     tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
#     for ix = 2: n-2
#         tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
#         tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
#     end
#     tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
#     tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
#     MvxBvx = spdiagm(C3)
#     MvxBp = spdiagm(denum) * kron(tmp, speye(m))
#     return MvxBvx, MvxBp
# end

function Mvx(vx_pml::SparseMatrixCSC{Float64,Int64}, dx::Float64, dt::Float64, ext::Int64, iflag::Int64)
    (m,n) = size(vx_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dx; c2 = a2/dx;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt + vx_pml[iz,ix]/2
            b = 1/dt - vx_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(n, n)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for ix = 2: n-2
        tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
        tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
    end
    tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
    tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
    MvxBvx = spdiagm(C3)
    MvxBp  = spdiagm(denum) * kron(tmp, speye(m))
    return MvxBvx, MvxBp
end

function Mvx_back(vx_pml::SparseMatrixCSC{Float64,Int64}, dx::Float64, dt::Float64, ext::Int64, iflag::Int64)
    (m,n) = size(vx_pml)
    a1 = 9/8;   a2 = -1/24;
    c1 = a1/dx; c2 = a2/dx;
    C3 = zeros(m*n)
    denum = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a = 1/dt - vx_pml[iz,ix]/2
            b = 1/dt + vx_pml[iz,ix]/2
            denum[(ix-1)*m+iz] = 1 / a
            C3[(ix-1)*m+iz]= b / a
        end
    end
    tmp = spzeros(n, n)
    tmp[1,1] = -c1; tmp[1,2] = c1; tmp[1,3] = c2;
    for ix = 2: n-2
        tmp[ix,ix  ] = -c1; tmp[ix, ix+1] = c1;
        tmp[ix,ix-1] = -c2; tmp[ix, ix+2] = c2;
    end
    tmp[n-1,n-1] = -c1; tmp[n-1, n  ] = c1; tmp[n-1,n-2] = -c2;
    tmp[n  ,n  ] = -c1; tmp[n  , n-1] =-c2;
    MvxBvx = spdiagm(C3)
    MvxBp  = -spdiagm(denum) * kron(tmp, speye(m))
    return MvxBvx, MvxBp
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
function petk05()
#
# M. Petkovic testing polynomials, p142
#
    p = [1.0;-7;20;-28;-18;110;-92;-44;345;225];
    z = [3.00000000000000 + 0.00000000000000im
         1.00000000000000 + 2.00000000000000im
         1.00000000000000 - 2.00000000000000im
        -1.00000000000000 - 0.00000000000000im];
    z = [z [2;2;2;3]];

    p, PolyZeros(z)
end


text:
#! /usr/bin/julia

# Make a pretty graphic for a Rosetta Code task

using Color, Images, FixedPointNumbers

const w = 10
const h = 16
const cw = 4w
const ch = 4h
const cnum = 11
const fw = (cnum + 1)*cw
const snum = 3
const hskip = 2h
const fh = snum*ch + (snum + 1)*hskip
const bh = 3h

chips = zeros(RGB{Ufixed8}, cnum, snum)
chips[:,1] = colormap("Reds", cnum)
chips[:,2] = colormap("Greens", cnum)
chips[:,3] = colormap("Blues", cnum)

a = fill(RGB{Ufixed8}(0.5, 0.5, 0.5), fh, fw)
img = Image(a)

bwr = w:(1+fw-w)
bhr = 1:bh
cwr = 1:cw
chr = 1:ch

for i in 1:snum
    hoff = ch + (i-1)*(ch + hskip)
    img["x", bwr, "y", bhr+hoff-bh] = RGB(1, 1, 1)
    img["x", bwr, "y", bhr+hoff] = RGB(0, 0, 0)
    for j in 1:cnum
        woff = 2w + (j-1)*cw
        img["x", cwr+woff, "y", chr+hoff-2h] = chips[j, i]
    end
end

imwrite(img, "color_chips.ppm")
imwrite(img, "color_chips.png")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code