text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 17fc4346-a462-11eb-0484-7b4029b01f4d
using Images

# ╔═╡ 3ca7a180-e6ad-4734-b19f-4480de52d13c
md"
# Dataset: Water Clusters as Graphs
"

# ╔═╡ 30339b9c-95a9-4a0e-b7d6-abd65f958e3d
md"
The database we will be working with primarily comes from my own research. In the past, we have generated a very large database of molecular structures called **water clusters**. Water clusters are gas-phase blobs of water molecules which are held together by **hydrogen bonds**. They can consist of essentially any number of water molecules. The actual database consists of the cartesian coordinates of each atom in 3D space and can be accessed at [sites.uw.edu/wdbase](https://sites.uw.edu/wdbase/database-of-water-clusters/). Because of the complexity of analyzing this data, I will be working with a simplified version of this database which treats each water cluster as a graph.
"

# ╔═╡ 5bc0d0b1-0078-4f84-baa0-d0267f8bd7b6
load("./assets/Figure_1_a_W3_25_image_combined.png")

# ╔═╡ 2761102e-78c2-4a4a-a7e2-c972c1bd1f4a
md"
The above image shows the **lowest energy** structures of $\mathrm{(H_2O)_n}$ where $\mathrm{n}$ is the number of water molecules present in the cluster. As you can see, the hydrogen bonds connecting each water molecule can be thought of as defining a graph. In these graphs, each node represents a water molecule and each edge represents a hydorgen bond. We can treat these as either simple graphs (no edge direction) or directed graphs in which each edge is directed along the hydrogen bonds, pointing from a hydrogen atom to an oxygen atom.
"

# ╔═╡ 420184cc-1258-4f98-880e-40dc47914011
load("./assets/cluster_to_graph_illustration_annotated.png")

# ╔═╡ 6b3f36e0-efa5-4010-8baa-758b4a72e75d
md"
The above image illustrates this transformation. On the left is a visualization of each atom in a $\mathrm{(H_2O)_{10}}$ cluster. In the middle is a depiction of a directed graph. The edge directions go from the hydrogen atom (white atoms in the left picture) to an oxygen atom (red atoms). Finally, the right is a simple graph representing the structure. This undirected version will sometimes be referred to as a \"family\" because it defines the connectivity of molecules, but there are many different ways the molecules can be conected by hydrogen bonds, and each of these clusters will have a different energy.
"

# ╔═╡ 78efba80-47b8-4c9c-9435-a9034e7d4053
md"
With that context out of the way, let's describe the actual data I've extracted from these graphs. First, I have enumerated the number of primitive cycles in the simple graphs from size three up to size ten. As an example, the above simple graph contains two 5-cycles (the pentagons) and five 4-cycles (the squares which connect the pentamers). However, we do not count any 8-cycles which would be formed by two neighboring 4-cycles. In this sense, the 4- and 5-cycles are the primitives of this graph.

Additionally, we have information from the directed graphs which describes the connectivity of each molecule in a structure. Namely, molecule 1 in the directed graph above is conventionally described as an AAD molecule. This means it accepts two hydrogen bonds and donates one. Molecule 2 would be called ADD since it accepts one hydrogen bond and donates two.

Finally, we have the total energy of each molecule. The energy is a measure of the stability of each water cluster. Roughly speaking, the lower the energy, the more physically relevant the structure is, as it is much more likely to be present in experiments or the atmosphere, etc.
"

# ╔═╡ 21a82ca1-6b0b-489a-82e8-9a19804829e7
md"
## Questions of Interest:
"

# ╔═╡ 8e6b5906-b80a-45d7-aa8b-b73b670f65da
md"
Some questions I would like to be able to answer are the following:
1. How does the distribution of rings change with cluster size? We should expect 6-membered rings to outpace 3-, 4-, and 5-membered rings as clusters get larger because 6-membered rings are dominant in solid ice.

2. Can we confirm the conventional wisdom that more hydrogen bonds will lead to more stable clusters?

3. Are there any interesting trends when relating different hydrogen bonding configurations (e.g. AAD vs. ADD) to the total energy of a cluster?
"

# ╔═╡ 7c8a60f4-896a-4fdc-9971-1fa48d9df27a
md"
### How are rings distributed with cluster size?
"

# ╔═╡ 9e29427e-4967-4d27-a885-0ac6d0319da8
md"
Let's begin by looking at the distribution of these ring cycles with the number of molecules in each cluster.
"

# ╔═╡ 948c2330-ea9b-4350-a276-ebe940340412
load("assets/rings_vs_cluster_size_lowest_family.png")

# ╔═╡ 47f26447-a892-49cc-926f-ca1db8ee64b9
md"
I have plotted the number of each size cycle against the size of the water cluster. The points represent an average across all graphs of that size, and the shaded area is the standard deviation of each distribution.

We can see that as the clusters get larger, the larger rings becomes much more common. This is to be expected just from the combinatorial nature of the number of paths of a given size.

Typically, chemists think of the molecular structures containing rings in terms of 3-, 4-, 5-, and 6-membered rings. For this reason, I have plotted these distributions in the inset. We can see that at smaller cluster sizes, the smaller rings are more prevalent. However, at 21 molecules, 6-cycles become more common and seem to stay more prevalent than 3-, 4-, and 5-cycles. This is the behavior we would hope to see because bulk ice consists of a lattice of 6-cycles, so we should approach that limit, which seems to be the case.
"

# ╔═╡ d6bcbd4e-d4e6-4e7f-a7eb-41fcb7cac2d7
load("assets/ring_size_correlation.png")

# ╔═╡ 7d873195-eeed-4582-a344-9e7c90eae459
md"
Let's see if these cycle sizes have any strong correlations amongst each other. Above is the pair correlations of the counts of each cycle size from 3 to 10. Roughly speaking, it seems that clusters which are more similar in size are more likely to appear at similar rates. Perhaps this is unsurprising.
"

# ╔═╡ 54f8cd88-f42e-42f2-a62f-8149926a21de
md"
### How are different water bonding arrangements distributed with respect to cluster size?
"

# ╔═╡ 1fe20c54-4c54-4f66-893b-89ff3aef08b4
load("assets/population_of_ring_types_with_cluster_size.png")

# ╔═╡ a5a29ab5-fb59-44d7-9db5-356010f00b23
md"
Let's examine the molecule-level data we have for each graph. As a reminder, AAADD means that a molecule would have three acceptors and two donor molecules. In terms of graphs, this would mean a particular node in the graph has in-degree 3 and out-degree two.

Since there are many more nodes in larger clusters, the counts of each occurence are heavily skewed towards towards the larger water clusters. To a certain extent, this is good because these clusters are better models of liquid water and/or ice. On the other hand, we still care about the smaller clusters, so I have plotted the data on a log-plot on the left and on linear axes on the right. The color bar indicates the number of water molecules present in the cluster. It is commonly said that water is tetrahedrally coordinated. This would correspond to the AADD configuration, which is indeed the most common type of node, and it seems to be getting more common as clusters get larger. Nonetheless, AAD and ADD configurations are still fairly common

Somewhat surprisingly, there is a non-negligible number of AAADD molecules. This is surprising because people usually think of water as only being capable of accepting two hydrogen bonds and donating two hydrogen bonds. Further work needs to be done to understand the importance of these structures.
"

# ╔═╡ b35efb02-713a-4795-a249-84cb4d506abf
md"
### How are hydrogen bonds related to cluster stability?
"

# ╔═╡ 153c2888-1cdf-41d3-bafb-8a1bbeaab03a
load("assets/strength_of_hbond_with_cluster_size.png")

# ╔═╡ 6fe153a6-5262-4baa-b0f1-19cf140c3f0a
md"
When we plot the energy of each cluster against the number of hydrogen bonds (sum of degrees of each node), we get the above plot. First of all, the banding occures because there are discrete jumps in energy when the cluster size increases by one. Since hydrogen bonds are the stabilizing interaction between water molecules, conventional wisdom would say that the more hydrogen bonds are present, the more stable the cluster will be. The above plot seems to show that at lower cluster sizes this is the case, but as the cluster gets bigger, there is a wider range of number of hdyrogen bonds which correspond to low energies. This can be seen by the growing width of each band.

It may be illuminating to normalize these numbers in some way. Let's try that.
"

# ╔═╡ 34627d00-e9fc-4357-b626-4c5c0b1ed6b4
load("assets/strength_of_hbond_with_cluster_size_normalized.png")

# ╔═╡ 7c5a19cd-9952-4ecc-adef-676031112e99
md"
In the above two visualizations, I have normalized the x-axis by the number of molecules in the cluster (this is a common normalization in the literature). The plot on the left normalizes the y-axis by the energy, while the plot on the right normalizes by the number of hydrogen bonds.

These plots are both very illuminating in slightly different ways. The figure on the left shows that the binding energy per molecule increases until there are about three hydrogen bonds per molecule. For reference, liquid water has an average of about 3.6 hydrogen bonds per molecule, and ice has exactly 4 per molecule. Once clusters are large enough to have this many hydrogen bonds per molecule, the binding energy seems to level off. This can be thought of as the clusters reaching some bulk limit.

One might be led to believe that the hydrogen bonds are getting stronger with cluster size by looking at the plot on the left, but in fact, the plot on the right indicates that if anything, there is a slight decrease in the average hdyrogen bond strength with cluster size. That is, just as we have seen above, the number of hydrogen bonds in a cluster is not necessarily a good preictor of the stability of a cluster because accomodating more hydrogen bonds might require adopting sub-optimal configurations and hence result in weaker hydrogen bonds.
"

# ╔═╡ b920bfbb-82ad-4f57-8e1c-f0001e7d77a0
md"
Having looked at the correlation between hydrogen bonds and energy, let's see if there's any similar correlation between the type of molecules in a cluster and that cluster's energy.
"

# ╔═╡ 0eda0aaf-63b7-4540-af84-667ca53d57b0
md"
### How are bonding configurations of water related to cluster stability?
"

# ╔═╡ e7895abf-ad9b-4153-9927-2c5215b52562
load("assets/bonding_configuration_correlation.png")

# ╔═╡ b3527c13-bda8-4b8d-a48d-8e20631eec52
md"
Let's begin by explorig the correlations between different bonding types. Above, I have plotted the five most prominent ones correlated against one another. There are a few intersting trends here. One is that AAD and ADD configurations are strongly correlated to one another. This makes sense I think, as these configurations are most likely to be found on the surface of a cluster, and there are many surface molecules in these size clusters generally.

Another trend is that for larger clusters (N > 20), the presence of AD molecules greatly decreases the number of AAD, ADD, and AADD molecules. Finally, a large number of AAADD molecules tends to be anti-correlated to AADD molecules in particular.

Now that we have a sense of the data, let's look at how different bonding arrangements contribute to molecule stability. In particular, let's focus on clusters with 20-30 molecules since this data skews towwards large clusters quite dramatically.
"

# ╔═╡ c9b0f626-65c3-4912-9b66-e76f7decc0f6
load("assets/rings_vs_energy_20_to_30.png")

# ╔═╡ aa2ee694-5f1e-47ac-93e7-f3f99023103f
md"
For visual clarity, we do this analysis for clusters of 20 to 30 molecules. Once again, we get distinct bands corresponding to cluster size. The most obvious trend is that having zero AD molecules is definitely the most favorable situation for finding a low-energy cluster.

Note that the dashed line is the lowest energy structure at each cluster size.

For ADD and AAD molecules, it seems to be that being somewhere in the middle of the distribution is certainly the most favorable. This makes sense given the trade-off between number of hydrogen bonds and overall hydrogen bond strength which we just saw.

Perhaps the most interesting trend is that the number of AADD molecules doesn't seem to matter nearly as much as having an even number of AADD molecules. This is much easier to see if we plot this by itself.
"

# ╔═╡ ef7eb7e2-492e-4203-9106-90342735e47c
load("assets/rings_vs_energy_20_to_30_AADD_only.png")

# ╔═╡ 5e708228-6075-47f2-9153-991e8dcfb3f3
md"
While it's not the most pronounced trend ever, it's now clearer that having 4, 6, 8, or 10 AADD molecules is likely to be more stable than having 2, 3, 5, 7, 9, or 11. Two is likely just too few hydrogen bonds overall. The odd-numbered structures must bias the structures in some way as to make the connectivity less favorable. This is an intriguing but subtle trend which will have to be investigated further in the future.
"

# ╔═╡ e065e882-a5cb-4c5f-bead-434f839d909a
md"
## Conclusions:
"

# ╔═╡ 48d3550d-c965-47c3-995e-d8c6fc69606e
md"
To answer the questions we set out with:
1. We were able to show that 6-cycle rings do indeed outpace 3-, 4-, and 5-cycles as clusters get larger. Even larger rings seem to get more important. Further investigation of these rings is required though because these are not rings which chemists usually think of as being important.

2. We did see the trend that more hydrogen bonds generally correlates with a more stable structure, but we found that as clusters get larger, the number of hydrogen bonds can vary quite a bit while still giving a similar total energy. In the future, identifying the soecific arrangement which are most stable would be quite interesting.

3. Finally, we were able to tease out interesting trends associated with the number of different hydrogen bonding configurations and the stability of the cluster. Namely, AADD configurations seem to be most stabilizing when there are an even number of configurations. This was unexpected and merits further investigation. It was also clear that having AD molecules is generally quite bad for having the lowest energy at a given cluster size.
"

# ╔═╡ 7df79f4b-14e0-483a-afa9-30e66cc4d537
md"
Some comments on data quality:

Fortunately, because this data comes from my own research, I was able to guarantee that there were not any missing values in the dataset. On the other hand, the definition of rings within a cluster is somewhat arbitrary, so the abundance of large ring sizes may be an artifact of the way the analysis is carried out, and hence this may change over time.
"

# ╔═╡ Cell order:
# ╟─17fc4346-a462-11eb-0484-7b4029b01f4d
# ╟─3ca7a180-e6ad-4734-b19f-4480de52d13c
# ╟─30339b9c-95a9-4a0e-b7d6-abd65f958e3d
# ╟─5bc0d0b1-0078-4f84-baa0-d0267f8bd7b6
# ╟─2761102e-78c2-4a4a-a7e2-c972c1bd1f4a
# ╟─420184cc-1258-4f98-880e-40dc47914011
# ╟─6b3f36e0-efa5-4010-8baa-758b4a72e75d
# ╟─78efba80-47b8-4c9c-9435-a9034e7d4053
# ╟─21a82ca1-6b0b-489a-82e8-9a19804829e7
# ╟─8e6b5906-b80a-45d7-aa8b-b73b670f65da
# ╟─7c8a60f4-896a-4fdc-9971-1fa48d9df27a
# ╟─9e29427e-4967-4d27-a885-0ac6d0319da8
# ╟─948c2330-ea9b-4350-a276-ebe940340412
# ╟─47f26447-a892-49cc-926f-ca1db8ee64b9
# ╟─d6bcbd4e-d4e6-4e7f-a7eb-41fcb7cac2d7
# ╟─7d873195-eeed-4582-a344-9e7c90eae459
# ╟─54f8cd88-f42e-42f2-a62f-8149926a21de
# ╟─1fe20c54-4c54-4f66-893b-89ff3aef08b4
# ╟─a5a29ab5-fb59-44d7-9db5-356010f00b23
# ╟─b35efb02-713a-4795-a249-84cb4d506abf
# ╟─153c2888-1cdf-41d3-bafb-8a1bbeaab03a
# ╟─6fe153a6-5262-4baa-b0f1-19cf140c3f0a
# ╟─34627d00-e9fc-4357-b626-4c5c0b1ed6b4
# ╟─7c5a19cd-9952-4ecc-adef-676031112e99
# ╟─b920bfbb-82ad-4f57-8e1c-f0001e7d77a0
# ╟─0eda0aaf-63b7-4540-af84-667ca53d57b0
# ╟─e7895abf-ad9b-4153-9927-2c5215b52562
# ╟─b3527c13-bda8-4b8d-a48d-8e20631eec52
# ╟─c9b0f626-65c3-4912-9b66-e76f7decc0f6
# ╟─aa2ee694-5f1e-47ac-93e7-f3f99023103f
# ╟─ef7eb7e2-492e-4203-9106-90342735e47c
# ╟─5e708228-6075-47f2-9153-991e8dcfb3f3
# ╟─e065e882-a5cb-4c5f-bead-434f839d909a
# ╟─48d3550d-c965-47c3-995e-d8c6fc69606e
# ╟─7df79f4b-14e0-483a-afa9-30e66cc4d537


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--gpu"
        help = "Use GPU (cuda)"
        action = :store_true
    "--profile"
        help = "Use profiler"
        action = :store_true
    "--minibatch"
        help = "Minibatch size"
        arg_type = Int
        default = 50
    "--epochs"
        help = "Epochs to train"
        arg_type = Int
        default = 100
    "--graphviz"
        help = "Dump the network structure to mnist.gv"
        action = :store_true
    "--64"
        help = "Use 64 bits"
        action = :store_true
end

const args = parse_args(ARGS, s)
const gpu = args["gpu"]
const profile = args["profile"]
const minibatch_size = args["minibatch"]
const epochs = args["epochs"]

if !gpu
  ENV["SEEP_NO_GPU"]=1
end

using MNIST, Seep, JLD

const T = args["64"] ? Float64 : Float32

const ftrain = map(T, traindata()[1]/255)

if gpu
    using CUDArt
    CUDArt.init(0)
    device(0)
    make_node(n, a::Array) = ANode(n, CudaArray(map(T, a)))
    make_node(n, x::Int...) = ANode(n, x...)
    get_data(x) = to_host(x)
    const pool = BuddyPool(CudaArray{T}(1<<24))
else
    make_node(n, a::Array) = ANode(n, map(T, a))
    make_node(n, x::Int...) = ANode(n, x...)
    get_data(x) = x
    const pool = BuddyPool(Array{T}(1<<24))
end

const n0 = 28*28
n = n0

const input = ANode("input", n, minibatch_size)
x = input

weights = ANode[]
parameters = ANode[]
for neurons = [50, 50, 2, 50, 50]
  A = make_node(randn(neurons, n)/sqrt(n))
  b = make_node(zeros(neurons, 1))
  x = A*x .+ b
  if neurons > 10
    x = tanh(x)
  end
  n = neurons

  push!(weights, A)
  push!(parameters, A)
  push!(parameters, b)
end

let neurons=28*28
  A = make_node(randn(neurons, n)/sqrt(n))
  b = make_node(zeros(neurons, 1))
  x = A*x .+ b
  n = neurons

  push!(weights, A)
  push!(parameters, A)
  push!(parameters, b)
end

const e = sum((x-input).^2)./2
const obj = e + 1e-3*reduce(+, [sum(A.^2) for A in weights])
const grad = gradients(obj)

upd = map(x->Seep.adam_fast(x, grad[x]), parameters)

const E = Float64[]

if args["graphviz"]
  #open(f->graphviz(f,  e, upd...), "mnist.gv", "w")
  open(f->graphviz(f,  e), "mnist.gv", "w")
end

Seep.devectorize(e, upd...)

const run,dict = instance(pool, e, upd...)
#const run,dict = instance(e, upd...)
const xa = dict[input]

#run() = global dict = evaluate(e, upd...)

Seep.print_stats(STDOUT, pool)

n = 0

function minibatch(ii)
    copy!(xa, ftrain[:,ii])
    run()
end

function epoch()
  perm = randperm(60000)
  ee = Float64[]
  for i in 1:minibatch_size:60000
    minibatch(perm[i+(1:minibatch_size)-1])
    push!(ee, sum(get_data(dict[e]))/minibatch_size)
  end
  push!(E, mean(ee))
end

epoch()

Profile.clear_malloc_data()

if profile
    Profile.init(delay=0.01)
    for i in 1:epochs
        @profile epoch()
    end

    data,lidict = Profile.retrieve()
    @save "mnist-profile.jld" data lidict
else
    for i in 1:epochs
        @time epoch()
        println(E[end])
    end
end

if gpu
  free(pool.pool)
  CUDArt.close(0)
end

ext: jl


text:
module Day6

using SparseArrays

M = spdiagm(1 => ones(Int32, 8))
M[7, 1] = M[9, 1] = 1

test = "3,4,3,1,2"

sim(input::String, day::Int) = begin
    pops = foldl(split(input, ","), init=zeros(Int, 9)) do acc, age
        acc[parse(Int, age)+1] += 1
        acc
    end
    sum(M^day * pops)
end

@assert sim(test, 80) == 5934
@assert sim(test, 256) == 26984457539

input = read("resources/day6_input.txt", String)
sim(input, 80) |> display
sim(input, 256) |> display

@time M^256

function do_it_really_really_fast()
    M2 = M^2
    M4 = M2^2
    M8 = M4^2
    M16 = M8^2
    M32 = Array(M16^2)
    M64 = M32^2
    M128 = M64^2
    M256 = M128^2
end

end


text:
using BacktestEngine
using Test

@testset "BacktestEngine.jl" begin
    # Write your tests here.
end


text:
"""
    Post(opinion, weight, source_agent, published_at)

Provide data structure for posts of agents within agent-based simulation.

# Examples
```julia-repl
julia>Post(0, 0, 1, 1)
Post(0.0, 0.0 , 1, 1, Int[], 0, 0, 0)
```

# Arguments
- `opinion`: opinion value between -1 and 1
- `weight`: weight value for importance of post
- `source_agent`: ID of origin agent of post
- `published_at`: publication time of post
- `seen_by`: List of agents that saw it in their feed
- `like_count`: Number of post likes
- `dislike_count`: Number of post dislikes
- `share_count`: Number of post shares

See also: [Agent](@ref), [update_feed!](@ref)
"""
mutable struct Post
    opinion::Float64
    weight::Float64
    source_agent::Int64
    published_at::Int64
    seen_by::Array{Int64, 1}
    like_count::Int64
    dislike_count::Int64
    share_count::Int64

    function Post(
        opinion,
        weight,
        source_agent,
        published_at
    )

        if opinion < -1 || opinion > 1
            error("invalid opinion value")
        end

        if weight < 0
            error("invalid weight value")
        end

        new(
            opinion,
            weight,
            source_agent,
            published_at,
            Int64[],
            0,
            0,
            0
        )
    end
end

Base.:isless(x::Post, y::Post) = x.weight < y.weight

# suppress output of include()
;


text:
using Test
using SymDiff

@test simplify(1) == 1
@test simplify(:a) == :a


text:
module CSV

using Sequoia

include("csv/utils.jl")

# Cursor

function Sequoia.Cursor(filepath)
    file = CSV.File(filepath)
    ncols = file.cols
    header = file.header
    io = open(file)
    seek(io, io.file.datapos + 1)
    coltypes = Tuple{file.types...}
    fields = Symbol[]
    # fix issue where newline char is currently captured into name of last field
    for i in 1:length(header)-1
        push!(fields, header[i])
    end
    push!(fields, header[end][1:end-1])
    fields = Tuple{fields...}
    return Sequoia.Cursor(ncols, fields, coltypes, io, 100, 0)
end

function Base.show(io::IO, curs::Sequoia.Cursor{CSV.Stream})
    println(io, summary(curs))
    row = read(curs)
    print(io, "CURRENT: ")
    for (i, val) in enumerate(row)
        isnull(val) ? print(io, "#NULL") : print(io, val.value)
        i < length(row) && print(", ")
    end
    print('\n')
    info("status: ", curs.status)
end

function Base.read(curs::Sequoia.Cursor{CSV.Stream})
    io = curs.stmt
    types = collect(typeof(curs).parameters[2].parameters)
    startpos = curs.stmt.pos
    res = Array(Any, length(types))
    for (i, typ) in enumerate(types)
        v = CSV.readfield(io, typ, curs.row, i)
        res[i] = Nullable(v[1], v[2])
    end
    curs.stmt.pos = startpos
    return tuple(res...)
end

function Sequoia.reset!(curs::Sequoia.Cursor{CSV.Stream})
    curs.io.pos = curs.io.file.datapos + 1
    curs.row = 0
end

function Sequoia.step!(curs::Sequoia.Cursor{CSV.Stream})
    io = curs.stmt
    CSV.skipn!(io, 1, io.file.quotechar, io.file.escapechar)
    curs.row += 1
    return curs
end

end # module CSV


text:
using CentroidBased
using Test

@testset "CentroidBased.jl" begin

end


text:
module AudioSchedules

import Base:
    eltype,
    empty!,
    iterate,
    IteratorEltype,
    IteratorSize,
    push!,
    show,
    write
using Base: EltypeUnknown, SizeUnknown
using Base.FastMath: sin_fast
using DataStructures: SAIterationState, SortedSet
using InfiniteArrays: Fill, ∞
using MacroTools: @capture
using PortAudio: PortAudioStream, write_buffer
import SampledSignals: nchannels, SampleBuf, samplerate
using Unitful: dB, Hz, s, ustrip

export Hz, s

const FLOAT_HERTZ = typeof(1.0Hz)
const FLOAT_PER_SECOND = typeof(1.0 / s)
const FLOAT_SECONDS = typeof(0.0s)

include("series.jl")

"""
    function Scale(ratio)

A simple wrapper that will multiply inputs by the ratio.

```jldoctest
julia> using AudioSchedules

julia> Scale(3)(2)
6
```
"""
struct Scale{Ratio}
    ratio::Ratio
end

(scale::Scale)(thing) = thing * scale.ratio

precompile(Scale{Float64}, (Float64,))

export Scale

# TODO: use FFT here?
struct SawTooth{overtones} end

"""
    SawTooth(overtones)

Build a saw-tooth wave from its partials, starting with the fundamental (1), up to
`overtones`.

To increase richness but also buziness, increase `overtones`.

```jldoctest
julia> using AudioSchedules

julia> SawTooth(3)(π / 4)
0.9185207636218614
```
"""
SawTooth(overtones) = SawTooth{overtones}()

export SawTooth

const ADJUST = 2 / pi

function (saw::SawTooth{overtones})(an_angle) where {overtones}
    ADJUST * sum(ntuple(let an_angle = an_angle
        function (overtone)
            sin_fast(overtone * an_angle) / overtone
        end
    end, Val{overtones}()))
end

precompile(SawTooth{1}, (Float64,))
precompile(SawTooth{2}, (Float64,))
precompile(SawTooth{3}, (Float64,))
precompile(SawTooth{4}, (Float64,))
precompile(SawTooth{5}, (Float64,))
precompile(SawTooth{6}, (Float64,))
precompile(SawTooth{7}, (Float64,))
precompile(SawTooth{8}, (Float64,))
precompile(SawTooth{9}, (Float64,))
precompile(SawTooth{10}, (Float64,))

export AudioSchedule

struct AudioSchedule
    sample_rate::FLOAT_HERTZ
    # instruments and the time they are first scheduled
    orchestra::Vector{Tuple{Any, FLOAT_SECONDS}}
    # which instruments in the orchestra are currently active
    activated::Vector{Bool}
    # times to turn instruments on and off
    triggers::SortedSet{Tuple{FLOAT_SECONDS, Int}}
end

"""
    AudioSchedule(; sample_rate = 44100Hz)

Create an empty `AudioSchedule`.
You can [`push!`](@ref) new synthesizers to the audio_schedule.
Provide 4 arguments to `push!`: the schedule, the synthesizer, the start time,
and an envelope that you create with [`@envelope`](@ref).

```jldoctest audio_schedule
julia> using AudioSchedules

julia> audio_schedule = AudioSchedule()
0.0 s 44100.0 Hz AudioSchedule

julia> push!(audio_schedule, Map(sin, Cycles(440Hz)), 0s, @envelope(
           0,
           Line => 1s,
           0.2,
           Line => 1s,
           0,
       ))

julia> push!(audio_schedule, Map(sin, Cycles(660Hz)), 2s, @envelope(
            0,
            Line => 1s,
            0.2,
            Line => 1s,
            0,
        ))

julia> audio_schedule
4.0 s 44100.0 Hz AudioSchedule
```

You can iterate over an `AudioSchedule`. Each element will just be a vector
of amplitudes.

```jldoctest audio_schedule
julia> length(first(audio_schedule))
44100

julia> collect(AudioSchedule())
Any[]
```

You can use write the audio_schedule directly to a `PortAudio.PortAudioStream`.
The `PortAudioStream` must have exactly 1 output channel, and a matching
sample rate.

```jldoctest audio_schedule
julia> using PortAudio: PortAudioStream

julia> PortAudioStream(0, 1, warn_xruns = false) do stream
           write(stream, audio_schedule)
       end

julia> PortAudioStream(0, 2) do stream
            write(stream, audio_schedule)
        end
ERROR: ArgumentError: PortAudioStream does not have 1 output channel
[...]

julia> PortAudioStream(0, 1, samplerate = 48000) do stream
            write(stream, audio_schedule)
        end
ERROR: ArgumentError: Sample rates of PortAudioStream (48000.0) and AudioSchedule (44100.0) do not match
[...]
```

You can save an `AudioSchedule` as a `SampledSignals.SampleBuf`.

```jldoctest audio_schedule
julia> using SampledSignals: SampleBuf

julia> saved = SampleBuf(audio_schedule)
176400-frame, 1-channel SampleBuf{Float64, 1}
4.0s sampled at 44100.0Hz
▃▄▄▄▄▅▅▅▅▅▅▅▅▆▆▆▆▆▆▆▆▆▆▆▆▆▆▅▅▅▅▅▅▅▅▄▄▄▄▃▃▄▄▄▄▅▅▅▅▅▅▅▅▆▆▆▆▆▆▆▆▆▆▆▆▆▆▅▅▅▅▅▅▅▅▄▄▄▄▃
```

You can `empty!` an `AudioSchedule` and reuse it.

```jldoctest audio_schedule
julia> empty!(audio_schedule)

julia> audio_schedule
0.0 s 44100.0 Hz AudioSchedule
```
"""
function AudioSchedule(;
    sample_rate = 44100Hz,
    orchestra = Tuple{Any, FLOAT_SECONDS}[],
    activated = Bool[],
    triggers = SortedSet{Tuple{FLOAT_SECONDS, Int}}()

)
    AudioSchedule(sample_rate, orchestra, activated, triggers)
end

precompile(AudioSchedule, ())

IteratorEltype(::Type{AudioSchedule}) = EltypeUnknown()
IteratorSize(::Type{AudioSchedule}) = SizeUnknown()

function empty!(audio_schedule::AudioSchedule)
    empty!(audio_schedule.orchestra)
    empty!(audio_schedule.activated)
    empty!(audio_schedule.triggers)
    nothing
end

precompile(empty!, (AudioSchedule,))

function samplerate(audio_schedule::AudioSchedule)
    ustrip(Hz, audio_schedule.sample_rate)
end

function nchannels(audio_schedule::AudioSchedule)
    1
end

"""
    duration(audio_schedule::AudioSchedule)

Find the duration of an `AudioSchedule` in seconds.

```jldoctest
julia> using AudioSchedules

julia> audio_schedule = AudioSchedule();

julia> push!(audio_schedule, Map(sin, Cycles(440Hz)), 0s, @envelope(
            0,
            Line => 1s,
            1,
            Line => 1s,
            0,
        ))

julia> duration(audio_schedule)
2.0 s
```
"""
function duration(audio_schedule::AudioSchedule)
    triggers = audio_schedule.triggers
    if length(triggers) > 0
        # time from first instrument turned on to last turned off
        last(triggers)[1] - first(triggers)[1]
    else
        # zero for an empty schedule
        0.0s
    end
end

precompile(duration, (AudioSchedule,))
export duration

function show(io::IO, audio_schedule::AudioSchedule)
    print(io, "$(duration(audio_schedule)) $(audio_schedule.sample_rate) AudioSchedule")
end

function envelope_macro_pieces(start_level_expression, pair_expression, end_level_expression, the_rest...)
    if !(@capture pair_expression envelopefunction_ => duration_)
        throw(ArgumentError("$pair_expression is not a pair"))
    end
    Expr(:call, (=>),
        Expr(:call, envelopefunction,
            start_level_expression,
            duration,
            end_level_expression
        ),
        Expr(:call, float, duration)
    ),
    envelope_macro_pieces(end_level_expression, the_rest...)...
end

function envelope_macro_pieces(_)
    ()
end

function envelope_macro(arguments...)
    Expr(:tuple, envelope_macro_pieces(arguments...)...)
end

"""
    @envelope(arguments...)

Create an envelope.
Start with the start amplitude (a number between 0 and 1).
Then, specify a pair, `segment_function` => `duration`, where `duration` is the duration of the segment.
Then, specify the end amplitude (again, a number between 0 and 1).
So for example,

For example, `@envelope(0, Line => 1s, 1)` will create an envelope with 1 segment.
The segment is created with `segment_function(start_level, duration, end_level) => duration`.
So, for this example, the segment will be `Line(0, 1s, 1) => 1s`.

After you finished your first segment, you can add as many more segments as you'd like.
The end level of the previous segment will be the start level of the next segment.

For example, `@envelope(0, Line => 1s, 1, Line => 1s, 0)` will create an envelope with 2 segments:

1. `Line(0, 1s, 1) => 1s`
2. `Line(1, 1s, 0) => 1s`

```jldoctest
julia> using AudioSchedules

julia> @envelope(0, Line => 1s, 1, Line => 1s, 0)
(Line(0.0, 1.0 s⁻¹) => 1.0 s, Line(1.0, -1.0 s⁻¹) => 1.0 s)

julia> @envelope(1, 2, 3)
ERROR: LoadError: ArgumentError: 2 is not a pair
[...]
```
"""
macro envelope(arguments...)
    esc(envelope_macro(arguments...))
end

export @envelope

"""
    push!(audio_schedule::AudioSchedule, synthesizer, start_time,
        start_level, shape => duration, end_level, more_segments...
    )

Add a synthesizer to a [`AudioSchedule`](@ref), where `synthesizer` is anything that supports [`make_series`](@ref),
`start_time` has units of time (like `s`), and the rest of the arguments specify the shape of the envelope.

For all envelope [`segment`](@ref), call

```
segment(shape, start_level, duration, end_level)
```

`duration` should have units of time (like `s`). For example,

```
push!(audio_schedule, synthesizer, start_time, @envelope(0, Line => 1s, 1, Line => 1s, 0))
```

will call segment twice:

```
segment(Line, 0, 1s, 1)
segment(Line, 1, 1s, 0)
```
"""
function push!(audio_schedule::AudioSchedule, (@nospecialize wave), start_time, (@nospecialize envelope))
    triggers = audio_schedule.triggers
    orchestra = audio_schedule.orchestra
    activated = audio_schedule.activated
    time_in_envelope = 0.0s
    for (shape, duration) in envelope
        synthesizer = Map(
            *,
            Skip(wave, time_in_envelope),
            shape,
        )
        push!(orchestra, (synthesizer, start_time + time_in_envelope))
        push!(activated, false)
        instrument_index = length(orchestra)
        push!(triggers, (start_time + time_in_envelope, instrument_index))
        time_in_envelope = time_in_envelope + duration
        push!(triggers, (start_time + time_in_envelope, instrument_index))
    end
    nothing
end

export push!

function iterate(audio_schedule::AudioSchedule)
    iteration = iterate(audio_schedule.triggers)
    # Done if no more triggers
    if iteration === nothing
        nothing
    else
        # Start at 0s
        trigger_iterate(audio_schedule, 0.0s, iteration...)
    end
end

precompile(iterate, (AudioSchedule,))

# TODO: precompile
function iterate(audio_schedule::AudioSchedule, (schedule_time, trigger_state))
    iteration = iterate(audio_schedule.triggers, trigger_state)
    if iteration === nothing
        nothing
    else
        trigger_iterate(audio_schedule, schedule_time, iteration...)
    end
end

precompile(iterate, (AudioSchedule, Tuple{FLOAT_SECONDS, SAIterationState}))

function trigger_iterate(
    audio_schedule,
    schedule_time,
    (trigger_time, instrument_index),
    trigger_state,
)
    sample_rate = audio_schedule.sample_rate
    activated = audio_schedule.activated
    triggers = audio_schedule.triggers
    # pull the rest of the triggers at the current time
    # that is, all triggers less than 1 sample away
    while trigger_time - schedule_time < 1 / sample_rate
        activated[instrument_index] = !activated[instrument_index]
        iteration = iterate(triggers, trigger_state)
        if iteration === nothing
            return nothing
        end
        (trigger_time, instrument_index), trigger_state = iteration
    end
    active = view(audio_schedule.orchestra, activated)
    combined = view((
        if isempty(active)
            Fill(0.0, ∞)
        else
            make_series(Map(
                +,
                Iterators.map(
                    let schedule_time = schedule_time
                        function ((instrument, instrument_start_time),)
                            Skip(instrument, schedule_time - instrument_start_time)
                        end
                    end,
                    active
                )...
            ), sample_rate)
        end
    ), 1:round(Int, (trigger_time - schedule_time) * sample_rate))
    # pull the next trigger
    activated[instrument_index] = !activated[instrument_index]
    (combined, (trigger_time, trigger_state))
end

function write_series!(buffer, series, buffer_at)
    buffer_data = buffer.data
    series_at = 0
    series_length = length(series)
    frames_per_buffer = buffer.frames_per_buffer
    left_in_buffer = frames_per_buffer - buffer_at
    if series_length >= left_in_buffer
        buffer_data[:, (buffer_at + 1):frames_per_buffer] .=
            transpose(view(series, (series_at + 1):(series_at + left_in_buffer)))
        write_buffer(buffer)
        buffer_at = 0
        series_at = series_at + left_in_buffer
        while (series_length - series_at) >= frames_per_buffer
            buffer_data .=
                transpose(view(series, (series_at + 1):(series_at + frames_per_buffer)))
            write_buffer(buffer)
            buffer_at = 0
            series_at = series_at + frames_per_buffer
        end
        number_left = series_length - series_at
        if series_length - series_at > 0
            buffer_data[:, (buffer_at + 1):number_left] .=
                transpose(view(series, (series_at + 1):series_length))
        end
        number_left
    else
        last_buffer_at = buffer_at + series_length
        buffer_data[:, (buffer_at + 1):last_buffer_at] .= transpose(series)
        last_buffer_at
    end
end

function write(
    stream::PortAudioStream,
    audio_schedule::AudioSchedule
)
    sink_messenger = stream.sink_messenger
    if nchannels(sink_messenger) != 1
        throw(ArgumentError("PortAudioStream does not have 1 output channel"))
    end
    if samplerate(stream) != samplerate(audio_schedule)
        throw(ArgumentError("Sample rates of PortAudioStream ($(samplerate(stream))) and AudioSchedule ($(samplerate(audio_schedule))) do not match"))
    end
    buffer = sink_messenger.buffer
    buffer_at = 0
    for series in audio_schedule
        buffer_at = write_series!(buffer, series, buffer_at)
    end
    if buffer_at > 0
        write_buffer(buffer, buffer_at)
    end
    nothing
end

precompile(write, (PortAudioStream, AudioSchedule))

function SampleBuf(audio_schedule::AudioSchedule)
    SampleBuf(
        reduce(vcat, audio_schedule),
        ustrip(Hz, audio_schedule.sample_rate)
    )
end

precompile(SampleBuf, (AudioSchedule,))

end


text:
module EchogramPyPlot

using PyPlot
using ColorSchemes
using EchogramColorSchemes

export echogram, eg

"""
    echogram(m; vmin = nothing, vmax = nothing, cmap=nothing, range=nothing)

Create publication quality echograms via Matplotlib.

"""
function echogram(m; vmin = nothing, vmax = nothing, cmap=nothing, range=nothing)

    y, x =size(m)

    if range==nothing
        top = 0.0
        bottom = y
    elseif isa(range, Number)
        top = 0.0
        bottom = range
    elseif isa(range, Tuple)
        top, bottom = range
    end

    stepr = 10^floor(log10(bottom-top))
    stepy = y * stepr / (bottom-top)
    yticks = 0:stepy:y
    yticklabels = top:stepr:bottom

    if length(yticks) < 3
        yticks = 0:stepy/5:y
        yticklabels = top:stepr/5:bottom
    end

    yticks = floor.(Int, yticks)
    yticklabels = floor.(Int, yticklabels)
    ax = gca()
    ax.invert_yaxis()
    ax.set_yticks(yticks)
    ax.set_yticklabels(yticklabels)

    if cmap == nothing
        cmap = addwhite(EK80) # Default
    end

    if isa(cmap, ColorScheme)
        cmap = cmap.colors
    end

    pcolormesh(m, vmin = vmin, vmax = vmax, cmap=ColorMap(cmap), rasterized=true)
end

"""
    eg(A; vmin = nothing, vmax = nothing, cmap=nothing, range=nothing)

Show a quick interactive echogram of `A`.

"""
function eg(A; vmin = nothing, vmax = nothing, cmap=nothing, range=nothing)
    figure()
    echogram(A; vmin = vmin, vmax = vmax, cmap=cmap, range=range)
    PyPlot.colorbar()
end

end # module


text:
# TODO: make a TransferMatrix type?

TransferMatrix(ψ::InfiniteMPS) = TransferMatrix(ψ, Cell(1))

function TransferMatrix(ψ::InfiniteMPS, c::Cell)
  N = nsites(ψ)
  ψᴴ = prime(linkinds, dag(ψ))
  ψᶜ = ψ[c]
  ψᶜᴴ = ψᴴ[c]
  r = unioninds(linkinds(ψ, N => N + 1), linkinds(ψᴴ, N => N + 1))
  l = unioninds(linkinds(ψ, 1 => 0), linkinds(ψᴴ, 1 => 0))
  return ITensorMap(ψᶜ, ψᶜᴴ; input_inds=r, output_inds=l)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct DedicatedHostApi <: SwaggerApi
    client::Swagger.Client
end

"""

Lists all of the dedicated hosts in the specified dedicated host group. Use the nextLink property in the response to get the next page of dedicated hosts.
Param: resourceGroupName::String (required)
Param: hostGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DedicatedHostListResult
"""
function _swaggerinternal_dedicatedHostsListByHostGroup(_api::DedicatedHostApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DedicatedHostListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/hostGroups/{hostGroupName}/hosts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "hostGroupName", hostGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function dedicatedHostsListByHostGroup(_api::DedicatedHostApi, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostsListByHostGroup(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dedicatedHostsListByHostGroup(_api::DedicatedHostApi, response_stream::Channel, resourceGroupName::String, hostGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dedicatedHostsListByHostGroup(_api, resourceGroupName, hostGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export dedicatedHostsListByHostGroup


text:
using AlphaVantage
using Test

TEST_SLEEP_TIME =  parse(Float64, get(ENV, "TEST_SLEEP_TIME", "15"))
MAX_TESTS = parse(Int64, get(ENV, "MAX_TESTS", "1"))

@testset "AlphaVantage" begin

include("stock_time_series_test.jl")
include("foreign_exchange_curency_test.jl")
include("sector_performance_test.jl")
include("technical_indicators_test.jl")
include("digital_currency_test.jl")
include("fundamentals_test.jl")

end


text:
module TestSources

using MLJBase
using Test

X = 7

Xs = source(X)
@test Xs() == X
@test Xs(8) == 8
@test elscitype(Xs) == Count
@test scitype(Xs) == MLJBase.CallableReturning{Count}
rebind!(Xs, nothing)
@test isempty(Xs)
@test Xs.scitype == Nothing

n = MLJBase.ErrorNode(ArgumentError("Stop!"))
@test_throws ArgumentError("Stop!") n(rows=1:3)
@test_throws ArgumentError("Stop!") n(1:3)

end
true