text:
# # Baffled piston in a half-sphere domain

# Source code: [`baffled_piston_tut.jl`](baffled_piston_tut.jl)

# ## Description

# Baffled piston in a half-sphere domain with Absorbing Boundary Condition (ABC).
# This is a steady-state simulation.

# ![](baffled_piston.png)

# ## Goals

# - Demonstrate the use of an algorithm to run the simulation.

##
# ## Definitions

# The finite element code relies on the basic functionality implemented in this
# package.
using FinEtools
# The acoustics functionality is brought in:
using FinEtoolsAcoustics

# We shall need some facilities from the linear algebra package
using LinearAlgebra

# The input quantities are provided including the units:
rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
omega =  7500*phun("rev/s");      # frequency of the piston
a_piston =  -1.0*phun("mm/s")     # amplitude of the piston acceleration
R = 50.0*phun("mm");# radius of the piston
Ro = 150.0*phun("mm"); # radius of the enclosure
nref = 4;#number of refinements of the sphere around the piston
nlayers = 35;                     # number of layers of elements surrounding the piston
tolerance = R/(2^nref)/100

# Hexahedral mesh
fens,fes  =  H8sphere(R,nref);
bfes  =  meshboundary(fes)
File  =   "baffledabc_boundary.vtk"
vtkexportmesh(File, connasarray(bfes), fens.xyz, FinEtools.MeshExportModule.VTK.Q4)
@async run(`"paraview.exe" $File`)

l = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
fens1, fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
fens, newfes1, fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Piston surface mesh
bfes  =  meshboundary(fes)
l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
piston_fes = subset(bfes,intersect(l1,l2));

# Outer spherical boundary
louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
outer_fes = subset(bfes,louter);

# The simulation is driven by setting up and then executing an algorithm. For
# this purpose we store the data in a "model data dictionary".

# Region of the fluid. Define the finite element machine for the fluid region.
material = MatAcoustFluid(bulk, rho)
region1 =  FDataDict("femm"=>FEMMAcoust(IntegDomain(fes, GaussRule(3, 2)), material))

# Define the finite element machine for the ABC (absorbing boundary condition).
abc1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(outer_fes, GaussRule(2, 2)), material))

# The normal flux is prescribed on the surface of the piston in terms of the
# known acceleration.
flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(piston_fes, GaussRule(2, 2)), material),  "normal_flux"=> -rho*a_piston+0.0im);

# Make model data dictionary. It completely defines the parameters of the
# problem.
modeldata =  FDataDict("fens"=>fens, "omega"=>omega, "regions"=>[region1], "flux_bcs"=>[flux1], "ABCs"=>[abc1])

# Call the solver. The model data is returned enriched of the solution
# parameters, such as the pressure field `P`.
modeldata = FinEtoolsAcoustics.AlgoAcoustModule.steadystate(modeldata)

# Extract geometry field and the pressure field for postprocessing.
geom = modeldata["geom"]
P = modeldata["P"]

# Write out the postprocessing data is a VTK file. Dump both the magnitude and
# the components of the solution (complex field).
File  =   "baffledabc.vtk"
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.VTK.H8; scalars = [("absP", abs.(P.values)), ("realP", real.(P.values)), ("imagP", imag.(P.values))])
# If the `paraview` program is installed, run it.
@async run(`"paraview.exe" $File`)

true


text:
#=
This example demonstrates the use of callback functions working with a neighborhood struct.
It is a 2D advection-reaction equation.
It is solved in two ways:
 - u1 using a symbolic upwinding operator for the flux.
 - u2 using a callback function for the flux that uses the neighborhood of a face to do the same upwinding.

The execution time is also compared, but doesn't account for compilation/first run time.
Also, the number of flops using the symbolic version is much higher because it is actually doing
    1/2*dot(a,n)*(side1+side2) + (1-alpha)/2*abs(dot(a,n))*(side1-side2)
with a parameter alpha rather than
    dot(a,n) > 0 ? dot(a,n)*side1 : dot(a,n)*side2
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVneighborhood");

useLog("FVneighborhoodlog", level=3)

# Configuration setup
domain(2)
solverType(FV)

timeStepper(RK4)
n = 20
mesh(QUADMESH, elsperdim=n, bids=4)
# finiteVolumeOrder(2);

# Variables and BCs
u1 = variable("u1", location=CELL) # will use symbolic upwind operator
u2 = variable("u2", location=CELL) # will use callback with neighborhood

# This boundary provides a pulsing input through a small window on one side
boundary(u1, 1, FLUX, "(abs(y-0.2) < 0.11 && sin(6*pi*t)>0) ? 1 : 0") # x=0
boundary(u1, 2, NO_BC) # x=1
boundary(u1, 3, NO_BC) # y=0
boundary(u1, 4, NO_BC) # y=1
boundary(u2, 1, FLUX, "(abs(y-0.2) < 0.11 && sin(6*pi*t)>0) ? 1 : 0") # x=0
boundary(u2, 2, NO_BC) # x=1
boundary(u2, 3, NO_BC) # y=0
boundary(u2, 4, NO_BC) # y=1

# Time interval and initial condition
T = 1;
timeInterval(T)
initial(u1, "0")
initial(u2, "0")

coefficient("a", ["cos(pi*x/2)","sin(pi*x/2)"], type=VECTOR) # advection velocity
coefficient("s", "sin(pi*x)^4 * sin(pi*y)^4") # source

# Here is the callback function for the flux that simply does first order upwinding.
@callbackFunction(
    function flux_fun(neighborhood, a, normal)
        # neighborhood is a struct containing:
        # - neighborhood.cells -> an array of two arrays of integers like [[L1, L2,...], [R1, R2,...]]
        #                         where L1 is the nearest neighbor on the left, L2 is next-nearest, etc.
        # - neighborhood.centers -> similar, but holding cell center coordinates like [ [L1x L2x ...  , [R1x R2x ...
        #                                        *note these are 2d arrays->             L1y L2y ...]    R1y R2y ...] ]
        # - neighborhood.values -> holds the values of the requested variable(s) like [[u_L1, u_L2,...], [u_R1, u_R2,...]]
        #                          *vector valued variables or multiple variables will have columns of values

        a_dot_n = a[1]*normal[1] + a[2]*normal[2]; # a and normal are 2D vectors
        if a_dot_n > 0
            result = neighborhood.values[1][1] * a_dot_n; # use the first value on the left: values[1][1]
        else
            result = neighborhood.values[2][1] * a_dot_n; # use the first value on the right: values[2][1]
        end

        return result;
    end
)

# The PDE expression input
# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
flux(u1, "upwind(a, u1)")
source(u1, "0.1 * s")

flux(u2, "flux_fun(neighborhood(u2), a, normal())")
source(u2, "0.1 * s")

exportCode("fvneighborhoodcode") # uncomment to export generated code to a file
#importCode("fvneighborhoodcode") # uncomment to import code from a file

t1 = @elapsed(solve(u1))
t2 = @elapsed(solve(u2))

println("time for u1: "*string(t1))
println("time for u2: "*string(t2))

finalize_finch()

# uncomment to output to file
# output_values([u1,u2], "fvad2d", format="vtk");

##### Uncomment below to plot

# xy = Finch.fv_info.cellCenters

# using Plots
# pyplot();
# p1 = plot(xy[1,:], xy[2,:], u1.values[:], st=:surface)
# p2 = plot(xy[1,:], xy[2,:], u2.values[:], st=:surface)
# display(plot(p1, p2, layout=2))

ext: jl

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using Word2Vec
using LinearAlgebra

const stopWords = Set{String}([
    "a",
	"an",
	"and",
	"as",
	"based",
	"by",
	"for",
	"from",
	"in",
	"on",
	"of",
    "over",
	"that",
	"the",
	"this",
	"to",
	"via",
	"with",
	"without",
	"is",
	"are",
	"be",
	"you",
	"we",
	"they",
	"it",
	"your",
	"our",
	"their",
	"its",
	"what",
	"when",
	"where",
	"how",
	"do",
	"use",
	"here",
	"there",
	"using",
    "than",
	"、",
	",",
	"，",
	":",
	"：",
	".",
	"。",
	"‧",
	"!",
	"！",
	"?",
	"？",
	";",
	"；",
	"(",
	"（",
	")",
	"）",
	"'",
	"‘",
	"’",
	"\"",
	"「",
	"」",
	"“",
	"”",
	"`",
	"…",
    "</s>",
])

function getNonStopWords(wordVectors::WordVectors)::Tuple{Vector{String}, Dict{String,Int}}
    nonStopWords = Vector{String}()
    nonStopWordIDs = Dict{String,Int}()
    for (wordID, word) in enumerate(wordVectors.vocab)
        if !(word in stopWords)
            push!(nonStopWords, word)
            nonStopWordIDs[word] = wordID
        end
    end
    (nonStopWords, nonStopWordIDs)
end

function normalizeVectors(wordVectors::WordVectors, nonStopWords::Vector{String})::Array{Float64,2}
    numDims, numWords = size(wordVectors)
    numVectors = length(nonStopWords)
    result = Array{Float64,2}(undef, (numDims, numVectors))
    for i in 1:numVectors
        vector = get_vector(wordVectors, nonStopWords[i])
        w = 1.0 / norm(vector)
        result[:,i] = vector * w
    end
    result
end

struct IDPair
    i::Int
    j::Int
end

mutable struct Biterm
    pair::IDPair
    topicID::Int
    gpuFlag::Bool
end

struct NBTMWEParameters
    numTopics::Int
    maxIter::Int
    simThres::Float64
    alpha::Float64
    beta::Float64
    phrases::Vector{String}
    phraseIDs::Dict{String,Int}
    docs::Vector{Vector{Int}}
end

mutable struct NBTMWE
    parameters::NBTMWEParameters
    biterms::Vector{Biterm}
    similarPhraseIDs::Dict{IDPair, Vector{Int}}
    numBitermsOfTopic::Vector{Int}
    numRelatedOfTopic::Vector{Int}
    phraseFreqsOfTopic::Array{Float64,2}
end

function saveNBTMWE(fileName::String, btm::NBTMWE)
    file = open(fileName, "w")
    param = btm.parameters
    numPhrases = length(param.phrases)
    numDocs = length(param.docs)
    numBiterms = length(btm.biterms)
    numSimilar = length(btm.similarPhraseIDs)
    println(file, "$(param.numTopics) $(param.maxIter) $(param.simThres) $(param.alpha) $(param.beta) $numPhrases $numDocs $numBiterms $numSimilar")
    for phraseID in 1:numPhrases
        println(file, param.phrases[phraseID])
    end
    for docID in 1:numDocs
        doc = param.docs[docID]
        docSize = length(doc)
        for i in 1:docSize
            if i > 1
                print(file, " ")
            end
            print(file, doc[i])
        end
        print(file, "\n")
    end
    for bitermID in 1:numBiterms
        biterm = btm.biterms[bitermID]
        i = biterm.pair.i
        j = biterm.pair.j
        topicID = biterm.topicID
        gpuFlag = biterm.gpuFlag
        println(file, "$i $j $topicID $gpuFlag")
    end
    for (pair, relatedIDs) in btm.similarPhraseIDs
        print(file, "$(pair.i) $(pair.j)")
        for id in relatedIDs
            print(file, " $id")
        end
        print(file, "\n")
    end
    for topicID in 1:param.numTopics
        if topicID > 1
            print(file, " ")
        end
        print(file, btm.numBitermsOfTopic[topicID])
    end
    print(file, "\n")
    for topicID in 1:param.numTopics
        if topicID > 1
            print(file, " ")
        end
        print(file, btm.numRelatedOfTopic[topicID])
    end
    print(file, "\n")
    for topicID in 1:param.numTopics
        for phraseID in 1:numPhrases
            if phraseID > 1
                print(file, " ")
            end
            print(file, btm.phraseFreqsOfTopic[topicID,phraseID])
        end
        print(file, "\n")
    end
    close(file)
end

function newNBTMWE(fileName::String, nonStopWords::Vector{String}, vectors::Array{Float64,2}, numTopics::Int, maxIter::Int, simThres::Float64, alpha::Float64, beta::Float64)::NBTMWE
    phrases = nonStopWords
    phraseIDs = Dict{String,Int}()
    for (phraseID, phrase) in enumerate(phrases)
        phraseIDs[phrase] = phraseID
    end

    docs = Vector{Vector{Int}}()
    file = open(fileName)
    for line in readlines(file)
        phrasesAtLine = split(line)
        doc = Vector{Int}()
        for phrase in phrasesAtLine
            if haskey(phraseIDs, phrase)
                push!(doc, phraseIDs[phrase])
            end
        end
        if length(doc) > 0
            push!(docs, doc)
        end
    end

    numPhrases = length(phrases)
    numBitermsOfTopic = fill(0, numTopics)
    numRelatedOfTopic = fill(0, numTopics)
    phraseFreqsOfTopic = fill(0.0, (numTopics, numPhrases))
    biterms = Vector{Biterm}()
    for doc in docs
        docSize = length(doc)
        for i in 1:docSize
            for j in 1:docSize
                if i != j
                    topicID = rand(1:numTopics)
                    push!(biterms, Biterm(IDPair(doc[i], doc[j]), topicID, false))
                    numBitermsOfTopic[topicID] += 1
                    phraseFreqsOfTopic[topicID, doc[i]] += 1.0
                    phraseFreqsOfTopic[topicID, doc[j]] += 1.0
                end
            end
        end
    end

    similarPhraseIDs = Dict{IDPair, Vector{Int}}()
    numBiterms = length(biterms)
    println("building similar phrases list...")
    simMat = 0.5 .+ 0.5 .* (vectors' * vectors)
    for (bitermID, biterm) in enumerate(biterms)
        if haskey(similarPhraseIDs, biterm.pair)
            continue
        end
        relatedIDs = Vector{Int}()
        for phraseID in 1:numPhrases
            if phraseID == biterm.pair.i || phraseID == biterm.pair.j
                continue
            end
            sim1 = simMat[biterm.pair.i, phraseID]
            sim2 = simMat[biterm.pair.j, phraseID]
            if sim1 >= simThres && sim2 >= simThres
                push!(relatedIDs, phraseID)
            end
        end
        similarPhraseIDs[biterm.pair] = relatedIDs
        if bitermID % 10000 == 0
            println("similar phrase list $bitermID of $numBiterms")
        end
    end
    println("finish building similar phrases list")

    NBTMWE(NBTMWEParameters(numTopics, maxIter, simThres, alpha, beta, phrases, phraseIDs, docs), biterms, similarPhraseIDs, numBitermsOfTopic, numRelatedOfTopic, phraseFreqsOfTopic)
end

function loadNBTMWE(fileName::String)::NBTMWE
    file = open(fileName)
    firstLine = readline(file)
    fieldsOfFirstLine = split(firstLine)
    @assert length(fieldsOfFirstLine) == 9
    numTopics = parse(Int, fieldsOfFirstLine[1])
    maxIter = parse(Int, fieldsOfFirstLine[2])
    simThres = parse(Float64, fieldsOfFirstLine[3])
    alpha = parse(Float64, fieldsOfFirstLine[4])
    beta = parse(Float64, fieldsOfFirstLine[5])
    numPhrases = parse(Int, fieldsOfFirstLine[6])
    numDocs = parse(Int, fieldsOfFirstLine[7])
    numBiterms = parse(Int, fieldsOfFirstLine[8])
    numSimilar = parse(Int, fieldsOfFirstLine[9])
    phrases = Vector{String}()
    phraseIDs = Dict{String,Int}()
    for phraseID in 1:numPhrases
        phrase = readline(file)
        push!(phrases, phrase)
        phraseIDs[phrase] = phraseID
    end
    docs = Vector{Vector{Int}}()
    for docID in 1:numDocs
        doc = Vector{Int}()
        lineOfDoc = readline(file)
        fieldsOfDoc = split(lineOfDoc)
        docSize = length(fieldsOfDoc)
        for i in 1:docSize
            phraseID = parse(Int, fieldsOfDoc[i])
            push!(doc, phraseID)
        end
    end
    biterms = Vector{Biterm}(undef, numBiterms)
    for bitermID in 1:numBiterms
        lineOfBiterm = readline(file)
        fieldsOfBiterm = split(lineOfBiterm)
        @assert length(fieldsOfBiterm) == 4
        i = parse(Int, fieldsOfBiterm[1])
        j = parse(Int, fieldsOfBiterm[2])
        topicID = parse(Int, fieldsOfBiterm[3])
        gpuFlag = parse(Bool, fieldsOfBiterm[4])
        biterm = Biterm(IDPair(i, j), topicID, gpuFlag)
        biterms[bitermID] = biterm
    end
    similarPhraseIDs = Dict{IDPair, Vector{Int}}()
    for iSimilar in 1:numSimilar
        lineOfSimilar = readline(file)
        fieldsOfSimilar = split(lineOfSimilar)
        @assert length(fieldsOfSimilar) >= 2
        i = parse(Int, fieldsOfSimilar[1])
        j = parse(Int, fieldsOfSimilar[2])
        relatedIDs = Vector{Int}()
        for k in 3:length(fieldsOfSimilar)
            id = parse(Int, fieldsOfSimilar[k])
            push!(relatedIDs, id)
        end
        similarPhraseIDs[IDPair(i,j)] = relatedIDs
    end
    numBitermsOfTopic = Vector{Int}(undef, numTopics)
    lineOfNumBiterms = readline(file)
    fieldsOfNumBiterms = split(lineOfNumBiterms)
    @assert length(fieldsOfNumBiterms) == numTopics
    for topicID in 1:numTopics
        numBitermsOfTopic[topicID] = parse(Int, fieldsOfNumBiterms[topicID])
    end
    numRelatedOfTopic = Vector{Int}(undef, numTopics)
    lineOfNumRelated = readline(file)
    fieldsOfNumRelated = split(lineOfNumRelated)
    @assert length(fieldsOfNumRelated) == numTopics
    for topicID in 1:numTopics
        numRelatedOfTopic[topicID] = parse(Int, fieldsOfNumRelated[topicID])
    end
    phraseFreqsOfTopic = Array{Float64,2}(undef, (numTopics, numPhrases))
    for topicID in 1:numTopics
        lineOfTopic = readline(file)
        fieldsOfTopic = split(lineOfTopic)
        @assert length(fieldsOfTopic) == numPhrases
        for phraseID in 1:numPhrases
            freq = parse(Float64, fieldsOfTopic[phraseID])
            phraseFreqsOfTopic[topicID, phraseID] = freq
        end
    end
    close(file)

    NBTMWE(NBTMWEParameters(numTopics, maxIter, simThres, alpha, beta, phrases, phraseIDs, docs), biterms, phraseFreqsOfTopic, numBitermsOfTopic, numRelatedOfTopic, phraseFreqsOfTopic)
end

function computeProbTopicBiterm(btm::NBTMWE, biterm::Biterm)::Vector{Float64}
    numTopics = btm.parameters.numTopics
    numPhrases = length(btm.parameters.phrases)
    numBiterms = length(btm.biterms)
    numTotalRelated = sum(btm.numRelatedOfTopic)
    numBitermRelated = if biterm.gpuFlag
        length(btm.similarPhraseIDs[biterm.pair])
    else
        0
    end
    alpha = btm.parameters.alpha
    beta = btm.parameters.beta
    prob = Vector{Float64}(undef, numTopics)
    for topicID in 1:numTopics
        wRelated = 0.5 * btm.parameters.simThres * (numTotalRelated - numBitermRelated)
        if topicID != biterm.topicID
            w = 1.0 / (2 * btm.numBitermsOfTopic[topicID] + numPhrases * beta)
            prob1 = (btm.phraseFreqsOfTopic[topicID, biterm.pair.i] + beta) * w
            prob2 = (btm.phraseFreqsOfTopic[topicID, biterm.pair.j] + beta) * w
            wTopicRelated = 0.5 * btm.parameters.simThres * btm.numRelatedOfTopic[topicID]
            probTopic = (btm.numBitermsOfTopic[topicID] + wTopicRelated + alpha) / (numBiterms - 1 + wRelated + numTopics * alpha)
            prob[topicID] = prob1 * prob2 * probTopic
        else
            w = 1.0 / (2 * btm.numBitermsOfTopic[topicID] - 2 + numPhrases * btm.parameters.beta)
            prob1 = (btm.phraseFreqsOfTopic[topicID, biterm.pair.i] - 1 + btm.parameters.beta) * w
            prob2 = (btm.phraseFreqsOfTopic[topicID, biterm.pair.j] - 1 + btm.parameters.beta) * w
            wTopicRelated = 0.5 * btm.parameters.simThres * (btm.numRelatedOfTopic[topicID] - numBitermRelated)
            probTopic = (btm.numBitermsOfTopic[topicID] - 1 + wTopicRelated + alpha) / (numBiterms - 1 + wRelated + numTopics * alpha)
            prob[topicID] = prob1 * prob2 * probTopic
        end
    end
    c = 1.0 / sum(prob)
    prob .*= c
    prob
end

function sampleProb(prob::Vector{Float64})::Int
    x = rand()
    partialSum = 0.0
    for i in 1:length(prob)
        partialSum += prob[i]
        if partialSum >= x
            return i
        end
    end
    length(prob)
end

function sampleGPUFlag(prob::Vector{Float64}, topicID::Int)::Bool
    maxProb = max(prob...)
    ratio = prob[topicID] / maxProb
    x = rand()
    ratio > x
end

function trainNBTMWE!(btm::NBTMWE, fileName::String)
    numBiterms = length(btm.biterms)
    numCPUs = length(Sys.cpu_info())
    Segments = Vector{Tuple{Int,Int}}
    TopFlag = Tuple{Int, Bool}
    for iter in 1:btm.parameters.maxIter
        jobs = Channel{Tuple{Int,Int}}(numCPUs)
        jobOutputs = Channel{Tuple{Segments,Vector{Vector{TopFlag}}}}(numCPUs)

        function makeJobs(batchSize::Int)
            for i in 1:batchSize:numBiterms
                put!(jobs, (i, min(i + batchSize - 1, numBiterms)))
            end
        end

        function runJob(iCPU::Int)
            mySegments = Segments()
            myOutputs = Vector{Vector{TopFlag}}()
            for job in jobs
                push!(mySegments, job)
                jobOutput = Vector{TopFlag}(undef, job[2] - job[1] + 1)
                for bitermID in job[1]:job[2]
                    prob = computeProbTopicBiterm(btm, btm.biterms[bitermID])
                    newTopID = sampleProb(prob)
                    newFlag = sampleGPUFlag(prob, newTopID)
                    jobOutput[bitermID - job[1] + 1] = (newTopID,newFlag)
                end
                push!(myOutputs, jobOutput)
            end
            put!(jobOutputs, (mySegments, myOutputs))
        end

        bind(jobs, @async makeJobs(100))
        for iCPU in 1:numCPUs
            Threads.@spawn runJob(iCPU)
        end

        prevPhraseFreqsOfTopic = copy(btm.phraseFreqsOfTopic)
        simThres = btm.parameters.simThres
        for iCPU in 1:numCPUs
            segmentedOutputs = take!(jobOutputs)
            for (segID, segment) in enumerate(segmentedOutputs[1])
                segOutput = segmentedOutputs[2][segID]
                for bitermID in segment[1]:segment[2]
                    biterm = btm.biterms[bitermID]
                    topicFlag = segOutput[bitermID - segment[1] + 1]
                    newTopicID = topicFlag[1]
                    newFlag = topicFlag[2]
                    if newTopicID == biterm.topicID && newFlag == biterm.gpuFlag
                        continue
                    end
                    btm.numBitermsOfTopic[biterm.topicID] -= 1
                    btm.phraseFreqsOfTopic[biterm.topicID, biterm.pair.i] -= 1
                    btm.phraseFreqsOfTopic[biterm.topicID, biterm.pair.j] -= 1
                    if biterm.gpuFlag
                        relatedIDs = btm.similarPhraseIDs[biterm.pair]
                        btm.numRelatedOfTopic[biterm.topicID] -= length(relatedIDs)
                        for id in relatedIDs
                            btm.phraseFreqsOfTopic[biterm.topicID, id] -= simThres
                        end
                    end
                    btm.numBitermsOfTopic[newTopicID] += 1
                    btm.phraseFreqsOfTopic[newTopicID, biterm.pair.i] += 1
                    btm.phraseFreqsOfTopic[newTopicID, biterm.pair.j] += 1
                    btm.biterms[bitermID].topicID = newTopicID
                    if newFlag
                        relatedIDs = btm.similarPhraseIDs[biterm.pair]
                        btm.numRelatedOfTopic[newTopicID] += length(relatedIDs)
                        for id in relatedIDs
                            btm.phraseFreqsOfTopic[newTopicID, id] += simThres
                        end
                    end
                end
            end
        end
        deltaFreq = sum(abs.(btm.phraseFreqsOfTopic - prevPhraseFreqsOfTopic))
        println("iter $iter: deltaFreq = $deltaFreq")

        if iter % 100 == 0
            saveNBTMWE(fileName, btm)
            println("iter $iter: output saved")
        end
    end
    if maxIter % 100 != 0
        saveNBTMWE(fileName, btm)
    end
    println("final output saved")
end

numTopics = 100
maxIter = 300
simThres = 0.9
wordVectors = wordvectors("DualAspectTopicRepresentation/workspace/NIPS.w2v.model.txt")
nonStopWords, nonStopWordIDs = getNonStopWords(wordVectors)
vectors = normalizeVectors(wordVectors, nonStopWords)
btm = newNBTMWE("DualAspectTopicRepresentation/workspace/NIPS.w2v.training.txt", nonStopWords, vectors, numTopics, maxIter, simThres, 2.0, 0.005)
#saveBTM("DualAspectTopicRepresentation/workspace/NIPS.btm.txt", btm)
trainNBTMWE!(btm, "DualAspectTopicRepresentation/workspace/NIPS.nbtmwe.txt")


text:
using DSGE, Test

# Test switch
original = zeros(1,4)
default  = ones(1,4)

@testset "Test scenario switch" begin
    @test DSGE.switch(original, default, zeros(4), zeros(4))[1] == default
    @test DSGE.switch(original, default, ones(4), zeros(4))[1] == original
    @test DSGE.switch(original, default, ones(4), ones(4))[1] == [0. 1. 1. 1.]
end

# Test choose_last_period
@testset "Test scenario switch in the last period" begin
    @test DSGE.choose_switch_period(ones(3), 1) == 1
    @test DSGE.choose_switch_period(ones(3), 2) == 2
    @test DSGE.choose_switch_period(ones(3), 5) == 4
    @test DSGE.choose_switch_period(zeros(3), 1) == 4
end


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const FALLBACK_URL = "mosek://solve.mosek.com:30080"

using MosekTools
const optimizer = Mosek.Optimizer()
MOI.set(optimizer, MOI.RawOptimizerAttribute("fallback"), FALLBACK_URL)
MOI.set(optimizer, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Mosek"
end

@testset "Parameters" begin
    optimizer = Mosek.Optimizer()
    MOI.set(optimizer, MOI.RawOptimizerAttribute("fallback"), FALLBACK_URL)
    @testset "Double Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("INTPNT_CO_TOL_DFEAS"), 1e-7)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-7
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1e-8)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-8
        @testset "with integer value" begin
            MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1)
            @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1
        end
    end
    @testset "Integer Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS"), 100)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 100
        MOI.set(optimizer, MOI.RawOptimizerAttribute("INTPNT_MAX_ITERATIONS"), 200)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 200
        @testset "with enum value" begin
            MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_OPTIMIZER"), MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
            @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_IPAR_OPTIMIZER")) == convert(Int32, MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
        end
    end
    @testset "String Parameter" begin
        MOI.set(optimizer, MOI.RawOptimizerAttribute("PARAM_WRITE_FILE_NAME"), "foo.txt")
        # Needs https://github.com/JuliaOpt/Mosek.jl/pull/174
        #@test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "foo.txt"
        MOI.set(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME"), "bar.txt")
        #@test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "bar.txt"
    end
    @testset "TimeLimitSec" begin
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
        MOI.set(optimizer, MOI.TimeLimitSec(), 1.0)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == 1.0
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === 1.0
        MOI.set(optimizer, MOI.TimeLimitSec(), nothing)
        @test MOI.get(optimizer, MOI.RawOptimizerAttribute("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
    end
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(optimizer, false)
    @test MOIU.supports_default_copy_to(optimizer, true)
end

const config = MOIT.Config(
    Float64, atol=1e-3, rtol=1e-3,
    exclude=Any[MOI.ConstraintName, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus], # result in errors for now
)

@testset "Direct optimizer tests" begin
    MOIT.runtests(optimizer, config,
        exclude=[
            # issues/74,
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
            "test_linear_Indicator_ON_ONE",
            "test_model_ListOfConstraintAttributesSet", # TODO implement
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            "test_model_ModelFilter_AbstractVariableAttribute",
            "test_model_VariableName",
            "test_model_VariablePrimalStart",
            "test_model_copy_to_UnsupportedConstraint",
            "test_model_duplicate_VariableName",
            "test_modification_incorrect_VariableIndex",
            "test_solve_ObjectiveBound_MAX_SENSE_LP",
            "test_objective_set_via_modify",
            "test_model_VariableIndex_ConstraintName",
        ],
    )
end

@testset "Bridged and cached" begin
    model = MOIB.full_bridge_optimizer(Mosek.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)

    # linear and basic tests
    MOIT.runtests(model, config,
        exclude=[
            "test_basic_ScalarQuadraticFunction_EqualTo", # non-PSD quadratic
            "test_basic_ScalarQuadraticFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_Interval",
            "test_basic_ScalarQuadraticFunction_Semi",
            "test_basic_ScalarQuadraticFunction_ZeroOne",
            "test_basic_ScalarQuadraticFunction_Integer",
            "test_basic_VectorQuadraticFunction_Nonnegatives",
            "test_basic_VectorQuadraticFunction_Zeros",
            "test_basic_VectorQuadraticFunction_DualExponentialCone",
            "test_basic_VectorQuadraticFunction_DualPowerCone",
            "test_basic_VectorQuadraticFunction_ExponentialCone",
            "test_basic_VectorQuadraticFunction_GeometricMeanCone",
            "test_basic_VectorQuadraticFunction_LogDetConeTriangle",
            "test_basic_VectorQuadraticFunction_NormInfinityCone",
            "test_basic_VectorQuadraticFunction_NormNuclearCone",
            "test_basic_VectorQuadraticFunction_NormOneCone",
            "test_basic_VectorQuadraticFunction_NormSpectralCone",
            "test_basic_VectorQuadraticFunction_PositiveSemidefiniteConeSquare",
            "test_basic_VectorQuadraticFunction_PositiveSemidefiniteConeTriangle",
            "test_basic_VectorQuadraticFunction_PowerCone",
            "test_basic_VectorQuadraticFunction_RelativeEntropyCone",
            "test_basic_VectorQuadraticFunction_RootDetConeTriangle",
            "test_basic_VectorQuadraticFunction_RotatedSecondOrderCone",
            "test_basic_VectorQuadraticFunction_SecondOrderCone",
            "test_quadratic_nonconvex_constraint_basic",
            "test_quadratic_nonconvex_constraint_integration",
            "test_basic_VectorAffineFunction_PositiveSemidefiniteConeSquare", # AssertionError: (m.x_sd[ref2id(vi)]).matrix == -1 src/variable.jl:173
            "test_basic_VectorOfVariables_PositiveSemidefiniteConeSquare",
            "test_basic_VectorAffineFunction_LogDetConeTriangle",
            "test_basic_VectorAffineFunction_NormNuclearCone",
            "test_basic_VectorOfVariables_NormNuclearCone",
            "test_basic_VectorAffineFunction_NormSpectralCone",
            "test_basic_VectorOfVariables_NormSpectralCone",
            "test_basic_VectorAffineFunction_RootDetConeTriangle",
            "test_basic_VectorOfVariables_RootDetConeTriangle",
            "test_basic_VectorAffineFunction_PositiveSemidefiniteConeTriangle", # TODO: implement get ConstraintSet for SAF
            "test_basic_VectorOfVariables_PositiveSemidefiniteConeTriangle",
            "test_quadratic_SecondOrderCone_basic",
            "test_basic_VectorOfVariables_LogDetConeTriangle", # Mosek.MosekError(1307, "Variable '' (1) is a member of cone '' (0).") src/msk_functions.jl:477
            "test_conic_LogDetConeTriangle_VectorOfVariables",
            "test_constraint_ZeroOne_bounds", # Cannot put multiple bound sets of the same type on a variable
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
            "test_basic_ScalarQuadraticFunction_ZeroOne", # non-PSD quadratic
            "test_model_ListOfConstraintAttributesSet", # list not properly set
            "BoundAlreadySet", # TODO throw error if bound already set
            "test_model_ModelFilter_AbstractVariableAttribute",
            "test_model_VariableName", # Mosek currently throws when setting twice, not when getting names
            "test_model_Name_VariableName_ConstraintName",
            "test_model_duplicate_VariableName",
            "test_model_VariablePrimalStart", # able to set but not to get VariablePrimalStart
            "test_objective_qp_ObjectiveFunction_zero_ofdiag", # MOI.ListOfModelAttributesSet
            "test_objective_set_via_modify",
            "test_quadratic_nonconvex_constraint_integration",
            "test_solve_ObjectiveBound_MAX_SENSE_LP", # ObjectiveBound invalid
        ],
    )

    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)

    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                Mosek.Optimizer(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)

    MOI.Test.runtests(model, config,
        exclude=[
            "test_basic_ScalarQuadraticFunction_EqualTo", # non-PSD quadratic
            "test_basic_ScalarQuadraticFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_Interval",
            "test_basic_ScalarQuadraticFunction_Semi",
            "test_basic_ScalarQuadraticFunction_ZeroOne",
            "test_basic_ScalarQuadraticFunction_Integer",
            "test_basic_VectorQuadraticFunction_Nonnegatives",
            "test_basic_VectorQuadraticFunction_Zeros",
            "test_quadratic_nonconvex_constraint_basic",
            "test_quadratic_nonconvex_constraint_integration",
            "test_quadratic_SecondOrderCone_basic",
            "test_basic_VectorQuadraticFunction_", # not PSD because of equality
            "test_basic_VectorOfVariables_LogDetConeTriangle", # Mosek.MosekError(1307, "Variable '' (1) is a member of cone '' (0).") src/msk_functions.jl:477
            "test_conic_LogDetConeTriangle_VectorOfVariables",
            "test_constraint_ZeroOne_bounds",
            "BoundAlreadySet", # TODO throw error if bound already set
            "test_solve_ObjectiveBound_MAX_SENSE_LP",
            "test_variable_solve_ZeroOne_with_0_upper_bound",
            "test_variable_solve_ZeroOne_with_upper_bound",
        ],
    )
end


text:
#!/usr/bin/env julia

using LogMath, Base.Test

include("test-math.jl")
include("test-LogNum.jl")


text:
#!/usr/bin/env julia

using Base.Test
using Distributions
using ForwardDiff
using StaticArrays

using Celeste: Model, DeterministicVI

import Celeste: DeterministicVI, ParallelRun, Log
import Celeste: PSF, SDSSIO, SensitiveFloats, Transform
import Celeste.SDSSIO: RunCamcolField

include(joinpath(Pkg.dir("Celeste"), "test", "SampleData.jl"))

using SampleData

# Set logging level and timing reporting (TODO: ability to set on command line)
Log.LEVEL[] = Log.WARN
verbose_timing = false

# Check whether to run time-consuming tests.
long_running_flag = "--long-running"
test_long_running = long_running_flag in ARGS

test_files = setdiff(ARGS, [ long_running_flag ])
if length(test_files) > 0
    testfiles = ["test_$(arg).jl" for arg in test_files]
else
    testdir = joinpath(Pkg.dir("Celeste"), "test")
    testfiles = filter(r"^test_.*\.jl$", readdir(testdir))
end

timing_info = Any[]
for testfile in testfiles
        _, t, bytes, gctime, memallocs = @timed include(testfile)
        push!(timing_info, (t, bytes, gctime, memallocs))
end

println("\nTiming info:")
totaltime = 0.0
for i in eachindex(timing_info)
    t, bytes, gctime, memallocs = timing_info[i]
    totaltime += t
    @printf "%30s: " testfiles[i]
    if verbose_timing
        Base.time_print(1e9 * t, memallocs.allocd, memallocs.total_time,
                        Base.gc_alloc_count(memallocs))
    else
        @printf "%7.2f seconds\n" t
    end
end
@printf "%30s: %7.2f seconds\n" "Total time" totaltime


text:
"""
    ComtradeCfg

`ComtradeCfg` holds information from a COMTRADE "cfg" configuration file. Fields
include:

* `station_name::String` -- Identifier for the substation name.
* `rec_dev_id::String` -- Identifier for the recording device.
* `rev_year::Int` -- File format revision, normally 1991, 1999, or 2013.
* `tt::Int` -- Total number of channels.
* `nA::Int` -- Number of analog channels.
* `nD::Int` -- Number of digital channels.
* `A::DataFrame` -- Table of analog channel information, one row per channel.
* `D::DataFrame` -- Table of digital channel information, one row per channel.
* `lf::Float64` -- Nominal line frequency in hertz.
* `nrates::Int` -- Number of sample rates.
* `samp::Vector{Float64}` -- Sample rate in hertz.
* `endsamp::Vector{Int}` -- Numper of sample points at each sampling rate.
* `npts::Int` -- Number of sample points.
* `time` -- Time stamp of the first data point.
* `triggertime` -- Time stamp of the trigger time.
* `ft::ASCIIString` -- File type, options include "ASCII" or "BINARY".
* `timemult::Float64` -- Multiplication factor for the timestamp.
* `time_code::ASCIIString` -- Time difference between local time and UTC, including daylight savings (2013).
* `local_code::ASCIIString` -- Time difference between local time and UTC when time stamps are in UTC (2013).
* `tmq_code::ASCIIString` -- Time quality as defined in IEEE C37.118. F=not reliable, 4=time within 1 μs, 0=clock locked in (2013).
* `leapsec::Int` -- Leap second indicator, 0=no adjustment, 1=a leap second was added, 2=a leap second was subtracted, 3=no capability to address (2013).

`A` includes the following columns:

* `:An::Int` -- Analog channel index.
* `:ch_id::String` -- Channel identifier.
* `:ph::String` -- Channel phase identifier.
* `:ccbm::String` -- Circuit component.
* `:uu::String` -- Units.
* `:a::Float64` -- Channel multiplier.
* `:b::Float64` -- Channel offset adder.
* `:skew::Float64` -- Time skew in microseconds from the start of the sample period.
* `:min::Float64` -- Minimum data value.
* `:max::Float64` -- Maximum data value.
* `:primary::Float64` -- Voltage or current transformer ratio primary factor.
* `:secondary::Float64` -- Voltage or current transformer ratio secondary factor.
* `:PS::String` -- "P" means `ax+b` is a primary value; "S" means a secondary value.

`D` includes the following columns:

* `:Dn::Int` -- Digital channel index.
* `:ch_id::String` -- Channel identifier.
* `:ph::String` -- Channel phase identifier.
* `:ccbm::String` -- Circuit component.
* `:y::Int` -- Normal state, 0 or 1.

Note that the `:PS` column of `A` is not used to adjust the data in `read_comtrade`.
Likewise, `:skew` is not used to adjust the time stamp.

Some of the columns in `A` and `D` may be missing in older file revisions.
"""

using CSV

mutable struct ComtradeCfg
    station_name::String
    rec_dev_id::String
    rev_year::Int
    tt::Int
    nA::Int
    nD::Int
    A::DataFrame
    D::DataFrame
    lf::Float64
    nrates::Int
    samp::Vector{Float64}
    endsamp::Vector{Int}
    npts::Int
    time
    triggertime
    ft::String
    timemult::Float64
    time_code::String
    local_code::String
    tmq_code::String
    leapsec::Int
end

ComtradeCfg() = ComtradeCfg("", "", 1991, 0, 0, 0,
                            DataFrame(),
                            DataFrame(Any[Int[] for i in 1:4], [:Dn, :ch_id, :ph, :ccbm]),
                            0.0, 1, [0.0], [0], 0,
                            0, 0, "", 1.0,
                            "", "", "", 0)
asint(x) = parse(Int, x)
asfloat(x) = parse(Float64, x)

function read_cfg(fn)
    x = ComtradeCfg()
    a = readlines(fn)
    sa = split(a[1], ',')
    x.station_name = sa[1]
    if length(sa) > 1
        x.rec_dev_id = sa[2]
    end

    if length(sa) > 2
        x.rev_year = asint(sa[3])
    end
    sa = split(a[2], ',')
    x.tt = asint(sa[1])
    x.nA = asint(replace(sa[2], "A"=> ""))
    x.nD = asint(replace(sa[3], "D"=>""))
    x.A = CSV.read(fn, delim = ',', header = false, skipto = 3,
                    limit = x.nA)
    rename!(x.A, [:An, :ch_id, :ph, :ccbm, :uu, :a, :b, :skew, :min, :max, :primary, :secondary, :PS][1:ncol(x.A)])
    if x.nD > 0
        x.D = CSV.read(fn, delim = ',', header = false, skipto = 3 + x.nA,
                        limit = x.nD)
        if x.rev_year > 1991
            rename!(x.D, [:Dn, :ch_id, :ph, :ccbm, :y])
        else
            rename!(x.D, [:Dn, :ch_id, :y])
        end
    end
    x.nrates  = asint(a[x.tt + 4])
    for i in 1:max(x.nrates, 1)
        sa = split(a[x.tt + 4 + i], ',')
        x.samp[i]    = asfloat(sa[1])
        x.endsamp[i] = asint(sa[2])
    end
    nrow = x.tt + max(x.nrates, 1)
    x.npts = sum(x.endsamp)
    try
        x.time        = DateTimeMicro(strip(a[nrow + 5]))
    catch
        x.time        = strip(a[nrow + 5])
    end
    try
        x.triggertime = DateTimeMicro(strip(a[nrow + 6]))
    catch
        x.triggertime = strip(a[nrow + 6])
    end
    x.ft     = uppercase(strip(a[nrow + 7]))
    if x.rev_year >= 1999
        x.timemult = asfloat(a[nrow + 8])
    end
    if x.rev_year >= 2013
        sa = split(a[nrow + 9], ',')
        x.time_code = sa[1]
        x.local_code = strip(sa[2])
        sa = split(a[nrow + 10], ',')
        x.tmq_code = sa[1]
        x.leapsec = asint(sa[2])
    end
    return x
end


text:
# Methods for HMM for parsing reads of sequences.

"""
Column index for each nucleotide in observation matrix.
"""
const NUCLEOTIDE_COLS = Dict('A' => UInt8(1),
                             'C' => UInt8(2),
                             'G' => UInt8(3),
                             'T' => UInt8(4));

"""
    viterbi_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})

Return viterbi path and log probability for that path. Takes logs of matrices.
`observations_given_states` has # rows = # states, # columns = # steps of markov process.
"""
function viterbi_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})
    numstates, numsteps = size(observations_given_states)
    scores = zeros(numstates, numsteps)
    paths = zeros(scores, UInt16)
    # Forward Score propogation
    scores[:, 1] = initials + observations_given_states[:, 1]
    for t in 2:numsteps
        for j in 1:numstates
            transitioned_scores = scores[:, t-1] + transitions[:, j]
            paths[j, t] = indmax(transitioned_scores)
            scores[j, t] = transitioned_scores[paths[j, t]] + observations_given_states[j, t]
        end
    end
    # Backtrack
    path = zeros(Int, numsteps)
    probs = zeros(path, Float64)
    path[end] = indmax(scores[:, end])
    probs[end] = scores[path[end], end]
    for t in numsteps-1:-1:1
        path[t] = paths[path[t+1], t+1]
        probs[t] = scores[path[t+1], t+1]
    end
    return path, maximum(scores[:, end])
end

"""
    trans_mat(; uniform_cycle_prob = 0.9999999999, homopoly_cycle_prob = 0.98)

Return 5x5 transition matrix with given transition probabilities
State 1 -- uniform observation distribution.
State 2 -- High "A" observation likelihood.
State 3 -- High "C" observation likelihood.
State 4 -- High "G" observation likelihood.
State 5 -- High "T" observation likelihood.
Each state has high likelihood to return to state 1.
States 2 through 5 represent high likelihood of homopolymer regions of respective nuc."""
function trans_mat(; uniform_cycle_prob = 0.9999999999, homopoly_cycle_prob = 0.98)
    uh_prob = (1-uniform_cycle_prob) / 4
    hu_prob = 1 - homopoly_cycle_prob
    T = [
        uniform_cycle_prob uh_prob uh_prob uh_prob uh_prob;

        hu_prob homopoly_cycle_prob 0 0 0;

        hu_prob 0 homopoly_cycle_prob 0 0;

        hu_prob 0 0 homopoly_cycle_prob 0;

        hu_prob 0 0 0 homopoly_cycle_prob;
    ]
    return T
end

"""
    obs_mat(; homopoly_prob = 0.99)

Return 5x4 observation matrix with given probabilities.
State 1 -- uniform observation distribution.
State 2 -- High "A" observation likelihood.
State 3 -- High "C" observation likelihood.
State 4 -- High "G" observation likelihood.
State 5 -- High "T" observation likelihood.
Columns are in same order as columns 2 through 5
(the last four rows give a symmetric matrix).
"""
function obs_mat(; homopoly_prob = 0.99)
    uniform_prob = 0.25
    not_hpoly_prob = (1-homopoly_prob) / 3
    O = [
        uniform_prob uniform_prob uniform_prob uniform_prob;

        homopoly_prob  not_hpoly_prob not_hpoly_prob not_hpoly_prob;

        not_hpoly_prob homopoly_prob  not_hpoly_prob not_hpoly_prob;

        not_hpoly_prob not_hpoly_prob homopoly_prob  not_hpoly_prob;

        not_hpoly_prob not_hpoly_prob not_hpoly_prob homopoly_prob;
    ]
    return O
end

"""
    initial_dist(; uniform_state = 0.99)

Initial state distributions: see `trans_mat` for state descriptions. 5x1 vector.
"""
function initial_dist(; uniform_state = 0.99)
    not_unistate = (1 - uniform_state) / 4
    return [uniform_state; not_unistate; not_unistate; not_unistate; not_unistate]
end

"""
    get_obs_given_state(observation_matrix::Array{Float64,2}, observation_seq::String)

Populate 5xT matrix with likelihood of observation at each of T time steps given each
nucleotide in `observation_seq`.
"""
function get_obs_given_state(observation_matrix::Array{Float64,2}, observation_seq::String)
    obs_given_state = zeros(length(observation_seq), 5)
    for (i, nuc) in enumerate(observation_seq)
        obs_given_state[i, :] = observation_matrix[:, NUCLEOTIDE_COLS[nuc]]
    end
    return transpose(obs_given_state)
end

homopolymer_filter(seqs::Array{String}) = homopolymer_filter(seqs, [], [])[1]
markov_filter = homopolymer_filter

"""
    homopolymer_filter(seqs::Array{String,1}, phreds, names;
                       transmat = nothing, obsmat = nothing,
                       initialdist = nothing)

Filter sequences with "bad" sections in the middle -- abnormally long runs of a single base, using
viterbi alg inference. If this homopolymer occurs on one end of a sequence, keeps sequence and trims
homopolymer region off end. `phreds` and/or `names` may be `nothing`, in which case `nothing` is returned
for the respective field. If transition, observation, initial distribution matrices not provided (default)
then the default values from the respective constructors are used.
"""
function homopolymer_filter(seqs::Array{String,1}, phreds, names;
                       transmat = nothing, obsmat = nothing,
                       initialdist = nothing)
    newseqs = String[]
    # return empty arrays for newphreds, newnames if params were empty arrays
    # otherwise if params were nothing, return nothings
    newphreds= phreds==nothing ? nothing : Vector{Phred}[]
    newnames = names==nothing ? nothing : String[]
    phredsexist = (phreds != nothing) && (length(phreds) == length(seqs))
    namesexist = (names != nothing) && (length(names) == length(seqs))
    if (!phredsexist && phreds != nothing && length(phreds) > 0) || (!namesexist && names != nothing && length(names) > 0)
        error("Dimension mismatch in names or phreds array")
    end
    if transmat == nothing
        transmat = trans_mat()
    end
    if obsmat == nothing
        obsmat = obs_mat()
    end
    if initialdist == nothing
        initialdist = initial_dist()
    end
    T_logs, O_logs, i_logs = log.(transmat), log.(obsmat), log.(initialdist)
    for i in 1:length(seqs)
        seq = seqs[i]
        path = viterbi_logs(get_obs_given_state(O_logs, seq), T_logs, i_logs)[1]
        # remove bad parts
        lo = 1
        while lo < length(path) && path[lo] != 1
            lo += 1
        end
        hi = length(path)
        while hi > 0 && path[hi] != 1
            hi -= 1
        end
        if lo != 1 || hi != length(path)
            path = path[lo:hi]
        end
        # make sure no bad parts in middle
        ind = findnext(x->x != 1, path, lo)
        if ind > hi || ind == 0
            push!(newseqs, seq[lo:hi])
            if phredsexist
                push!(newphreds, phreds[i][lo:hi])
            end
            if namesexist
                push!(newnames, names[i])
            end
        end
    end
    # may return empty arrays for newphreds and newnames
    return newseqs, newphreds, newnames
end

"""
    homopolymer_filter(sourcepath::String, destpath::String;
                       transmat = nothing, obsmat = nothing,
                       initialdist = nothing, format="fastq")

Filter sequences with "bad" sections in the middle -- abnormally long runs of a single base, using
viterbi alg inference. If this homopolymer occurs on one end of a sequence, keeps sequence and trims
homopolymer region off end. Takes a file path for each of a source file of type `format` (which may be "fasta" or "fastq")
and a destination file is written which is the same file type.
"""
function homopolymer_filter(sourcepath::String, destpath::String;
                       transmat = nothing, obsmat = nothing,
                       initialdist = nothing, format="fastq")
    seqs, phreds, names = String[], Vector{Phred}[], String[]
    if format == "fastq"
        seqs, phreds, names = read_fastq(sourcepath, seqtype=String)
    elseif format == "fasta"
        names, seqs = read_fasta_with_names(sourcepath, seqtype=String)
    else
        println("Invalid format! File not written.")
    end
    println("Cleaning $(length(seqs)) sequences")
    newseqs, newphreds, newnames = homopolymer_filter(seqs, phreds, names,
                                                 transmat=transmat, obsmat=obsmat,
                                                 initialdist=initialdist)
    println("Writing $(length(newseqs)) sequences")
    if format == "fastq"
        write_fastq(destpath, newseqs, newphreds, names=newnames, DNASeqType=false)
    else
        write_fasta(destpath, newseqs, names=newnames, DNASeqType=false)
    end
end

"""
    forward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})

Compute logs of forward scores. Takes logs of matrices as inputs.
`observations_given_states` has # rows = # states, # columns = # steps of markov process.
"""
function forward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})
    numstates, numsteps = size(observations_given_states)
    fwds = zeros(numstates, numsteps)
    # Forward pass - (unnormalized?) probs of states at time t given observations up to t
    fwds[:, 1] = initials + observations_given_states[:, 1]
    for t in 2:numsteps
        for j in 1:numstates
            prob_state_given_prev = reduce(logsum, fwds[:, t-1] + transitions[:, j])
            fwds[j, t] = prob_state_given_prev + observations_given_states[j, t]
        end
    end
    return fwds
end

"""
    backward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})

Compute logs of backward scores and individual posterior probabilities. Takes logs of matrices as inputs.
`observations_given_states` has # rows = # states, # columns = # steps of markov process.
"""
function backward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})
    numstates, numsteps = size(observations_given_states)
    bwds = zeros(numstates, numsteps)
    # Backwards pass - (unnormalized?) probs of states at time t given observations after t
    bwds[:, end] = 1  # arbitrarily initialized
    log_posteriors = zeros(numstates, numsteps)
    for t in numsteps-1:-1:1
        for i in 1:numstates
            bwds[i, t] = reduce(logsum, bwds[:, t+1] + transitions[i, :] + observations_given_states[:, t+1])
        end
    end
    return bwds
end

"""
    forward_backward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})

Compute logs of forward-backward scores and individual probabilities. Takes logs of matrices as inputs.
`observations_given_states` has # rows = # states, # columns = # steps of markov process.
"""
function forward_backward_logs(observations_given_states::Array{Float64, 2}, transitions::Array{Float64, 2}, initials::Array{Float64})
    numstates, numsteps = size(observations_given_states)
    fwds = forward_logs(observations_given_states, transitions, initials)
    bwds = zeros(numstates, numsteps)
    # Backwards pass - Compute manually in order to compute individual probabilities simultaneously
    # - saves an extra pass through observations arrays
    bwds[:, end] = 1  # arbitrarily initialized
    log_posteriors = zeros(numstates, numsteps)
    for t in numsteps-1:-1:1
        for i in 1:numstates
            bwds[i, t] = reduce(logsum, bwds[:, t+1] + transitions[i, :] + observations_given_states[:, t+1])
        end
        # and combine with fwds to get posterior probability logs
        fb_probs = fwds[:, t] + bwds[:, t]
        norm = reduce(logsum, fb_probs)
        log_posteriors[:, t] = fb_probs - norm
    end
    # and for the end
    fb_probs = fwds[:, end] + bwds[:, end]
    norm = reduce(logsum, fb_probs)
    log_posteriors[:, end] = fb_probs - norm
    return log_posteriors
end

"""
    gen_seq_with_model(n::Int, trans_mat, obs_mat, initial_dists)

Generate a sequence given a markov model. May create bad reads (long runs of a single base).
"""
function gen_seq_with_model(n::Int, trans_mat, obs_mat, initial_dists)
    state_seq = zeros(Int, n)
    state_seq[1] = wsample(collect(1:5), trans_mat*initial_dists)
    obs_seq = wsample(["A", "C", "G", "T"], obs_mat[state_seq[1], :])
    for t in 2:n
        state_seq[t] = wsample(collect(1:5), trans_mat[state_seq[t-1], :])
        obs_seq *= wsample(["A", "C", "G", "T"], obs_mat[state_seq[t], :])
    end
    return obs_seq
end

"""
    viterbiprint(s::String)

Draw flagged sites with capital letters, safe sites with lowercase.
"""
function viterbiprint(s::String)
    obs = get_obs_given_state(obs_mat(), s)
    v = viterbi_logs(log.(obs),log.(trans_mat()),log.(initial_dist()))
    v = v[1]
    s = collect(s)
    for i in 1:length(v)
        if v[i]==1
            s[i] += 32
        end
    end
    return join(s)
end


text:
function shaker!(A)
    a = firstindex(A)
    b = length(A)
    while (a <= b)
        for i in b:a+1
            if A[i-1]>A[i]
                A[i-1],A[i]=A[i],A[i-1]
            end
        end
        a=a+1
        for j in a:b+
            if A[j-1]>A[j]
                A[j-1],A[j]=A[j],A[j-1]
            end
        end
        b=b-1
    end
    return A
end


text:
using NHST

bins = [3, 4, 8]
expected_bins = [8, 4, 3]

@assert isapprox(p_value(run(ChisqTest, bins)), 0.246597)
@assert isapprox(test_statistic(run(ChisqTest, bins)), 2.8)

@assert isapprox(p_value(ChisqTest, bins), 0.246597)
@assert isapprox(test_statistic(ChisqTest, bins), 2.8)

@assert isapprox(p_value(run(ChisqTest, bins, expected_bins)), 0.003249784)
@assert isapprox(test_statistic(run(ChisqTest, bins, expected_bins)), 11.45833)

@assert isapprox(p_value(ChisqTest, bins, expected_bins), 0.003249784)
@assert isapprox(test_statistic(ChisqTest, bins, expected_bins), 11.45833)


text:
CholType{T,S<:AbstractMatrix} = Cholesky{T,S}
chol_lower(a::Matrix) = cholesky(a).L

# always use the stored cholesky factor, not a copy
chol_lower(a::CholType) = a.uplo === 'L' ? a.L : a.U'
chol_upper(a::CholType) = a.uplo === 'U' ? a.U : a.L'

if HAVE_CHOLMOD
    CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf.L
end


text:
@doc raw"""
    StochasticGradientProblem <: Problem

A stochastic gradient problem consists of
* a `Manifold M`
* a(n optional) cost function ``f(x) = \displaystyle\sum_{i=1}^n f_i(x)
* an array of gradients, i.e. a function that returns and array or an array of functions ``\{∇f_i\}_{i=1}^n``.

# Constructors
    StochasticGradientProblem(M::Manifold, ∇::Function; cost=Missing())
    StochasticGradientProblem(M::Manifold, ∇::AbstractVector{<:Function}; cost=Missing())

Create a Stochastic gradient problem with an optional `cost` and the gradient either as one
function (returning an array) or a vector of functions.
"""
struct StochasticGradientProblem{MT<:Manifold,TCost,TGradient} <: Problem
    M::MT
    cost::TCost
    ∇::TGradient
end
function StochasticGradientProblem(
    M::TM, ∇::Function; cost::Union{Function,Missing}=Missing()
) where {TM<:Manifold}
    return StochasticGradientProblem{TM,typeof(cost),Function}(M, cost, ∇)
end
function StochasticGradientProblem(
    M::TM, ∇::AbstractVector{<:Function}; cost::Union{Function,Missing}=Missing()
) where {TM<:Manifold}
    return StochasticGradientProblem{TM,typeof(cost),typeof(∇)}(M, cost, ∇)
end

@doc raw"""
    get_gradients(P::StochasticGradientProblem, x)

Evaluate all summands gradients ``\{∇f_i\}_{i=1}^n`` at `x`.
"""
get_gradients(P::StochasticGradientProblem{<:Manifold,TC,<:Function}, x) where {TC} = P.∇(x)
function get_gradients(
    P::StochasticGradientProblem{<:Manifold,TC,<:AbstractVector}, x
) where {TC}
    return [∇i(x) for ∇i in P.∇]
end

@doc raw"""
    get_gradient(P::StochasticGradientProblem, k, x)

Evaluate one of the summands gradients ``∇f_k``, ``k\in \{1,…,n\}``, at `x`.
"""
function get_gradient(
    P::StochasticGradientProblem{<:Manifold,TC,<:Function}, k, x
) where {TC}
    return P.∇(x)[k]
end
function get_gradient(
    P::StochasticGradientProblem{<:Manifold,TC,<:AbstractVector}, k, x
) where {TC}
    return P.∇[k](x)
end

"""
    AbstractStochasticGradientDescentOptions <: Options

A generic type for all options related to stochastic gradient descent methods
"""
abstract type AbstractStochasticGradientProcessor <: DirectionUpdateRule end

"""
    StochasticGradientDescentOptions <: AbstractStochasticGradientDescentOptions

Store the following fields for a default stochastic gradient descent algorithm,
see also [`StochasticGradientProblem`](@ref) and [`stochastic_gradient_descent`](@ref).

# fields

# Fields
* `x` the current iterate
* `stopping_criterion` ([`StopAfterIteration`](@ref)`(1000)`)– a [`StoppingCriterion`](@ref)
* `stepsize` ([`ConstantStepsize`](@ref)`(1.0)`) a [`Stepsize`](@ref)
* `evaluation_order` – (`:Random`) – whether
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle permuted sequence (`:Linear`) or the default `:Random` one.
* `order` the current permutation
* `retraction_method` – (`ExponentialRetraction()`) a `retraction(M,x,ξ)` to use.

# Constructor
    StochasticGradientDescentOptions(x)

Create a [`StochasticGradientDescentOptions`](@ref) with start point `x`.
all other fields are optional keyword arguments.
"""
mutable struct StochasticGradientDescentOptions{
    TX,
    D<:DirectionUpdateRule,
    TStop<:StoppingCriterion,
    TStep<:Stepsize,
    RM<:AbstractRetractionMethod,
} <: AbstractGradientDescentOptions
    x::TX
    direction::D
    stop::TStop
    stepsize::TStep
    order_type::Symbol
    order::Vector{<:Int}
    retraction_method::RM
    k::Int # current iterate
end
function StochasticGradientDescentOptions(
    x;
    direction::DirectionUpdateRule=StochasticGradient(),
    order_type::Symbol=:RandomOrder,
    order::Vector{<:Int}=Int[],
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    stoping_criterion::StoppingCriterion=StopAfterIteration(1000),
    stepsize::Stepsize=ConstantStepsize(1.0),
)
    return StochasticGradientDescentOptions{
        typeof(x),
        typeof(direction),
        typeof(stoping_criterion),
        typeof(stepsize),
        typeof(retraction_method),
    }(
        x, direction, stoping_criterion, stepsize, order_type, order, retraction_method, 0
    )
end

"""
    StochasticGradient <: DirectionUpdateRule

The default gradient processor, which just evaluates the (stochastic) gradient or a subset
thereof.
"""
struct StochasticGradient <: AbstractStochasticGradientProcessor end
function (s::StochasticGradient)(
    p::StochasticGradientProblem, o::StochasticGradientDescentOptions, iter
)
    # for each new epoche choose new order if we are at random order
    ((o.k == 1) && (o.order_type == :Random)) && shuffle!(o.order)
    # i is the gradient to choose, either from the order or completely random
    j = o.order_type == :Random ? rand(1:length(o.order)) : o.order[o.k]
    return o.stepsize(p, o, iter), get_gradient(p, j, o.x)
end
function MomentumGradient(
    p::StochasticGradientProblem,
    x0::P,
    s::DirectionUpdateRule=StochasticGradient();
    ∇=zero_tangent_vector(p.M, x0),
    momentum=0.2,
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(∇),typeof(momentum),VTM}(
        deepcopy(x0), ∇, momentum, s, vector_transport_method
    )
end
function AverageGradient(
    p::StochasticGradientProblem,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=StochasticGradient();
    gradients=fill(zero_tangent_vector(p.M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, deepcopy(x0), s, vector_transport_method
    )
end


text:
using TerminalUI; using VT100

const thisdir = dirname(@__FILE__)

function _compare(output, outbuf)
    result = outbuf == output
    if !result
        println("Test failed. Expected result written to expected.out,
            actual result written to failed.out")
        open("failed.out","w") do f
            write(f,outbuf)
        end
        open("expected.out","w") do f
            write(f,output)
        end
    end
    return result
end

function compare(em, output, decorator = nothing)
    buf = IOBuffer()
    decoratorbuf = IOBuffer()
    VT100.dump(buf,decoratorbuf,em)
    outbuf = takebuf_array(buf)
    decoratorbuf = takebuf_array(decoratorbuf)
    println(String(decoratorbuf))
    _compare(output, outbuf) || return false
    decorator === nothing && return true
    _compare(decorator, decoratorbuf)
end

function test_against(em,file)
    output = open(read,joinpath(thisdir,file))
    compare(em, output)
end

function load_outputs(file)
    outputs = String[]
    decorators = String[]
    is_output = true
    for line in readlines(file)
        if line[1] == '+'
            push!(outputs,"")
            is_output = true
            continue
        elseif line[1] == '-'
            push!(decorators,"")
            is_output = false
            continue
        elseif line[1] == '|'
            array = is_output ? outputs : decorators
            array[end] = string(array[end],isempty(array[end])?"":"\n",line[2:end-1])
        else
            error("Unrecognized first character \"$(line[1])\"")
        end
    end
    outputs, decorators
end

function do_compare(em, i, s, output, decorator)
    TerminalUI.redraw(s, i.w)
    isa(output, String) && (output = output.data)
    isa(decorator, String) && (decorator = decorator.data)

    outbuf = IOBuffer()
    TerminalUI.render(outbuf, s)
    seekstart(outbuf); VT100.parseall!(em, outbuf)
    compare(em, output, decorator)
end

function interact(w, inputs, file)
    outputs, decorators = load_outputs(joinpath(thisdir,file))

    outbuf = IOBuffer()
    i = TerminalUI.InlineDialog(w,
        TerminalUI.IOBufferTerminal(outbuf))
    s = TerminalUI.setup_screen(i)
    s.offset = (1,1)
    s.fullsize = Base.displaysize(i.tty)
    em = VT100.ScreenEmulator()
    TerminalUI.focus(i.w)

    (length(decorators) < length(outputs)) &&
        (decorators = Any[decorators; Vector{Void}(length(outputs))])
    @test do_compare(em, i, s, outputs[1], decorators[1])
    @assert length(outputs) >= length(inputs)+1
    # Run over inputs
    for (input, output, decorator) in zip(inputs, outputs[2:end], decorators[2:end])
        inbuf = IOBuffer(input)
        seekstart(inbuf)
        result = :ok
        while result != :done && !eof(inbuf)
            result = Base.LineEdit.match_input(i.w.ctx.focuss.keymap,nothing,inbuf)(s,nothing)
        end
        @test do_compare(em, i, s, output, decorator)
    end
end

# Test cellblock
let
    cells = Vector{Cell}[
      [Cell(TerminalUI.dummycell(:blue,:red); content = 'x') for i = 1:5]
    for i = 1:5];
    w = TerminalUI.CellDisplay(cells);
    em = TerminalUI.debug(w, false)
    @test test_against(em,"outputs/cellblock.out")
end

# Test simple form
let
    myform = form(Pair[
         :name => ask("Enter your name:", AbstractString),
         :address => ask("Enter your address:", AbstractString),
         :age => ask("Enter your age:", Int)
       ],style = :incremental);
    interact(myform.chain, ["The"," Wizard\n","Oz\n","10"], "outputs/form.multiout")
end

# Test simple list
let
    fruit = ["Apple","Pear","Kiwi"]
    w = ListWidget(fruit)
    interact(w, ["\e[B","\e[B"], "outputs/list.multiout")
end


text:
using Test

left = flip_types([
  Option(1),
  Option(:b)])
right = Option([1, :b])

@test left == right


text:
using Test, Base.CoreLogging
import Base.CoreLogging: Info

using GPUCompiler

using LLVM, LLVM.Interop

using Pkg

include("util.jl")

@testset "GPUCompiler" begin

GPUCompiler.reset_runtime()

GPUCompiler.enable_timings()

include("native.jl")
include("ptx.jl")
if VERSION >= v"1.4"
  Pkg.add(["SPIRV_LLVM_Translator_jll", "SPIRV_Tools_jll"])
  include("spirv.jl")
  include("gcn.jl")
end

include("examples.jl")

haskey(ENV, "CI") && GPUCompiler.timings()

end


text:
# See [introduction](../docs/introduction.md) and [core interface](../docs/interfaces/core.md).

# # Types
#
# We start off by defining the abstract types that will be used for dispatch.

# ### [`LearningMethod`](#)
"""
    abstract type LearningMethod

Represents a concrete approach for solving a learning task.

See [core interface](../docs/interfaces/core.md) for more on
how to implement custom `LearningMethod`s
"""
abstract type LearningMethod end

# ### [`Context`](#) and concrete types
"""
    abstract type Context

Represents a context in which a data transformation
is made. This allows using dispatching for varying behavior,
for example, to apply augmentations only during training or
use non-destructive cropping during inference.

Available contexts are [`Training`](#), [`Validation`](#) and
[`Inference`](#).
"""
abstract type Context end

struct Training <: Context end
struct Validation <: Context end
struct Inference <: Context end

# # Core interface
#
# Next the core interface is defined:
# - [`encode`](#)
# - [`encodeinput`](#)
# - [`encodetarget`](#)
# - [`decodeŷ`](#)
# - [`decodey`](#)
# - [`shouldbatch`](#)
"""
    encodeinput(method, context, input) -> x

Encode `input` into a representation that a model for `method`
takes as input.

See also [`LearningMethod`](#), [`encode`](#), and [`encodetarget`](#).
"""
function encodeinput end

"""
    encodetarget(task, target; augment = false, inference = false) -> y

Encode `target` into a representation that a model for `task` outputs.
"""
function encodetarget end

"""
    encode(method, context, sample) -> (x, y)
    encode(method, context, (input, target)) -> (x, y)

Encode a `sample` containing both input and target.

If `sample` is a `Tuple` of (input, target), the default behavior is to
pass them to [`encodeinput`](#) and [`encodetarget`](#).

### Remarks

- *When should I implement `encode` vs. `encodeinput` and `encodetarget`?*

  In simple cases like image classification we can encode the inputs and
  targets separately and you should prefer `encodeinput` and `encodetarget`.
  The default implementation for `encode` when given an `(input, target)`-tuple
  is to delegate to `encodeinput` and `encodetarget`.

  In other cases like semantic segmentation, however, we want to apply
  stochastic augmentations to both image and segmentation mask. In that
  case you need to encode both at the same time using `encode`.

  Another situation where `encode` is needed is when `sample` is not a
  tuple of `(input, target)`, for example a `Dict` that includes additional
  information. `encode` still needs to return an `(x, y)`-tuple, though.
"""
function encode end

"""
    decodeŷ(method, context, ŷ) -> target

Decodes a model output into a target.
"""
function decodeŷ end

"""
    decodey(method, context, y) -> target

Decodes an encoded target back into a target.

Defaults to using [`decodeŷ`](#)
"""
decodey(method, context, y) = decodeŷ(method, context, y)

"""
    shouldbatch(method) = true

Whether models for `method` take in batches of inputs. Default
is `true`.
"""
shouldbatch(method) = true


text:
using Test
using Distributed
addprocs(1)

# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
@everywhere using MUMPSjInv

include("getDivGrad.jl")

Ar = getDivGrad(30, 31, 33)
n  = size(Ar, 1)
Ac  = Ar + im * spdiagm(0 => rand(n))

nrhs = 10
rhsr = randn(n, nrhs) + im * randn(n, nrhs)
rhsc = randn(n, nrhs)

println("Factorize real matrix on main instance and destroy")
fact = factorMUMPS(Ar, 1)
destroyMUMPS(fact)

println("Factorize real matrix on remote worker and destroy")
fact = remotecall_fetch(factorMUMPS, workers()[1], Ar, 1)
remotecall_fetch(destroyMUMPS, workers()[1], fact)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    GenericDensity{F<:Function} <: AbstractDensity

Constructors:

    GenericDensity(f)

Turns the density function `f` into a BAT-compatible [`AbstractDensity`](@ref).
The return type of `f(v)` must supported by [`logvalof`](@ref).

It must be safe to execute `f` in parallel on multiple threads and
processes.
"""
struct GenericDensity{F<:Function} <: AbstractDensity
    f::F
end

Base.convert(::Type{GenericDensity}, f::Function) = GenericDensity(f)
Base.convert(::Type{AbstractDensity}, f::Function) = GenericDensity(f)

Base.parent(density::GenericDensity) = density.f

function eval_logval_unchecked(density::GenericDensity, v::Any)
    logvalof(density.f(v))
end


text:
# Number Guessing Game in Julia

function play_number_guess_human()

    total_numbers = 25 #

    # Generate a random number within a certain range
    target_number = rand(1:total_numbers)
    guess = 0

    # While the number has not been guessed, keep prompting for guesses
    while guess != target_number
        print("Please guess a number between 1 and $total_numbers: ")
        guess = parse(Int64, readline())
        # Convert the string value input to a number

        # If we are within +/-2 of the target, give a hint
        if abs(target_number - guess) <= 2 && target_number != guess
            print("\nYou are getting closer!\n")
        end
    end

    print("Nice job, you got it!")
end


text:
# JIT compilation of Julia code to PTX

include(joinpath("compiler", "common.jl"))
include(joinpath("compiler", "irgen.jl"))
include(joinpath("compiler", "optim.jl"))
include(joinpath("compiler", "validation.jl"))
include(joinpath("compiler", "rtlib.jl"))
include(joinpath("compiler", "mcgen.jl"))
include(joinpath("compiler", "debug.jl"))
include(joinpath("compiler", "driver.jl"))

function __init_compiler__()
    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("--nvptx-fma-level=1")
end


text:
using Documenter, MultiscaleGraphSignalTransforms

makedocs(
    sitename="MultiscaleGraphSignalTransforms.jl",
    format = Documenter.HTML(collapselevel = 1),
    authors = "Jeff Irion, Haotian Li, Naoki Saito, Yiqun Shao",
    clean = true,
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "1D Path" => "examples/P64.md",
            "2D Lattice" => "examples/Grid7x3.md",
            "Sunflower Graph" => "examples/Sunflower.md",
        ],
        "Functions" => [
            "Recursive Graph Partitioning" => "functions/Partition.md",
            "Hierarchical Graph Laplacian Eigen Transform" => "functions/HGLET.md",
            "Generalized Haar-Walsh Transform" => "functions/GHWT.md",
            "extended Generalized Haar-Walsh Transform" => "functions/eGHWT.md",
            "Natural Graph Wavelet Dictionaries" => "functions/NGWD.md",
            "Utils" => "functions/Utils.md",
        ],
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/UCD4IDS/MultiscaleGraphSignalTransforms.jl.git"
)


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil: @require
using ParallelStencil.FiniteDifferences1D
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @init_parallel_stencil($package, Float64, 1)
        @require @is_initialized()
        nx  = 7
        A   =  @rand(nx  );
        Ax  =  @rand(nx+1);
        Axx =  @rand(nx+2);
        R   = @zeros(nx  );
        Rxx = @zeros(nx+2);
        @testset "1. compute macros" begin
            @testset "differences" begin
                @require @is_initialized()
                @parallel d!(R, Ax) = (@all(R) = @d(Ax); return)
                @parallel d2!(R, Axx) = (@all(R) = @d2(Axx); return)
                R.=0; @parallel d!(R, Ax);  @test all(R .== Ax[2:end].-Ax[1:end-1])
                R.=0; @parallel d2!(R, Axx);  @test all(R .== (Axx[3:end].-Axx[2:end-1]).-(Axx[2:end-1].-Axx[1:end-2]))
            end;
            @testset "selection" begin
                @require @is_initialized()
                @parallel all!(R, A) = (@all(R) = @all(A); return)
                @parallel inn!(R, Axx) = (@all(R) = @inn(Axx); return)
                R.=0; @parallel all!(R, A);  @test all(R .== A)
                R.=0; @parallel inn!(R, Axx);  @test all(R .== Axx[2:end-1])
            end;
            @testset "averages" begin
                @require @is_initialized()
                @parallel av!(R, Ax) = (@all(R) = @av(Ax); return)
                R.=0; @parallel av!(R, Ax);  @test all(R .== (Ax[1:end-1].+Ax[2:end]).*0.5)
            end;
            @testset "others" begin
                @require @is_initialized()
                @parallel maxloc!(R, Axx) = (@all(R) = @maxloc(Axx); return)
                R.=0; @parallel maxloc!(R, Axx);  @test all(R .== max.(max.(Axx[3:end],Axx[2:end-1]),Axx[1:end-2]))
            end;
        end;
        @testset "2. apply masks" begin
            @testset "selection" begin
                @require @is_initialized()
                @parallel inn_all!(Rxx, A) = (@inn(Rxx) = @all(A); return)
                @parallel inn_inn!(Rxx, Axx) = (@inn(Rxx) = @inn(Axx); return)
                Rxx.=0; @parallel inn_all!(Rxx, A);  @test all(Rxx[2:end-1] .== A)
                Rxx[2:end-1].=0; @test all(Rxx .== 0)  # Test that boundary values remained zero.
                Rxx.=0; @parallel inn_inn!(Rxx, Axx);  @test all(Rxx[2:end-1] .== Axx[2:end-1])
                Rxx[2:end-1].=0; @test all(Rxx .== 0)  # Test that boundary values remained zero.
            end;
            @testset "differences" begin
                @require @is_initialized()
                @parallel inn_d!(Rxx, Ax) = (@inn(Rxx) = @d(Ax); return)
                @parallel inn_d2!(Rxx, Axx) = (@inn(Rxx) = @d2(Axx); return)
                Rxx.=0; @parallel inn_d!(Rxx, Ax);  @test all(Rxx[2:end-1] .== Ax[2:end].-Ax[1:end-1])
                Rxx[2:end-1].=0; @test all(Rxx .== 0)  # Test that boundary values remained zero.
                Rxx.=0; @parallel inn_d2!(Rxx, Axx);  @test all(Rxx[2:end-1] .== (Axx[3:end].-Axx[2:end-1]).-(Axx[2:end-1].-Axx[1:end-2]))
                Rxx[2:end-1].=0; @test all(Rxx .== 0)  # Test that boundary values remained zero.
            end;
        end;
        @reset_parallel_stencil()
    end;
)) end == nothing || true;


text:
using NearestNeighbors, StaticArrays, LinearAlgebra
using Distances: Metric, Cityblock, Euclidean

export broomhead_king
export numericallyapunov
export Cityblock, Euclidean

#####################################################################################
#                    Numerical Lyapunov (from reconstruction)                       #
#####################################################################################
# Everything in this section is based on Ulrich Parlitz [1]

"""
    numericallyapunov(R::Dataset, ks;  refstates, w, distance, ntype)
Return `E = [E(k) for k ∈ ks]`, where `E(k)` is the average logarithmic distance
between states of a [`neighborhood`](@ref)
that are evolved in time for `k` steps (`k` must be integer).

## Keyword Arguments

* `refstates = 1:(length(R) - ks[end])` : Vector of indices
  that notes which
  states of the reconstruction should be used as "reference states", which means
  that the algorithm is applied for all state indices contained in `refstates`.
* `w::Int = 1` : The Theiler window, which determines
  whether points are separated enough in time to be considered separate trajectories
  (see [1] and [`neighborhood`](@ref)).
* `ntype::AbstractNeighborhood = FixedMassNeighborhood(1)` : The method to
  be used when evaluating the neighborhood of each reference state. See
  [`AbstractNeighborhood`](@ref) or [`neighborhood`](@ref) for more info.
* `distance::Metric = Cityblock()` : The distance function used in the
  logarithmic distance of nearby states. The allowed distances are `Cityblock()`
  and `Euclidean()`. See below for more info.

## Description
If the dataset/reconstruction
exhibits exponential divergence of nearby states, then it should clearly hold
```math
E(k) \\approx \\lambda\\Delta t k + E(0)
```
for a *well defined region* in the `k` axis, where ``\\lambda`` is
the approximated
maximum Lyapunov exponent. ``\\Delta t`` is the time between samples in the
original timeseries.
You can use [`linear_region`](@ref) with arguments `(ks .* Δt, E)` to
identify the slope
(= ``\\lambda``)
immediatelly, assuming you
have choosen sufficiently good `ks` such that the linear scaling region is bigger
than the saturated region.

The algorithm used in this function is due to Parlitz [1], which itself
expands upon Kantz [2]. In sort, for
each reference state a neighborhood is evaluated. Then, for each point in this
neighborhood, the logarithmic distance between reference state and neighborhood
state is
calculated as the "time" index `k` increases. The average of the above over
all neighborhood states over all reference states is the returned result.

If the `Metric` is `Euclidean()` then use the Euclidean distance of the
full `D`-dimensional points (distance ``d_E`` in ref. [1]).
If however the `Metric` is `Cityblock()`, calculate
the absolute distance of *only the first elements* of the `m+k` and `n+k` points
of the reconstruction `R` (distance ``d_F`` in ref. [1]).

## References

[1] : Skokos, C. H. *et al.*, *Chaos Detection and Predictability* - Chapter 1
(section 1.3.2), Lecture Notes in Physics **915**, Springer (2016)

[2] : Kantz, H., Phys. Lett. A **185**, pp 77–87 (1994)
"""
function numericallyapunov(R::AbstractDataset{D, T}, ks;
                           refstates = 1:(length(R) - ks[end]),
                           w = 1,
                           distance = Cityblock(),
                           ntype = FixedMassNeighborhood(1)) where {D, T}
    Ek = numericallyapunov(R, ks, refstates, w, distance, ntype)
end

function numericallyapunov(R::AbstractDataset{D, T},
                           ks::AbstractVector{Int},
                           ℜ::AbstractVector{Int},
                           w::Int,
                           distance::Metric,
                           ntype::AbstractNeighborhood) where {D, T}

    # ℜ = \Re<tab> = set of indices that have the points that one finds neighbors.
    # n belongs in ℜ and R[n] is the "reference state".
    # Thus, ℜ contains all the reference states the algorithm will iterate over.
    # ℜ is not estimated. it is given by the user. Most common choice:
    # ℜ = 1:(length(R) - ks[end])

    # ⩅(n) = \Cup<tab> = neighborhood of reference state n
    # which is evaluated for each n and for the given neighborhood ntype

    # Initialize:
    timethres = length(R) - ks[end]
    if maximum(ℜ) > timethres
        erstr = "Maximum index of reference states is > length(R) - ks[end] "
        erstr*= "and the algorithm cannot be performed on it. You have to choose "
        erstr*= "reference state indices of at most up to length(R) - ks[end]."
        throw(ArgumentError(erstr))
    end
    E = zeros(T, length(ks))
    E_n = copy(E); E_m = copy(E)
    data = R.data #tree data
    tree = KDTree(data, Euclidean()) # this creates a copy of `data`
    skippedm = 0; skippedn = 0

    for n in ℜ
        # The ⋓(n) can be evaluated on the spot instead of being pre-calculated
        # for all reference states. (it would take too much memory)
        # Since ⋓[n] doesn't depend on `k` one can then interchange the loops:
        # Instead of k being the outermost loop, it becomes the innermost loop!
        point = data[n]
        ⋓ = neighborhood(point, tree, ntype, n, w)
        for m in ⋓
            # If `m` is nearer to the end of the timeseries than k allows
            # is it completely skipped (and length(⋓) reduced).
            # If m is closer to n than the Theiler window allows, also skip.
            # The Theiler window defaults to τ
            if m > timethres || abs(m - n) <= w
                skippedm += 1
                continue
            end
            for (j, k) in enumerate(ks) #ks should be small (of order 10 to 100 MAX)
                E_m[j] = log(delay_distance(distance, R, m, n, k))
            end
            E_n .+= E_m # no need to reset E_m
        end
        if skippedm >= length(⋓)
            skippedn += 1
            skippedm = 0
            continue # be sure to continue if no valid point!
        end
        E .+= E_n ./ (length(⋓) - skippedm)
        skippedm = 0
        fill!(E_n, zero(T)) #reset distances for n reference state
    end

    if skippedn >= length(ℜ)
        ers = "skippedn ≥ length(R)\n"
        ers*= "Could happen because all the neighbors fall within the Theiler "
        ers*= "window. Fix: increase neighborhood size."
        error(ers)
    end
    E ./= length(ℜ) - skippedn
end

@inline function delay_distance(di::Cityblock, R, m, n, k)
    @inbounds abs(R[m+k][1] - R[n+k][1])
end

@inline function delay_distance(di::Euclidean, R, m, n, k)
    @inbounds norm(R[m+k] - R[n+k])
end

#####################################################################################
#                                  Broomhead-King                                   #
#####################################################################################
"""
    broomhead_king(s::AbstractVector, d::Int) -> U, S, Vtr
Return the Broomhead-King coordinates of a timeseries `s`
by performing `svd` on the so-called trajectory matrix with dimension `d`.

## Description
Broomhead and King coordinates is an approach proposed in [1] that applies the
Karhunen–Loève theorem to delay coordinates embedding with smallest possible delay.

The function performs singular value decomposition
on the `d`-dimensional trajectory matrix ``X`` of ``s``,
```math
X = \\frac{1}{\\sqrt{N}}\\left(
\\begin{array}{cccc}
x_1 & x_2 & \\ldots & x_d \\\\
x_2 & x_3 & \\ldots & x_{d+1}\\\\
\\vdots & \\vdots & \\vdots & \\vdots \\\\
x_{N-d+1} & x_{N-d+2} &\\ldots & x_N
\\end{array}
\\right) = U\\cdot S \\cdot V^{tr}.
```
where ``x := s - \\bar{s}``.
The columns of ``U`` can then be used as a new coordinate system, and by
considering the values of the singular values ``S`` you can decide how many
columns of ``U`` are "important". See the documentation page for example application.

## References
[1] :  D. S. Broomhead, R. Jones and G. P. King, J. Phys. A **20**, 9, pp L563 (1987)
"""
function broomhead_king(x::AbstractArray, d::Int)
    X = trajectory_matrix(x, d)
    F = LinearAlgebra.svd(X)
    return F.U, F.S, F.Vt
end

function trajectory_matrix(x::AbstractArray, d::Int)
    xdash = mean(x)
    N = length(x); sqN = √N
    X = zeros(N-d+1, d)
    for j in 1:d
        for i in 0:(N-d)
            @inbounds X[i+1, j] = (x[j+i] - xdash)/sqN
        end
    end
    return X
end


text:
import YAML

"""
    get_turb_loc_YAML(file_name)

read in turbine locations and related problem file names from .yaml

# Arguments
- `file_name::String`: path/and/name/of/location/file.yaml
"""
function get_turb_loc_YAML(file_name; returnaep=false)
    ### Retrieve turbine locations and auxiliary file names from <.yaml> file.
    ### Auxiliary (reference) files supply wind rose and turbine attributes.

    # Read in the .yaml file
    f = YAML.load(open(file_name))
    defs = f["definitions"]

    # Rip the (x,y) coordinates (Convert from <list> to <ndarray>)
    turb_coords = defs["position"]["items"]
    # println(turb_coords)
    nturbs = length(turb_coords)
    turbine_x = zeros(nturbs)
    turbine_y = zeros(nturbs)

    for i in 1:nturbs
        turbine_x[i] = turb_coords[i][1]
        turbine_y[i] = turb_coords[i][2]
    end

    # Rip the expected AEP, used for comparison
    AEP = defs["plant_energy"]["properties"]["annual_energy_production"]["default"]

    # Read the auxiliary filenames for the windrose and the turbine attributes (first one)
    fname_turb = string.(values(defs["wind_plant"]["properties"]["turbine"]["items"][1]))[1]
    fname_wr = string.(values(defs["plant_energy"]["properties"]["wind_resource"]["properties"]["items"][1]))[1]

    # Return turbine (x,y) locations, and the filenames for the others .yamls
    if returnaep
        return turbine_x, turbine_y, fname_turb, fname_wr, AEP
    else
        return turbine_x, turbine_y, fname_turb, fname_wr
    end
end

"""
    get_turb_atrbt_YAML(file_name)

read in turbine attributes from .yaml

# Arguments
- `file_name::String`: path/to/attribute/file.yaml
"""
function get_turb_atrbt_YAML(file_name)
    ###Retreive turbine attributes from the <.yaml> file###

    # Read in the .yaml file
    f = YAML.load(open(file_name))
    defs = f["definitions"]
    ops = defs["operating_mode"]
    turb = defs["wind_turbine"]
    rotor = defs["rotor"]
    hub = defs["hub"]

    # Rip the turbine attributes
    # (Convert from <list> to <float>)
    turb_ci = float(ops["cut_in_wind_speed"]["default"])
    turb_co = float(ops["cut_out_wind_speed"]["default"])
    rated_ws = float(ops["rated_wind_speed"]["default"])
    rated_pwr = float(turb["rated_power"]["maximum"])
    turb_diam = float(rotor["diameter"]["default"])
    turb_height = float(hub["height"]["default"])

    return turb_ci, turb_co, rated_ws, rated_pwr, turb_diam, turb_height
end

"""
    get_wind_rose_YAML(file_name)

read in wind resource information from .yaml

# Arguments
- `file_name::String`: path/to/wind/resource/file.yaml
"""
function get_wind_rose_YAML(file_name)
    ### Retrieve wind rose data (bins, freqs, speeds) from <.yaml> file.

    # Read in the .yaml file
    f = YAML.load(open(file_name))
    props = f["definitions"]["wind_inflow"]["properties"]

    # Rip wind directional bins, their frequency, and the windspeed parameters for each bin
    # (Convert from <list> to <ndarray>)
    wind_dir = props["direction"]["bins"]
    wind_dir_freq = props["direction"]["frequency"]
    ndirs = length(wind_dir)
    # (Convert from <list> to <float>)
    wind_speeds = props["speed"]["bins"]
    wind_speed_probs = props["speed"]["frequency"]
    # Get default number of windspeed bins per direction
    nspeeds = length(wind_speeds)
    min_speed = props["speed"]["minimum"]
    max_speed = props["speed"]["maximum"]

    # get ambient ti
    ti = props["turbulence_intensity"]["default"]

    # convert to flow farm standard format
    freq = zeros(ndirs*nspeeds)
    speed = zeros(ndirs*nspeeds)
    dir = zeros(ndirs*nspeeds)

    # calculate frequency for each direction/speed combination
    for i in 1:ndirs
        for j in 1:nspeeds
            freq[(i-1)*nspeeds+j] = wind_dir_freq[i]*wind_speed_probs[i][j]
            speed[(i-1)*nspeeds+j] = wind_speeds[j]
            dir[(i-1)*nspeeds+j] = wind_dir[i]
        end
    end

    return dir, speed, freq, ti
end

"""
    write_turb_loc_YAML(file_name, data)

write turbine locations and related information to .yaml

# Arguments
- `file_name::String`: path/and/name/of/location/file.yaml
"""
function write_turb_loc_YAML(filename, turbinex, turbiney; title="", titledescription="",
    turbinefile="", locunits="m", wakemodelused="", windresourcefile="", aeptotal=[],
    aepdirs=[], aepunits="MWh", baseyaml=string(@__DIR__, "/default.yaml"))

    ### Retrieve turbine locations and auxiliary file names from <.yaml> file.
    ### Auxiliary (reference) files supply wind rose and turbine attributes.

    # Read in the default .yaml file
    base = YAML.load(open(baseyaml))

    # get number of turbines
    nturbines = length(turbinex)

    # save the title and description to the yaml database
    base["title"] = title
    base["description"] = titledescription

    # save the title and description to the yaml database
    base["definitions"]["plant_energy"]["properties"]["wake_model"]["items"][1]["\$ref"] = wakemodelused

    # save positions in yaml database
    turb_coords = fill(zeros(2), nturbines)
    for i in 1:nturbines
        turb_coords[i] = [turbinex[i], turbiney[i]]
    end

    base["definitions"]["position"]["items"] = turb_coords

    # save the AEP in yaml database
    base["definitions"]["plant_energy"]["properties"]["annual_energy_production"]["default"] = aeptotal

    # save the directional AEPs in the yaml database
    base["definitions"]["plant_energy"]["properties"]["annual_energy_production"]["binned"] = aepdirs

    # save the directional AEP units in the yaml database
    base["definitions"]["plant_energy"]["properties"]["annual_energy_production"]["units"] = aepunits

    # save the auxiliary filenames for the windrose and the turbine attributes
    base["definitions"]["wind_plant"]["properties"]["turbine"]["items"][1]["\$ref"] = turbinefile
    base["definitions"]["plant_energy"]["properties"]["wind_resource"]["properties"]["items"][1]["\$ref"] = windresourcefile

    # write result to .yaml
    YAML.write_file(filename, base)
end

function get_reduced_wind_rose_YAML(file_name)
    ### Retrieve wind rose data (bins, freqs, speeds) from <.yaml> file.
    # Only the average wind speeds for each direction are reported.

    # Read in the .yaml file
    f = YAML.load(open(file_name))
    props = f["definitions"]["wind_inflow"]["properties"]

    # Rip wind directional bins, their frequency, and the windspeed parameters for each bin
    # (Convert from <list> to <ndarray>)
    wind_dir = props["direction"]["bins"]
    wind_dir_freq = props["direction"]["frequency"]
    ndirs = length(wind_dir)
    # (Convert from <list> to <float>)
    wind_speeds = props["speed"]["bins"]
    wind_speed_probs = props["speed"]["frequency"]
    # Get default number of windspeed bins per direction
    nspeeds = length(wind_speeds)
    min_speed = props["speed"]["minimum"]
    max_speed = props["speed"]["maximum"]

    # get ambient ti
    ti = props["turbulence_intensity"]["default"]

    # convert to flow farm standard format
    freq = zeros(ndirs)
    speed = zeros(ndirs)
    dir = zeros(ndirs)

    # calculate average speed for each direction
    for i in 1:ndirs
        freq[i] = wind_dir_freq[i]
        speed[i] = sum(wind_speeds.*wind_speed_probs[i])
        dir[i] = wind_dir[i]
    end

    return dir, speed, freq, ti
end

"""
    getNextFileName(directory, file_name, file_type)

Checks if a file of the given directory and name exists. If not, increments to
the next index so as not to overwrite previously written files. If it reaches
the max number of overwrites, it will default to <directory/file_name.file_type>
To default to this, set <max_check=0> in function call.

# Arguments
- `directory::String`: path/to/write/file/at/
- `file_name::String`: Whatever the name of the file desired
- `file_type::String`: ex "yaml", "txt", "csv", etc...
- `max_check::Int`: the maximum number of files to check
"""
function getNextFileName(directory, file_name, file_type, max_check=100)
    full_file_name = string(directory, file_name, ".", file_type)
    for i in 1:100
        if !isfile(string(directory, file_name, "-(", i, ").", file_type))
            full_file_name = string(directory, file_name, "-(", i, ").", file_type)
            break
        end
    end
    return full_file_name
end


text:
using Flux
using DiffEqSensitivity
using DifferentialEquations
using LinearAlgebra
using Statistics
using Distributions
using Plots
using Random

# set random seed for reproducible results
Random.seed!(1)

## ----- setup true solution ------------------------------------

# define function for RHS of lorenz system
function lorenz!(du, u, p, t)
    x,y,z = u
    σ,ρ,β = p

    du[1] = σ*(y - x)
    du[2] = x*(ρ - z) - y
    du[3] = x*y - β*z
end

θ = [10.0, 28.0, (8/3)]  # true parameter values
dt = 0.01  # time step for simulation (this may be unnecessary)
tend = 10.0  # max time for total simulation
ts_sim = (0.0, tend)
u0True = [1.0, 1.0, 1.0]  # true initial condition

# integrate to obtain the true solution
prob_true = ODEProblem(lorenz!, u0True, ts_sim, θ, saveat=dt)
sol_true = solve(prob_true)

## ----- setup observation data ------------------------------------

# set up observation uncertainties
σ_m = 0.15  # i.e. observation std
# set up observation error covariance matrix
R = σ_m^2 .* I(3)

# set up details for measurement points
dt_m = 0.2
tm_end = 2.0  # final observation time
ts_m = dt_m:dt_m:tm_end

# define observation operator
function h(u)
    w = u
    return w
end

# populate observation vector
w = hcat([h(sol_true(tᵢ))+rand(MvNormal(zeros(3), R)) for tᵢ ∈ ts_m]...)

## ----- perform the assimilation ------------------------------------

# initial condition guess
u0a = [2.0, 3.0, 4.0]

# set up problem for (bad) forecast
prob = ODEProblem(lorenz!, u0a, (0.0,tend), θ)
sol_forecast = solve(prob)

function timeindex(t, Δt)
    return Int(round(t/Δt)) + 1
end

# for some reason I need to collect the matrices first
Rinv = collect(inv(R))
u0b = copy(u0a)
σ_b = 0.1  # model error std
B = σ_b^2 .* I(3)  # initial model error covariance
Binv = collect(inv(B))

function loss()
    res = Array(solve(ODEProblem(lorenz!, u0a, (0.0,tm_end), θ), saveat=dt_m:dt_m:tm_end))
    # ℓ = 0.5*(u0a-u0b)'*Binv*(u0a-u0b) + 0.5*sum((w[:,i]-h(res[:,i]))'*Rinv*(w[:,i]-h(res[:,i])) for i ∈ 1:size(res,2))
    ℓ = 0.5*sum((w[:,i]-h(res[:,i]))'*Rinv*(w[:,i]-h(res[:,i])) for i ∈ 1:size(res,2))
    return ℓ
end

u0a = [2.0, 3.0, 4.0]
loss()

# let's do the optimization
data = Iterators.repeated((), 500)
opt = ADAM(0.1)

cb = function ()
    display(loss())
end

Flux.train!(loss, Flux.params(u0a), data, opt, cb=cb)

u0a

prob = ODEProblem(lorenz!, u0a, (0.0, tend), θ)
sol_analysis = solve(prob)

var_names = ["x(t)", "y(t)", "z(t)"]
plots = []
for i ∈ 1:3
    p = plot(sol_true, vars=(0,i), color=:black, alpha=0.75, linewidth=2, label="true system")
    plot!(p, ts_m, w[i, :], seriestype=:scatter, color="light green", label="observation" )
    plot!(p, sol_forecast, vars=(0, i), linestyle=:dashdot, color="purple", label="forecast")
    plot!(p, sol_analysis, vars=(0,i), color="royal blue", alpha=1, linewidth = 1, label="analysis")
    xlabel!(p, "t")
    ylabel!(p, var_names[i])
    push!(plots, p)
end

plot(plots..., layout=(3,1))

savefig("4dvar_using_flux.svg")
savefig("4dvar_using_flux.png")


text:
module MKL

    using Libdl
    using LinearAlgebra: BlasInt

    include(joinpath("..", "deps", "deps.jl"))

    @enum Threading begin
        THREADING_INTEL
        THREADING_SEQUENTIAL
        THREADING_PGI
        THREADING_GNU
        THREADING_TBB
    end
    @enum Interface begin
        INTERFACE_LP64
        INTERFACE_ILP64
        INTERFACE_GNU
    end

    function set_threading_layer(layer::Threading = THREADING_INTEL)
        err = ccall((:MKL_Set_Threading_Layer, libmkl_rt), Cint,
            (Cint,), layer)
        if err == -1
            throw(ErrorException("return value was -1"))
        end
        return nothing
    end

    function set_interface_layer(interface = BlasInt == Int64 ? INTERFACE_ILP64 : INTERFACE_LP64)
        err = ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint,
            (Cint,), interface)
        if err == -1
            throw(ErrorException("return value was -1"))
        end
        return nothing
    end

    function __init__()

        if Sys.iswindows()
            # On Windows, we have to open the threading library before calling libmkl_rt
            dlopen(libmkl_intel_thread)
        end

        set_threading_layer()
        set_interface_layer()
    end

    mklnorm(x::Vector{Float64}) = ccall((:dnrm2_, libmkl_rt), Float64,
        (Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
        length(x), x, 1)

    # Include installation routines as well
    include("install.jl")
end