text:
include("../../common/numder/src/numder_module.jl")

module nonlineq
using LinearAlgebra
using ..numder

include("../../opti/src/line_search.jl")

struct ToleranceOptions
    max_iter::Int
    step_size_tol::Float64
end

function ToleranceOptions(;max_iter = 100, step_size_tol = 1e-15)
    ToleranceOptions(max_iter, step_size_tol)
end

function broyden(
    fcn,
    x0::Vector{T};
    tol_options::ToleranceOptions = ToleranceOptions(),
    B0 = nothing,
    damping = 1.0) where {T}

    n_dim = length(x0)

    if isnothing(B0)
        B = Matrix( one(T) * I, n_dim, n_dim )
        B = numder.jacobian_fw(fcn, x0)
    else
        B = B0
    end

    x = x0
    f = fcn(x)

    x_sol = x

    # merit_fcn(x) = norm( fcn(x) )
    iter = 0
    for outer iter = 1:tol_options.max_iter

        p = -(B \ f)
        # line search? -> damping
        s = damping * p

        if norm(s) <= tol_options.step_size_tol
            x_sol = x
            break
        end

        x = x + s
        f1 = fcn(x)
        y = f1 - f
        f = f1
        B = B + (y - B * s) * (s ./ dot(s, s))'

    end
    if iter == tol_options.max_iter
        x_sol = x
    end
    return (x_sol, iter)

end

end # nonlineq

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# load code
include("momentum.jl")
include("disk.jl")

# load calculation of occupation number fluctuations
include("fluctuations.jl")

# load correlations for different symmetries
include("correlation_lib/correlation_lib.jl")

# load tests and timers
include("test.jl")
include("timers.jl")

ext: jl


text:
baremodule DocumentationOverview

"""
    DocumentationOverview.table(module::Module; ...) -> table
    DocumentationOverview.table(fullnames::Expr; ...) -> table
    DocumentationOverview.table(apis; ...) -> table

Show the table of APIs.

The `table` object returned from this function supports `show` method with `text/plain`,
`text/markdown`, and `text/html` MIME types.

To preprocess the list of APIs, get a vector of [`API`](@ref)s from [`find`](@ref), process
it, and pass it to `DocumentationOverview.table` as an iterable of `API`s.  See also
[`list`](@ref).

See the [Gallery](@ref gallery) for example outputs.

# Extended help

# Arguments

The first argument specifies the list of APIs to be shown:

* `module::Module`: APIs discovered with `find(module::Module; ...)`
* `fullnames::Expr`: APIs are specified by the list of expression of form
  `:[a.b.c, d.e.f, ...]`
* `apis`: an iterable that produces [`API`](@ref)s

# Keyword Arguments
- `apicolumn::AbstractString = "**API**"`: The title for the API column.
- `signature`: following values are supported:
  - `nothing` (default): use the default signature.
  - `:name`: use `api.name`; i.e., `A.B.C.f(x, y)` becomes `f`.
  - `:strip_namespace`: strip the namespace part; i.e., `A.B.C.f(x, y)` becomes `f(x, y)`.
  - a callabel object: it must map an `API` to a `String` which is used as a signature.

Keyword arguments for [`find`](@ref) can also be passed with the method that takes
`module::Module` .

## Examples

```JULIA
julia> using DocumentationOverview

julia> DocumentationOverview.table(DocumentationOverview)
... a table of API printed ...
```

`DocumentationOverview.table` also takes an expression (notice `:` before `[`) to manually
list the APIs:

```julia
julia> using DocumentationOverview

julia> table = DocumentationOverview.table(:[  # ⇐ notice the `:` here
           DocumentationOverview.table,
           DocumentationOverview.API,
       ]);
```

The list of APIs retrieved using
[PublicAPI.jl](https://github.com/JuliaExperiments/PublicAPI.jl) API can also be used with
`DocumentationOverview.table`:

```julia
julia> using PublicAPI

julia> table = DocumentationOverview.table(PublicAPI.of(DocumentationOverview));
```
"""
function table end

"""
    DocumentationOverview.list(module::Module; ...) -> list
    DocumentationOverview.list(fullnames::Expr; ...) -> list
    DocumentationOverview.list(apis; ...) -> list

Show the list of APIs.

The `list` object returned from this function supports `show` method with `text/plain`,
`text/markdown`, and `text/html` MIME types.

To preprocess the list of APIs, get a vector of [`API`](@ref)s from [`find`](@ref), process
it, and pass it to `DocumentationOverview.list` as an iterable of `API`s.  See also
[`table`](@ref).

See the [Gallery](@ref gallery) for example outputs.

# Extended help

# Arguments

The first argument specifies the list of APIs to be shown:

* `module::Module`: APIs discovered with `find(module::Module; ...)`
* `fullnames::Expr`: APIs are specified by the list of expression of form
  `:[a.b.c, d.e.f, ...]`
* `apis`: an iterable that produces [`API`](@ref)s

# Keyword Arguments
- `signature`: following values are supported:
  - `nothing` (default): use the default signature.
  - `:name`: use `api.name`; i.e., `A.B.C.f(x, y)` becomes `f`.
  - `:strip_namespace`: strip the namespace part; i.e., `A.B.C.f(x, y)` becomes `f(x, y)`.
  - a callabel object: it must map an `API` to a `String` which is used as a signature.

Keyword arguments for [`find`](@ref) can also be passed with the method that takes
`module::Module` .

## Examples

```JULIA
julia> using DocumentationOverview

julia> DocumentationOverview.list(DocumentationOverview)
... a list of API printed ...
```

`DocumentationOverview.list` also takes an expression (notice `:` before `[`) to manually
list the APIs:

```julia
julia> using DocumentationOverview

julia> list = DocumentationOverview.list(:[  # ⇐ notice the `:` here
           DocumentationOverview.list,
           DocumentationOverview.API,
       ]);
```

The list of APIs retrieved using
[PublicAPI.jl](https://github.com/JuliaExperiments/PublicAPI.jl) API can also be used with
`DocumentationOverview.list`:

```julia
julia> using PublicAPI

julia> list = DocumentationOverview.list(PublicAPI.of(DocumentationOverview));
```
"""
function list end

"""
    DocumentationOverview.table_md(args...; options...) -> md::Markdown.MD

A shorthand for `Markdown.MD(DocumentationOverview.table(args...; options...))`.
"""
function table_md end

"""
    DocumentationOverview.list_md(args...; options...) -> md::Markdown.MD

A shorthand for `Markdown.MD(DocumentationOverview.list(args...; options...))`.
"""
function list_md end

"""
    DocumentationOverview.find(module::Module; [include]) -> apis::Vector{<:API}

List APIs in `module`.

# Extended help

# Keyword Arguments
- `include`: a callable of form `include(api::API) -> should_include::Bool`.  The `api`s
  are included in the output if and only if this function returns `true`.  The default
  function is equivalent to `api -> api.hasdoc`.

## Examples
```julia
julia> using DocumentationOverview

julia> for api in DocumentationOverview.find(DocumentationOverview)
           @show api.name
       end
api.name = :API
api.name = :find
api.name = :table
```
"""
function find end

"""
    DocumentationOverview.API

An object passed to `include` callback of [`table`](@ref).

`API` has the following properties settable using `Accessors.@set`:

* `module::Module`: the module in which the API is discovered.
* `name::Symbol`: the name of the API.
* `exported::Bool`: `true` iff the API is exported from `api.module`.
* `signature::String`: "signature" of the API; if `api.hasdoc` (see below), it defaults to
  the first line of the docstring; otherwise, it defaults to the fully qualified name of the
  API.

It also has the following derived properties that are not settable with `Accessors.@set`:

* `hasdoc::Bool`: indicate if the API has a docstring.
* `doc::Union{Markdown.MD,Nothing}`: docstring of the API.
* `api.value`: equivalent to `getproperty(api.module, api.value)`; may throw.
* `api.maybevalue`: evaluates to `Some(api.value)` if it does not throw or `nothing`
  otherwise
"""
abstract type API end

module Internal

using Base.Docs: DocStr, MultiDoc, catdoc, mapany, parsedoc

import Markdown
import PublicAPI
using Accessors: @set

using ..DocumentationOverview: DocumentationOverview

include("foldl.jl")
include("pkgapi.jl")
include("find.jl")
include("overview.jl")
include("table.jl")
include("list.jl")

end  # module Internal

Internal.PublicAPI.@public table table_md find API

end  # baremodule DocumentationOverview


text:
"""
    rank_to_sortperm(::NTuple{N,Int}) -> NTuple{N,Int}
Returns the `sortperm` of the stride ranks.
"""
function rank_to_sortperm(R::NTuple{N,Int}) where {N}
    sp = ntuple(zero, Val{N}())
    r = ntuple(n -> sum(R[n] .≥ R), Val{N}())
    @inbounds for n = 1:N
        sp = Base.setindex(sp, n, r[n])
    end
    sp
end
rank_to_sortperm(R) = sortperm(R)

# @inline function gep_no_offset(ptr::VectorizationBase.AbstractStridedPointer, i::Tuple)
    # VectorizationBase.gep(pointer(ptr), VectorizationBase.tdot(ptr, i, VectorizationBase.strides(ptr), VectorizationBase.nopromote_axis_indicator(ptr)))
# end
# @inline function similar_with_offset(sptr::StridedPointer{T,N,C,B,R,X,O}, ptr::Ptr{T}) where {T,N,C,B,R,X,O}
#     StridedPointer{T,N,C,B,R,X}(ptr, sptr.strd, zerotuple(Val{N}()))
# end
function view_quote(i, K, S, D, T, N, C, B, R, X, O, zero_offsets::Bool = false)
    @assert ((K == N) || isone(K))

    inds = Expr(:tuple)
    Nnew = 0
    s = Expr(:tuple)
    x = Expr(:tuple)
    o = Expr(:tuple)
    Rnew = Expr(:tuple)
    Dnew = Expr(:tuple)
    Cnew = -1
    Bnew = ifelse(iszero(B), 0, -1)
    sortp = rank_to_sortperm(R)
    still_dense = true
    densev = Vector{Bool}(undef, K)
    for k ∈ 1:K
        iₖ = Expr(:ref, :i, k)
        if i[k] === Colon
            Nnew += 1
            push!(inds.args, Expr(:ref, :o, k))
            push!(s.args, Expr(:ref, :s, k))
            push!(x.args, Expr(:ref, :x, k))
            push!(o.args, zero_offsets ? :(Zero()) : :(One()))
            if k == C
                Cnew = Nnew
            end
            if k == B
                Bnew = Nnew
            end
            push!(Rnew.args, R[k])
        else
            push!(inds.args, Expr(:call, :first, iₖ))
            if i[k] <: AbstractRange
                Nnew += 1
                push!(s.args, Expr(:call, :static_length, iₖ))
                push!(x.args, Expr(:ref, :x, k))
                push!(o.args, zero_offsets ? :(Zero()) : :(One()))
                if k == C
                    Cnew = Nnew
                end
                if k == B
                    Bnew = Nnew
                end
                push!(Rnew.args, R[k])
            end
        end
        spₙ = sortp[k]
        still_dense &= D[spₙ]
        if still_dense# && (D[spₙ])
            ispₙ = i[spₙ]
            still_dense = (ispₙ <: AbstractUnitRange) || (ispₙ === Colon)
            densev[spₙ] = still_dense
            if still_dense
                still_dense = if ((ispₙ === Colon)::Bool || (ispₙ <: Base.Slice)::Bool)
                    true
                else
                    ispₙ_len = ArrayInterface.known_length(ispₙ)
                    if ispₙ_len !== nothing
                        _sz = getfield(S, :parameters)[spₙ]
                        if _sz <: StaticInt
                            ispₙ_len == getfield(_sz, :parameters)[1]
                        else
                            false
                        end
                    else
                        false
                    end
                    false
                end
            end
        else
            densev[spₙ] = false
        end
    end
    for k ∈ 1:K
        iₖt = i[k]
        if (iₖt === Colon) || (iₖt <: AbstractVector)
            push!(Dnew.args, densev[k])
        end
    end
    quote
        $(Expr(:meta,:inline))
        sp = A.ptr
        s = A.size
        x = sp.strd
        o = sp.offsets
        new_sp = StridedPointer{$T,$Nnew,$Cnew,$Bnew,$Rnew}(gep(sp, $inds), $x, $o)
        PtrArray(new_sp, $s, Val{$Dnew}())
    end
end

@generated function Base.view(A::PtrArray{S,D,T,N,C,B,R,X,O}, i::Vararg{Union{Integer,AbstractRange,Colon},K}) where {K,S,D,T,N,C,B,R,X,O}
    view_quote(i, K, S, D, T, N, C, B, R, X, O)
end
@generated function zview(A::PtrArray{S,D,T,N,C,B,R,X,O}, i::Vararg{Union{Integer,AbstractRange,Colon},K}) where {K,S,D,T,N,C,B,R,X,O}
    view_quote(i, K, S, D, T, N, C, B, R, X, O, true)
end

@inline function Base.vec(A::PtrArray{S,D,T,N,C,0}) where {S,D,T,N,C}
    @assert all(D) "All dimensions must be dense for a vec view. Try `vec(copy(A))` instead."
    sp = StridedPointer(pointer(A), (VectorizationBase.static_sizeof(T),), (One(),))
    PtrArray(sp, (static_length(A),), Val((true,)))
end


text:
module LinOps

include("banded.jl")
include("dense.jl")
include("cuda.jl")
include("linop.jl")

end # module


text:
function lrtManual(upids,orthos,dist)
	df = DataFrame(fill(Vector{String},3),Symbol.(keys(upids)))
	# df = []
	for l in upids[:ligand]
		for r in upids[:receptor]
			for t in upids[:target]
				ligs = filter(row->row.entId==l&&row.orthoDist==string(dist),orthos)
				recs = filter(row->row.entId==r&&row.orthoDist==string(dist),orthos)
				targs = filter(row->row.entId==t&&row.orthoDist==string(dist),orthos)
				if size(ligs,1) > 0 && size(recs,1) > 0 && size(targs,1) > 0
					newrow = Dict(:ligand=>collect(skipmissing(ligs.orthoHGNC)),
						 		 :receptor=>collect(skipmissing(recs.orthoHGNC)),
								 :target=>collect(skipmissing(targs.orthoHGNC)))
				    push!(df,newrow)
				end
			end
		end
	end
	df
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractWorldModelImpl

export AbstractWorldModel

abstract AbstractWorldModel

end

module AbstractWorldModelInterfaces

export updateObjectState,
       getAllObjectStates,
       getAll,
       updateAll

function updateObjectState()
end

function getAllObjectStates()
end

function getAll()
end

function updateAll()
end

end


text:
"""
    diffusion_time_sta(f)

Get the short term approximation of the diffusion time of the time profile `f`.
"""
function diffusion_time_sta end

function diffusion_time_sta(f::PGSE)
    (; δ, Δ) = f
    out =
        4 / 35 * ((Δ + δ)^(7 / 2) + (Δ - δ)^(7 / 2) - 2Δ^(7 / 2) - 2δ^(7 / 2)) /
        (δ^2 * (Δ - δ / 3))
    out^2
end

function diffusion_time_sta(f::DoublePGSE)
    (; δ, Δ, tpause) = f
    tm = δ + tpause
    out =
        2 / 35 * (
            (2Δ + tm + δ)^(7 / 2) +
            (2Δ + tm - δ)^(7 / 2) +
            (tm + δ)^(7 / 2) +
            (tm - δ)^(7 / 2) - 2(2Δ + tm)^(7 / 2) - 2(Δ + tm + δ)^(7 / 2) -
            2(Δ + tm - δ)^(7 / 2) +
            2(Δ + δ)^(7 / 2) +
            2(Δ - δ)^(7 / 2) - 2tm^(7 / 2) + 4(Δ + tm)^(7 / 2) - 4Δ^(7 / 2) - 4δ^(7 / 2)
        ) / (δ^2 * (Δ - δ / 3))
    out^2
end

function diffusion_time_sta(f::CosOGSE)
    (; δ, Δ) = f
    n = f.nperiod
    S = sin(2π * n * Δ / δ)
    C = cos(2π * n * Δ / δ)
    FS1 = fresnels(2√n)
    FS2 = fresnels(2√(Δ * n / δ))
    FS3 = fresnels(2√((Δ + δ) * n / δ))
    FS4 = fresnels(2√((Δ - δ) * n / δ))
    FC1 = fresnelc(2√n)
    FC2 = fresnelc(2√(Δ * n / δ))
    FC3 = fresnelc(2√((Δ + δ) * n / δ))
    FC4 = fresnelc(2√((Δ - δ) * n / δ))
    out =
        3 / 8 / √(n * δ) * (
            2Δ * FC2 * C + δ * FC4 * C - δ * FC3 * C - Δ * FC3 * C - Δ * FC4 * C +
            2Δ * FS2 * S +
            δ * FS4 * S - δ * FS3 * S - Δ * FS3 * S - Δ * FS4 * S + 2δ * FC1
        ) +
        9√δ / 32π / n^(3 / 2) *
          (+2FS2 * C - FS3 * C - FS4 * C - 2FC2 * S + FC3 * S + FC4 * S + 2FS1)
    out^2
end

function diffusion_time_sta(f::SinOGSE)
    (; δ, Δ) = f
    n = f.nperiod
    S = sin(2π * n * Δ / δ)
    C = cos(2π * n * Δ / δ)
    FS1 = fresnels(2√n)
    FS2 = fresnels(2√(Δ * n / δ))
    FS3 = fresnels(2√((Δ + δ) * n / δ))
    FS4 = fresnels(2√((Δ - δ) * n / δ))
    FC1 = fresnelc(2√n)
    FC2 = fresnelc(2√(Δ * n / δ))
    FC3 = fresnelc(2√((Δ + δ) * n / δ))
    FC4 = fresnelc(2√((Δ - δ) * n / δ))
    out =
        (
            64π * n^(3 / 2) * Δ^(3 / 2) +
            64π * n^(3 / 2) * δ^(3 / 2) +
            42δ^(3 / 2) * FS2 * C +
            42δ^(3 / 2) * FS1 - 42δ^(3 / 2) * FC2 * S + 32π * n^(3 / 2) * √(Δ - δ) * δ -
            32π * n^(3 / 2) * √(Δ + δ) * Δ - 32π * n^(3 / 2) * √(Δ + δ) * δ -
            32π * n^(3 / 2) * √(Δ - δ) * Δ +
            24π * n * Δ * √δ * FC2 * C +
            24π * n * Δ * √δ * FS2 * S +
            24π * n * δ^(3 / 2) * FC1 +
            21δ^(3 / 2) * FC3 * S +
            21δ^(3 / 2) * FC4 * S - 21δ^(3 / 2) * FS3 * C - 21δ^(3 / 2) * FS4 * C +
            12π * n * δ^(3 / 2) * FS4 * S +
            12π * n * δ^(3 / 2) * FC4 * C - 12π * n * δ^(3 / 2) * FC3 * C -
            12π * n * δ^(3 / 2) * FS3 * S - 12π * n * Δ * √δ * FC4 * C -
            12π * n * Δ * √δ * FC3 * C - 12π * n * Δ * √δ * FS3 * S -
            12π * n * Δ * √δ * FS4 * S
        ) / (96π * n^(3 / 2) * δ)
    out^2
end


text:
using Plots
using FFTW

function norm(a)
    sum = 0
    for i = 1:length(a)
        sum += a[i]
    end

    a=a/sum
    return a
end

function plot_fft(res, timesteps, max_freq)
    for i = 1:timesteps
        time_name = "images/time_" * string(lpad(i-1, 4, "0")) * ".png"
        freq_name = "images/freq_" * string(lpad(i-1, 4, "0")) * ".png"
        comb_name = "images/comb_" * string(lpad(i-1, 4, "0")) * ".png"
        a = []
        if (i-1 <= div(timesteps,4))
            freq = (2+max_freq*(i-1)/(timesteps/4))
            a = [sin(freq*pi*Float64(j)/res) for j = 1:res]
        elseif (i-1 > div(timesteps,4) && i-1 <= div(timesteps,2))
            freq = (2+max_freq*(timesteps/2-(i-1))/(timesteps/4))
            a = [sin(freq*pi*Float64(j)/res) for j = 1:res]
        elseif (i-1 > div(timesteps,2) && i-1 < 3*div(timesteps,4))
            freq = (2 + max_freq*((i-1)-timesteps/2)/(timesteps/4))
            a = [sin(freq*pi*Float64(j)/res) for j = 1:res]
            b = [sin(2*freq*pi*Float64(j)/res) for j = 1:res]
            c = [sin(4*freq*pi*Float64(j)/res) for j = 1:res]
            a = a+b+c
        else
            freq = (2 + max_freq*(timesteps - (i-1))/(timesteps/4))
            a = [sin(freq*pi*Float64(j)/res) for j = 1:res]
            b = [sin(2*freq*pi*Float64(j)/res) for j = 1:res]
            c = [sin(4*freq*pi*Float64(j)/res) for j = 1:res]
            a = a+b+c
        end
        plot(a, label="Time Domain", legend = :top, linewidth = 2,
             xlabel = "Time", ylabel = "Amplitude" )
        savefig(time_name)

        b = abs2.(fft(a))
        b = norm(b)
        plot((b)[1:div(res,16)+2], label="Frequency Domain",
             legend = :top, linewidth = 2, xlabel = "Frequency",
             ylabel = "Normalized Amplitude")
        savefig(freq_name)

        run(`convert +append $time_name $freq_name $comb_name`)

    end
end

plot_fft(2000, 120, 60)


text:
# This file contains a type for interaction with result data

"""
Module containing the type to interact with result data.
"""
module Result

# Export
export ResultStruct

# Precisions and formats of numbers
using ..Prec

"""
    ResultStruct()

Instantiate this type to interact with result data.

# Data
- `Δt::`[`RT`](@ref)`=0`: sample step;
- `t::Vector{`[`RT`](@ref)`}=[]`: time array;
- `x::Vector{`[`RT`](@ref)`}=[]`: values array;
- `q::`[`RT`](@ref)`=0`: significance level;
- `threshold::`[`RT`](@ref)`=0`: signal threshold;
- `X_FFT_ABS::Vector{`[`RT`](@ref)`}=[]`: values array;
- `ν::Vector{`[`RT`](@ref)`}=[]`: periodogram frequencies array;
- `D::Vector{`[`RT`](@ref)`}=[]`: periodogram values array;
- `c::Vector{`[`RT`](@ref)`}=[]`: correlogram values array;
- `cw::Vector{`[`RT`](@ref)`}=[]`: weighted correlogram values array;
- `Dw::Vector{`[`RT`](@ref)`}=[]`: smoothed periodogram values array.

# Methods
- [`reset!`](@ref)`()`: reset an instance to default values.

"""
mutable struct ResultStruct

    # Data
    Δt::RT                # Sample step
    t::Vector{RT}         # Time array
    x::Vector{RT}         # Values array
    q::RT                 # Significance level
    threshold::RT         # Signal threshold
    X_FFT_ABS::Vector{RT} # Transformed values array
    ν::Vector{RT}         # Periodogram frequencies array
    D::Vector{RT}         # Periodogram values array
    c::Vector{RT}         # Correlogram values array
    cw::Vector{RT}        # Weighted correlogram values array
    Dw::Vector{RT}        # Smoothed periodogram values array

    # Methods
    reset!::Function # Reset an instance to default values

    # Construct an object of this type
    function ResultStruct()
        this = new(0, [], [], 0, 0, [], [], [], [], [], [])
        this.reset! = function () reset!(this) end
        return this
    end

end

# Sources
include("result_reset.jl") # Reset an instance to default values

end


text:
module FASTQDemultiplexerQC

import FASTQDemultiplexer:
    Output, Protocol, Barcode, Barcodes, InterpretedRecord,
    mergeoutput, cellidtoname, idtoname, umiid, cellid
using DataFrames
using Weave
using DataStructures

include("utils.jl")
include("polycounter.jl")

# TODO: why does it have to be mutable?
mutable struct OutputQC{N,C,U} <: Output
    results::DataFrame
    poly::Vector{PolyCounter}
    nreads::Float64
    countpoly::Bool
    protocol::Protocol{N,C,U}
end

function OutputQC{N,C,U}(protocol::Protocol{N,C,U},
                         barcodes::Barcodes;
                         outputdir::String = ".",
                         maxreads = Inf,
                         countpoly::Bool = true,
                         hamming::Bool = true,
                         polyspecs=[('A',1),('T',1)],
                         kwargs...)

    results = DataFrame(cellid = UInt[],
                        umiid = UInt[],
                        groupid = Int8[])
    poly = map(polyspecs) do specs
        PolyCounter(specs...,C)
    end
    return OutputQC{N,C,U}(results,poly,Float64(maxreads),countpoly,protocol)
end

function Base.write{N}(o::OutputQC{N},ir::InterpretedRecord{N})

    if o.nreads <= 0
        return
    else
        o.nreads-=1
    end

    push!(o.results,(cellid(ir),umiid(ir),ir.groupid))

    if o.countpoly && ir.groupid > -1
        for p in o.poly
            push!(p,ir)
        end
    end
end

function mergeoutput{N,C,U}(outputs::Vector{OutputQC{N,C,U}};
                            outputdir::String = ".",
                            hamming::Bool = false,
                            writestats::Bool = true,
                            writereport::Bool = true,
                            kwargs...)
    protocol = outputs[1].protocol
    if length(outputs) > 1
        results = vcat((o.results for o in outputs)...)

        poly = map(zip((o.poly for o in outputs)...)) do p
            merge(p...)
        end
    else
        results = first(outputs).results
        poly = first(outputs).poly
    end

    mkpath(outputdir)

    if writestats
        for p in poly
            write(outputdir,p)
        end
    end

    if writereport
        report(results,poly,outputdir,hamming,C,protocol.readnames)
    end
end

function report(results,poly,outputdir,hamming,celllen,readnames)

    args = Dict(:results=>results,
                :poly=>poly,
                :readnames=>readnames)

    if hamming
        cellids = unique(results[:cellid])
        distances = getminhammingdistance(cellids,celllen)
        args[:distances] =
            DataFrame(cellid=cellids,
                      mindistance=distances)
    end

    weave(joinpath(@__DIR__,"weave","QC.jmd"),
          args = args, out_path = outputdir)
end

end # module


text:
function set_seriescolor(seriescolor::Array, vars::Array)
    color_length = length(seriescolor)
    var_length = length(vars)
    n = Int(ceil(var_length / color_length))
    colors = repeat(seriescolor, n)[1:var_length]
    return colors
end

function _empty_plot(backend::Any)
    return Plots.plot()
end

function _dataframe_plots_internal(
    plot::Union{Plots.Plot, Nothing},
    variable::DataFrames.DataFrame,
    time_range::Array,
    backend::Any;
    kwargs...,
)
    existing_colors = ones(length(plot.series_list))
    seriescolor = permutedims(
        set_seriescolor(
            get(kwargs, :seriescolor, GR_DEFAULT),
            vcat(existing_colors, DataFrames.names(variable)),
        )[(length(existing_colors) + 1):end],
    )

    save_fig = get(kwargs, :save, nothing)
    title = get(kwargs, :title, " ")
    bar = get(kwargs, :bar, false)
    stack = get(kwargs, :stack, false)
    nofill = get(kwargs, :nofill, false)

    time_interval =
        IS.convert_compound_period(length(time_range) * (time_range[2] - time_range[1]))
    interval =
        Dates.Millisecond(Dates.Hour(1)) / Dates.Millisecond(time_range[2] - time_range[1])

    data = convert(Matrix, variable)
    labels = DataFrames.names(variable)

    isnothing(plot) && _empty_plot()
    plot_kwargs =
        Dict{Symbol, Any}(((k, v) for (k, v) in kwargs if k in SUPPORTED_EXTRA_PLOT_KWARGS))

    if stack
        plot_data = cumsum(data, dims = 2)
        if !nofill
            plot_kwargs[:fillrange] = hcat(zeros(length(time_range)), plot_data)
        end
    else
        plot_data = data
    end

    plot_kwargs[:seriescolor] = seriescolor
    plot_kwargs[:title] = title

    maxnan(a) = maximum(x -> isnan(x) ? -Inf : x, a)
    minnan(a) = minimum(x -> isnan(x) ? Inf : x, a)

    series_min, series_max = minnan(plot_data), maxnan(plot_data)
    for s in plot.series_list
        series_min = min(minnan(s[:y]), series_min)
        series_max = max(maxnan(s[:y]), series_max)
    end
    ymin = series_min <= 0.0 ? -Inf : 0.0
    ymax = series_max >= 0.0 ? Inf : 0.0

    plot_kwargs[:ylims] = get(kwargs, :ylims, (ymin, ymax))
    plot_kwargs[:ylabel] = get(kwargs, :y_label, "")
    plot_kwargs[:xlabel] = "$time_interval"
    plot_kwargs[:grid] = false

    if bar
        plot_data = sum(plot_data, dims = 1) ./ interval
        if stack
            x = nothing
            plot_data = plot_data[end:-1:1, end:-1:1]
            plot_kwargs[:lab] = hcat(string.(labels)...)[end:-1:1, end:-1:1]
            plot_kwargs[:seriescolor] = seriescolor[:, length(labels):-1:1]
            plot_kwargs[:legend] = :outerright
        else
            x = labels
            plot_data = permutedims(plot_data)
            plot_kwargs[:lab] = hcat(string.(labels)...)
            plot_kwargs[:seriescolor] = permutedims(seriescolor)
            plot_kwargs[:legend] = false
        end
        plot_func = nofill ? Plots.hline! : Plots.bar!
        p = plot_func(x, plot_data; plot_kwargs...)
    else
        plot_kwargs[:lab] = hcat(string.(labels)...)
        plot_kwargs[:linetype] = get(kwargs, :stair, false) ? :steppost : :line
        plot_kwargs[:xtick] = [time_range[1], last(time_range)]
        plot_kwargs[:legend] = :outerright

        p = Plots.plot!(time_range, plot_data; plot_kwargs...)
    end
    get(kwargs, :set_display, false) && display(p)
    title = title == " " ? "dataframe" : title
    !isnothing(save_fig) &&
        save_plot(p, joinpath(save_fig, "$(title).png"), backend; kwargs...)
    return p
end

function save_plot(plot::Plots.Plot, filename::String, backend::Any; kwargs...)
    Plots.savefig(plot, filename) # TODO: add kwargs support
end


text:
mutable struct SimpleRooms{S<:Integer, A<:Integer, R<:Real} <: DiscEnv{S, A, R}
    current::S
    observationspace::ObservationSpace{S}
    actionspace::ActionSpace{Symbol}
    transitions::Dict{S, Vector{S}}
    rewards::Vector{R}
    function SimpleRooms{S, A, R}() where {S<:Integer, A<:Integer, R<:Real}
        observationspace = ObservationSpace{S}(16)
        actionspace = ActionSpace([:up, :down, :left, :right])
        transitions = Dict{S, Vector{S}}(
            1 => [2, 5], 2 => [1, 3, 6], 3 => [2, 4, 7], 4 => [3, 8],
            5 => [1, 6, 9], 6 => [2, 5], 7 => [3, 8], 8 => [4, 7, 12],
            9 => [5, 10, 13], 10 => [9, 14], 11 => [12, 15], 12 => [8, 11, 16],
            13 => [9, 14], 14 => [10, 13, 15], 15 => [11, 14, 16], 16 => [12, 15]
        )
        rewards = zeros(R, 16)
        rewards[16] = one(R)
        new{S, A, R}(1, observationspace, actionspace, transitions, rewards)
    end
end

SimpleRooms() = SimpleRooms{Int64, Int64, Float64}()

function singleStep(s, a, t, r)
    if a == :up && (s-4 in t[s])
        return s-4, r[s-4]
    elseif a == :down && (s+4 in t[s])
        return s+4, r[s+4]
    elseif a == :left && (s-1 in t[s])
        return s-1, r[s-1]
    elseif a == :right && (s+1 in t[s])
        return s+1, r[s+1]
    end
    return s+0, zero(eltype(r))
end

step!(env::SimpleRooms, action::Integer) = step!(env, env.actionspace.actions[action])

function step!(env::SimpleRooms, action::Symbol)
    env.current, r = singleStep(env.current, action, env.transitions, env.rewards)
    return env.current, r, r == 1.0
end

reset!(env::SimpleRooms) = (env.current = 1)


text:
# simple discrete hmm with testable output
A = eye(2)
B_matrix = [0.0  1.0  0.0 ;
            1.0  0.0  0.0 ]

hmm = HMM(A,B_matrix)
hmm.p = vec([1 0])

# state sequence s, observation sequence o
n_obs = 10
s,o = generate(hmm,n_obs)
@test(length(s) == length(o) == n_obs)
@test(all(s .== 1))
@test(all(o .== 2))

# check forward algorithm on perfect sequence
alpha,p_obs = forward(hmm,o; scaling=false)
@test(p_obs == 1)
@test(all(alpha[:,1] .== 1))
@test(all(alpha[:,2] .== 0))

# test forward with scaling
alpha,log_p_obs,coeff = forward(hmm,o)
@test(log_p_obs == 0.0)
@test(all(alpha[:,1] .== 1))
@test(all(alpha[:,2] .== 0))
@test(all(coeff .== 1.0))

# check backward algorithm on perfect sequence
beta = backward(hmm,o)
@test(all(beta[end,:] .== 1))
@test(all(beta[:,1] .== 1))

# test backward with scaling
beta = backward(hmm,o;scale_coeff=coeff)
@test(all(beta[end,:] .== 1))
@test(all(beta[:,1] .== 1))

# check viterbi algorithm on perfect sequence
@test(all(viterbi(hmm,o) .== s))

# check forward-backward algorithm on impossible sequence
alpha,p_obs = forward(hmm,o+1; scaling=false)
beta = backward(hmm,o+1)
@test(p_obs == 0)
@test(all(alpha .== 0))
@test(all(beta[1:end-1,:] .== 0))

######
# Compare output to Michael Hamilton's "dishonest casino" example
# see Hamilton's HMM implementation in python:
# http://www.cs.colostate.edu/~hamiltom/code.html#python-hidden-markov-model

A = [ 0.99 0.01 ;
      0.01 0.99 ]
B = [ 1/6   1/6   1/6   1/6   1/6   1/6 ;
      1/10  1/10  1/10  1/10  1/10  1/2 ]
p = [0.5, 0.5]
hmm = HMM(A,B,p)

@test(hmm.n == 2)
for i = 1:hmm.n
  @test(typeof(hmm.B[i]) == Categorical)
  @test(hmm.B[i].K == 6)
end
@test(length(hmm.p) == hmm.n)

o = [1, 2, 1, 6, 6] # observation sequence

# check forward/backward (without scaling)
true_alpha = [ 8.3333e-02   1.3833e-02   2.2909e-03   3.7885e-04   6.295e-05  ;
               5.0000e-02   5.0333e-03   5.1213e-04   2.6496e-04   1.3305e-04 ]

true_beta = [ 8.9630e-04   5.2841e-03   3.0533e-02   1.7000e-01   1.0000e+00 ;
              2.4262e-03   2.4418e-02   2.4613e-01   4.9667e-01   1.0000e+00 ]

alpha,p_obs = forward(hmm,o; scaling=false)
beta = backward(hmm,o)

@test(round(p_obs,6) == 0.000196)
@test(all(round(alpha,6) .== round(true_alpha',6)))
@test(all(round(beta,5) .== round(true_beta',5)))

# check forward/backward (with scaling)
alpha_scaled = [ 0.625   0.7332  0.8173  0.5884  0.3212 ;
                 0.375   0.2668  0.1827  0.4116  0.6788 ]

beta_scaled = [  4.5729   3.5946   2.9391   2.4312  3.2847 ;
                 12.3785  16.6108  23.6922  7.1028  3.2847 ]

alpha,log_p_obs,coeff = forward(hmm,o; scaling=true)
beta = backward(hmm,o; scale_coeff=coeff)

@test(log_p_obs == log(p_obs))
@test(all(round(alpha,4) .== round(alpha_scaled',4)))
@test(all(round(beta,4) .== round(beta_scaled',4)))

# check viterbi
true_vit = [2,2,2,2,2]
@test(all(viterbi(hmm,o) .== true_vit))
@test(all(viterbi(hmm,[1,2,3,4,5]) .== ones(Int,5)))

# another check for viterbi
hmm.A = [ 0.7 0.3 ;
          0.3 0.7 ]
o = [ 1, 2, 1, 6, 6, 3, 3, 2, 1 ]
true_vit = [ 1, 1, 1, 2, 2, 1, 1, 1, 1 ]
@test(all(viterbi(hmm,o) .== true_vit))

# check Baum-Welch algorithm
o = [ 1, 2, 1, 6, 6, 3, 3, 2, 1 ]
expected_result_A = [ 0.75  0.25 ;
                      0.25  0.75 ]
expected_result_B = [ 0.6  0.4  0.0  0.0  0.0  0.0 ;
                      0.0  0.0  0.5  0.0  0.0  0.5 ]
expected_result_p = [1.0, 0.0]

hmm1 = deepcopy(hmm) # check without scaling alpha/beta
fit!(hmm1, o; max_iter=60, scaling=false)

@test(all(round(hmm1.A,4) .== expected_result_A))
@test(all(round(hmm1.B[1].p,4) .== expected_result_B[1,:]'))
@test(all(round(hmm1.B[2].p,4) .== expected_result_B[2,:]'))
@test(all(round(hmm1.p,4) .== expected_result_p))

hmm2 = deepcopy(hmm) # check with scaling alpha/beta
fit!(hmm2, o; max_iter=60, scaling=true)
@test(all(round(hmm2.A,4) .== expected_result_A))
@test(all(round(hmm2.B[1].p,4) .== expected_result_B[1,:]'))
@test(all(round(hmm2.B[2].p,4) .== expected_result_B[2,:]'))
@test(all(round(hmm2.p,4) .== expected_result_p))

# hmm3 = deepcopy(hmm) # check with multiple sequences
# oo = (Vector{Int})[o;o] # two copies of the same sequence
# baum_welch!(hmm3, oo; max_iter=60, tol=NaN, scaling=true)
# @test(all(round(hmm3.A,4) .== expected_result_A))
# @test(all(round(hmm3.B,4) .== expected_result_B))
# @test(all(round(hmm3.p,4) .== expected_result_p))

# ## Check Baum-Welch from random initialization
# # Model used to generate observation sequences:
# # A = [ 0.7 0.3 ;
# #       0.4 0.6 ]
# # B = [ 1/2        1/3   1/6+eps() ;
# #       1/6+eps()  1/3   1/2       ]
# # p = [0.5, 0.5]
# hmm = dHMM(2,3)
# o = [1,3,1,1,1,2,3,1,2,1,3,2,1,3,2,2,2,2,2,2,2,1,1,1,2,1,1,2,3,2,1,2,3,3,2,1,1,3,2,3,2,2,2,3,2,3,3,2,1,2]

# ch = baum_welch!(hmm, o; max_iter=1000, tol=1e-7)

# @test(all(round(sum(hmm.A,2),15) .== 1.0))
# @test(all(round(sum(hmm.B,2),15) .== 1.0))
# @test(round(sum(hmm.p),15) == 1.0)

# # log-liklihood should always increase under baum_welch
# @test(all(diff(ch) .>= 0.0))

# # log-liklihood should plateau (test tolerance parameter)
# if length(ch)<1000
#   @test(ch[end] - ch[end-1] < 1e-7)
# end

# ## Check Baum-Welch on multiple sequences
# o1 = [3,3,3,1,1,1,1,3,2,2]
# o2 = [1,3,3,1,3,2,2,3,1,1]
# o3 = [2,3,3,3,3,1,2,1,3,2]
# o4 = [2,1,1,2,3,1,3,2,3,1]
# o5 = [3,3,3,1,3,2,1,3,1,1]
# o6 = [2,1,1,2,2,2,2,3,2,3,1,1,1,2,3,3,2,3,3,1]
# o7 = [1,1,3,1,3,3,2,3,2,3,2,1,1,2,2,2,2,3,1,2]
# o8 = [1,3,2,1,3,2,2,2,2,1,2,1,3,2,2,2,2,1,2,2]
# o9 = [3,2,2,1,3,2,1,2,3,3,1,3,1,2,2,3,2,1,3,3]
# o10 = [2,2,1,3,2,2,1,3,1,1,2,3,3,1,1,3,3,1,1,3]

# seq1 = (Vector{Int})[o1; o2; o3; o4; o5]
# seq2 = (Vector{Int})[o6; o7; o8; o9; o10]
# seq3 = (Vector{Int})[o1; o2; o3; o4; o5; o6; o7; o8; o9; o10]

# # First set of sequences
# hmm = dHMM(2,3) # re-initialize with random params
# ch = baum_welch!(hmm, seq1; max_iter=1000, tol=1e-7)
# @test(all(round(sum(hmm.A,2),15) .== 1.0))
# @test(all(round(sum(hmm.B,2),15) .== 1.0))
# @test(round(sum(hmm.p),15) == 1.0)
# @test(all(diff(ch) .>= -1e-10))
# if length(ch)<1000
#   @test(ch[end] - ch[end-1] < 1e-7)
# end

# # Second set of sequences
# hmm = dHMM(2,3) # re-initialize with random params
# ch = baum_welch!(hmm, seq2; max_iter=1000, tol=1e-7)
# @test(all(round(sum(hmm.A,2),15) .== 1.0))
# @test(all(round(sum(hmm.B,2),15) .== 1.0))
# @test(round(sum(hmm.p),15) == 1.0)
# #@test(all(diff(ch) .>= -1e-10))
# if length(ch)<1000
#   @test(ch[end] - ch[end-1] < 1e-7)
# end

# # All sequences
# hmm = dHMM(2,3) # re-initialize with random params
# ch = baum_welch!(hmm, seq3; max_iter=100, tol=NaN)
# @test(all(round(sum(hmm.A,2),15) .== 1.0))
# @test(all(round(sum(hmm.B,2),15) .== 1.0))
# @test(round(sum(hmm.p),15) == 1.0)

# These tests fail...
#    @test(all(diff(ch) .>= 0.0))
#    @test(ch[end] - ch[end-1] < 1e-7)
# potentially because having different length observation sequences
# violates some assumptions?


text:
using GaussQuadrature, LinearAlgebra

# First construct the mapping for a spherical sector

# parameters:
nn = 12  # number of sectors
φ = π / nn # half angle of sector
r1 = 0.8 # inner radius
r2 = 1.0 # outer radius

# mappings
# annulus mappings
R(ξ¹; r1 = r1, r2 = r2) = (r2 - r1) * (ξ¹ + 1) / 2 + r1
Θ(ξ²; θ1 = π / 2 - φ, θ2 = π / 2 + φ) = (θ2 - θ1) * (ξ² + 1) / 2 + θ1

# annulus to cartesian mappings
x_mapping(r, θ) = r * cos(θ)
y_mapping(r, θ) = r * sin(θ)

# Now construct the discretization
# Construct Gauss-Lobatto points
N1 = 1
ξ1vec, ω1 = GaussQuadrature.legendre(N1 + 1, GaussQuadrature.both)
ω1 = reshape(ω1, (N1 + 1, 1))
N2 = 2
ξ2vec, ω2 = GaussQuadrature.legendre(N2 + 1, GaussQuadrature.both)
ω2 = reshape(ω2, (1, N2 + 1))

# Construct Differentiation Matrices
vec = ξ1vec
ξp = [vec[i]^(j - 1) for i in 1:length(vec), j in 1:length(vec)]
dξp = [(j - 1) * vec[i]^(j - 2) for i in 1:length(vec), j in 1:length(vec)]
∂ξ¹ = dξp / ξp

vec = ξ2vec
ξp = [vec[i]^(j - 1) for i in 1:length(vec), j in 1:length(vec)]
dξp = [(j - 1) * vec[i]^(j - 2) for i in 1:length(vec), j in 1:length(vec)]
∂ξ² = dξp / ξp

# Calculate the grid points in polar coordinates
rvec = R.(ξ1vec)
θvec = Θ.(ξ2vec)
# convert to cartesian
x_positions = [x_mapping(r, θ) for r in rvec, θ in θvec]
y_positions = [y_mapping(r, θ) for r in rvec, θ in θvec]

# Now do the discrete Calculus
∂x∂ξ¹ = copy(x_positions)
∂y∂ξ¹ = copy(x_positions)
for j in 1:size(∂x∂ξ¹)[2]
    ∂x∂ξ¹[:, j] = ∂ξ¹ * x_positions[:, j]
    ∂y∂ξ¹[:, j] = ∂ξ¹ * y_positions[:, j]
end

∂x∂ξ² = copy(x_positions)
∂y∂ξ² = copy(x_positions)
for i in 1:size(∂x∂ξ²)[1]
    ∂x∂ξ²[i, :] = ∂ξ² * x_positions[i, :]
    ∂y∂ξ²[i, :] = ∂ξ² * y_positions[i, :]
end
jacobian = zeros(size(x_positions)..., (2, 2)...)

# Construct convenient object for entries, The columns are the covariant vectors
jacobian[:, :, 1, 1] .= ∂x∂ξ¹
jacobian[:, :, 1, 2] .= ∂x∂ξ²
jacobian[:, :, 2, 1] .= ∂y∂ξ¹
jacobian[:, :, 2, 2] .= ∂y∂ξ²

detJ = [det(jacobian[i, j, :, :]) for i in 1:length(ξ1vec), j in 1:length(ξ2vec)]
M = detJ .* ω1 .* ω2
exact_area = (r2^2 - r1^2) * φ
approx_area = sum(M)

wrongness = (exact_area - approx_area) / exact_area
println("The error in computing the area is ", wrongness)

# Construct contravariant basis numerically, these are the face normals
ijacobian = copy(jacobian) # the rows are the contravariant vectors
for i in 1:length(ξ1vec), j in 1:length(ξ2vec)
    tmp = inv(jacobian[i, j, :, :])
    ijacobian[i, j, :, :] .= tmp
end
ijacobian[1, 1, :, :]
x_positions[1, 1], y_positions[1, 1]

# face 2 is the linear side
approx_vec = ijacobian[1, 1, 2, :] ./ norm(ijacobian[1, 1, 2, :])
exact_vec = [-y_positions[1, 1], x_positions[1, 1]] ./ norm([y_positions[1, 1], -x_positions[1, 1]])
println("angle face error ", norm(approx_vec - exact_vec) / norm(exact_vec))
# face 1 is the curvy side
approx_vec = ijacobian[1, 1, 1, :] ./ norm(ijacobian[1, 1, 1, :])
exact_vec = [x_positions[1, 1], y_positions[1, 1]] ./ norm([x_positions[1, 1], y_positions[1, 1]])
println("radial face error ", norm(approx_vec - exact_vec) / norm(exact_vec))

# Analytic Metrics
M, N = size(detJ)
a₁ᴬ = zeros(M, N, 2)
a₂ᴬ = zeros(M, N, 2)
jacobianᴬ = zeros(M, N, 2, 2)
ijacobianᴬ = zeros(M, N, 2, 2)
detJᴬ = zeros(M, N)
a¹ᴬ = zeros(M, N, 2)
a²ᴬ = zeros(M, N, 2)
for i in 1:M, j in 1:N
    a₁ᴬ[i, j, 1] = (rvec[2] - rvec[1]) / 2 * cos((θvec[end] - θvec[1]) * (ξ2vec[j] + 1) / 2 + θvec[1])
    a₁ᴬ[i, j, 2] = (rvec[2] - rvec[1]) / 2 * sin((θvec[end] - θvec[1]) * (ξ2vec[j] + 1) / 2 + θvec[1])

    a₂ᴬ[i, j, 1] = -((rvec[2] - rvec[1]) / 2 * (ξ1vec[i] + 1) + rvec[1]) * sin((θvec[end] - θvec[1]) * (ξ2vec[j] + 1) / 2 + θvec[1]) * (θvec[end] - θvec[1]) / 2
    a₂ᴬ[i, j, 2] = ((rvec[2] - rvec[1]) / 2 * (ξ1vec[i] + 1) + rvec[1]) * cos((θvec[end] - θvec[1]) * (ξ2vec[j] + 1) / 2 + θvec[1]) * (θvec[end] - θvec[1]) / 2

    jacobianᴬ[i, j, 1, 1] = a₁ᴬ[i, j, 1]
    jacobianᴬ[i, j, 2, 1] = a₁ᴬ[i, j, 2]
    jacobianᴬ[i, j, 1, 2] = a₂ᴬ[i, j, 1]
    jacobianᴬ[i, j, 2, 2] = a₂ᴬ[i, j, 2]
    detJᴬ[i, j] = det(jacobianᴬ[i, j, :, :])
    ijacobianᴬ[i, j, :, :] = inv(jacobianᴬ[i, j, :, :])

    a¹ᴬ[i, j, 1] = ijacobianᴬ[i, j, 1, 1]
    a¹ᴬ[i, j, 2] = ijacobianᴬ[i, j, 1, 2]

    a²ᴬ[i, j, 1] = ijacobianᴬ[i, j, 2, 1]
    a²ᴬ[i, j, 2] = ijacobianᴬ[i, j, 2, 2]
end

# Free-Stream Property
i = 1
j = 1
ijacobian[i, j, :, :]
jacobian[i, j, :, :]

a₁ = jacobian[:, :, :, 1];
a₂ = jacobian[:, :, :, 2];

a¹ = ijacobian[:, :, 1, :];
a² = ijacobian[:, :, 2, :];

dot(a¹[1, 1, :], a₁[1, 1, :])
dot(a¹[1, 1, :], a₂[1, 1, :])
dot(a²[1, 1, :], a₁[1, 1, :])
dot(a²[1, 1, :], a₂[1, 1, :])
J = reshape(detJ, (size(a¹)[1:2]..., 1))

# choose vector field in cartesian coordinates
c = [1 1]
# project onto contravariant components
u¹ = [(c*a¹[i, j, :])[1] for i = 1:M, j = 1:N]
u² = [(c*a²[i, j, :])[1] for i = 1:M, j = 1:N]

divu = zeros(size(u¹))
detiJ = 1 ./ detJ
for i = 1:M, j = 1:N
    for ii = 1:M
        divu[i, j] += detiJ[i, j] * ∂ξ¹[i, ii] * detJ[ii, j] * u¹[ii, j]
    end
end
for j = 1:N, i = 1:M
    for jj = 1:N
        divu[i, j] += detiJ[i, j] * ∂ξ²[j, jj] * detJ[i, jj] * u²[i, jj]
    end
end
println("numerically computed divergence with numerical metrics")
display(divu)

# choose vector field in cartesian coordinates
c = [1 1]
# project onto contravariant components
u¹ = [(c*a¹ᴬ[i, j, :])[1] for i = 1:M, j = 1:N]
u² = [(c*a²ᴬ[i, j, :])[1] for i = 1:M, j = 1:N]

divu = zeros(size(u¹))
detiJᴬ = 1 ./ detJᴬ
for i = 1:M, j = 1:N
    for ii = 1:M
        divu[i, j] += detiJᴬ[i, j] * ∂ξ¹[i, ii] * detJᴬ[ii, j] * u¹[ii, j]
    end
end
for j = 1:N, i = 1:M
    for jj = 1:N
        divu[i, j] += detiJᴬ[i, j] * ∂ξ²[j, jj] * detJᴬ[i, jj] * u²[i, jj]
    end
end
println("numerically computed divergence with analytic metrics")
display(divu)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
"""
Code for combining spectra (currently just sums them)

Author: Eric Ford
Created: August 2020
Contact: https://github.com/eford/
"""

"""   bin_times_consecutive( times, n )
Computes mean times from conseuctive bins of n times (to go with bin_consecutive_spectra).
Returns floor(length(times)/n) elements.

WARNING:  Simply takes consecutive times, so some bins may be from spectra that weren't taken close together.
TODO:  Create version that pays attention to timestamps.
"""
function bin_times_consecutive(times::AT, n::Integer) where { T<:Real, AT<:AbstractVector{T} }
    local num_in = length(times)
    @assert num_in >= n
    local num_out = floor(Int,num_in//n)
    local times_out = Vector{eltype(times)}(undef,num_out)
    local time_idx_out  = Array{UnitRange,1}(undef,num_out)
    local idx_start = 1
    for i in 1:num_out
      idx_stop = min(idx_start+n-1,num_in)
      times_out[i] = mean(times[idx_start:idx_stop])
      time_idx_out[i] = idx_start:idx_stop
      idx_start = idx_stop + 1
    end
    return times_out
end

bin_rvs_consecutive(rvs::AT, n::Integer) where { T<:Real, AT<:AbstractVector{T} } = bin_times_consecutive(rvs,n)

"""   bin_spectra_nightly
Bins spectra from a SpectralTimeSeriesCommonWavelengths object with a maximum spacing between observation times
"""
function bin_spectra_nightly(spectra::AbstractSpectralTimeSeriesCommonWavelengths, times::AT) where { T<:Real, AT<:AbstractVector{T} }
    Δt_threshold = 0.5
    bin_spectra_max_Δt(spectra,times,Δt_threshold)
end

function bin_times_nightly(times::AT)  where { T<:Real, AT<:AbstractVector{T} }
    Δt_threshold = 0.5
    bin_times_max_Δt(times,Δt_threshold)
end

function bin_rvs_nightly(;times::AT2, rvs::AT1)  where { T1<:Real, AT1<:AbstractVector{T1}, T2<:Real, AT2<:AbstractVector{T2}  }
    Δt_threshold = 0.5
    bin_rvs_max_Δt(times=times,rvs=rvs,Δt_threshold=Δt_threshold)
end

function bin_times_and_rvs_nightly(;times::AT1, rvs::AT2 )  where { T1<:Real, AT1<:AbstractVector{T1}, T2<:Real, AT2<:AbstractVector{T2}  }
    Δt_threshold = 0.5
    bin_times_and_rvs_max_Δt(times=times, rvs=rvs,  Δt_threshold=Δt_threshold)
end

function bin_times_max_Δt(times::AT, Δt_threshold::Real = 0.5)  where { T<:Real, AT<:AbstractVector{T} }
    local bin_idx = make_bin_indices_for_binning_max_Δt(times,Δt_threshold=Δt_threshold)
    local bin_labels = unique(bin_idx)
    local binned_times = zeros(length(bin_labels))
    for (i,label) in enumerate(bin_labels)
        binned_times[i] = mean(times[findall(isequal(label),bin_idx)])
    end
    return binned_times
end

function bin_rvs_max_Δt(;times::AT1, rvs::AT2, Δt_threshold::Real = 0.5)  where { T1<:Real, AT1<:AbstractVector{T1}, T2<:Real, AT2<:AbstractVector{T2}  }
    local bin_idx = make_bin_indices_for_binning_max_Δt(times,Δt_threshold=Δt_threshold)
    local bin_labels = unique(bin_idx)
    local binned_rvs = zeros(length(bin_labels))
    for (i,label) in enumerate(bin_labels)
        binned_rvs[i] = mean(rvs[findall(isequal(label),bin_idx)])
    end
    return binned_rvs
end

function bin_times_and_rvs_max_Δt(;times::AT1, rvs::AT2, Δt_threshold::Real = 0.5)  where { T1<:Real, AT1<:AbstractVector{T1}, T2<:Real, AT2<:AbstractVector{T2}  }
    local bin_idx = make_bin_indices_for_binning_max_Δt(times,Δt_threshold=Δt_threshold)
    local bin_labels = unique(bin_idx)
    local binned_times = zeros(length(bin_labels))
    local binned_rvs = zeros(length(bin_labels))
    local num_obs_in_bin = zeros(length(bin_labels))
    for (i,label) in enumerate(bin_labels)
        local idx = findall(isequal(label),bin_idx)
        binned_times[i] = mean(times[idx])
        binned_rvs[i] = mean(rvs[idx])
        num_obs_in_bin[i] = length(idx)
    end
    return (times=binned_times, rvs=binned_rvs, num_obs_in_bin=num_obs_in_bin)
end

function make_bin_indices_for_binning_max_Δt(times::AT; Δt_threshold::Real = 0.5)  where { T<:Real, AT<:AbstractVector{T} }
    @assert length(times) >= 1
    n = length(times)
    bin_indices = ones(Int,n)
    bin = 1
    bin_start_time = times[1]
    for i in 2:n
        if times[i]-bin_start_time > Δt_threshold
            bin += 1
            bin_start_time = times[i]
        end
        bin_indices[i] = bin
    end
    return bin_indices
end

"""   bin_spectra_max_Δt
Bins spectra from a SpectralTimeSeriesCommonWavelengths object with a maximum spacing between observation times
"""
function bin_spectra_max_Δt(spectra::AbstractSpectralTimeSeriesCommonWavelengths, times::AT, Δt::Real) where { T<:Real, AT<:AbstractVector{T} }
  local num_in = size(spectra.flux,2)
  @assert length(times) >= 1
  @assert Δt>=zero(Δt)
  local bin_idx = make_bin_indices_for_binning_max_Δt(times,Δt_threshold=Δt)
  local bin_labels = unique(bin_idx)
  local num_out = length(bin_labels)
  @assert num_in >= num_out
  local binned_times = zeros(num_out)
  local flux_out = Array{eltype(spectra.flux),2}(undef,size(spectra.flux,1),num_out)
  local var_out  = Array{eltype(spectra.flux),2}(undef,size(spectra.var ,1),num_out)
  local time_ave = Array{eltype(times),1}(undef,num_out)
  local time_idx_out = Array{UnitRange,1}(undef,num_out)
  for (i,label) in enumerate(bin_labels)
      idx_to_bin = findall(isequal(label),bin_idx)
      time_ave[i] = mean(times[idx_to_bin])
      flux_out[:,i] .= vec(sum(spectra.flux[:,idx_to_bin],dims=2))
      var_out[:,i]  .= vec(sum((spectra.var[:,idx_to_bin]),dims=2))
      if length(idx_to_bin) > 1
          @assert all(idx_to_bin[2:end].-idx_to_bin[1:end-1] .== 1)
      end
      time_idx_out[i] = first(idx_to_bin):last(idx_to_bin)
  end
  metadata_out = MetadataT(:time_idx=>time_idx_out, :bjd=>time_ave)
  return typeof(spectra)(spectra.λ,flux_out,var_out,spectra.chunk_map,spectra.inst,metadata_out)
end

"""   bin_spectra_consecutive
Bins consecutive spectra from a SpectralTimeSeriesCommonWavelengths object

WARNING:  Simply takes consecutive spectra, so some bins may be from spectra that weren't taken close together.
TODO:  Create version that pays attention to timestamps.
"""
function bin_spectra_consecutive(spectra::AbstractSpectralTimeSeriesCommonWavelengths, n::Integer)
  local num_in = size(spectra.flux,2)
  @assert num_in >= n
  local num_out = floor(Int,num_in//n)
  local flux_out = Array{eltype(spectra.flux),2}(undef,size(spectra.flux,1),num_out)
  local var_out  = Array{eltype(spectra.flux),2}(undef,size(spectra.flux,1),num_out)
  local time_idx_out  = Array{UnitRange,1}(undef,num_out)
  #idx_binned = map(i->1+(i-1)*obs_per_bin:obs_per_bin*i,1:num_obs_binned)

  local idx_start = 1
  for i in 1:num_out
    idx_stop = min(idx_start+n-1,num_in)
    flux_out[:,i] .= vec(sum(spectra.flux[:,idx_start:idx_stop],dims=2))
    var_out[:,i]  .= vec(sum((spectra.var[:,idx_start:idx_stop]),dims=2))
    time_idx_out[i] = idx_start:idx_stop
    idx_start = idx_stop + 1
  end
  metadata_out = MetadataT(:time_idx=>time_idx_out)
  return typeof(spectra)(spectra.λ,flux_out,var_out,spectra.chunk_map,spectra.inst,metadata_out)
end

"""   bin_times()
Bins times from a SpectralTimeSeriesCommonWavelengths object using the groupings from time_idx in the metadata.
"""
function bin_times(spectra::AbstractSpectralTimeSeriesCommonWavelengths, times::AT, n::Integer) where { T<:Real, AT<:AbstractVector{T} }
    @assert haskey(spectra.metadata,:time_idx)
    map(idx->mean(times[idx]), spectra.metadata[:time_idx])
end

"""  rms_rv_within_night(times, rvs)
Return RMS of RVs taken within the same night
"""
function rms_rvs_within_night(;times::AbstractVector{T1}, rvs::AbstractVector{T2}) where { T1<:Real, T2<:Real }
    Δt_threshold = 0.5
    bin_idx = make_bin_indices_for_binning_max_Δt(times,Δt_threshold=Δt_threshold)
    bin_labels = unique(bin_idx)
    sum_rms = 0
    sum_weight = 0
    for (i,label) in enumerate(bin_labels)
        idx = findall(isequal(label),bin_idx)
        if length(idx) <= 1    continue    end
        weight = length(idx)-1
        sum_rms += std(rvs[idx])*weight
        sum_weight += weight
    end
    @assert sum_weight > 0
    within_night_rms = sum_rms/sum_weight
    return within_night_rms
end

""" make_template_spectra( chunk_list_timeseries, [ options ] )
Combine portions of spectra in a ChunkListTimeseries into a template stored as a SpectralTimeSeriesCommonWavelengths.
Optionally remove radial velocities before combining spectra based on the :rv_est field in metadata.
# Inputs:
- chunk_list_timeseries
# Optional Parametes:
- remove_rv_est (true)
- oversample_factor (1)
- smooth_factor (1)
- alg (:TemporalGP)
# Outputs
- matrix:  timeseries interpolated to each wavelength in the form of a SpectralTimeSeriesCommonWavelengths
- f_mean:  mean flux
- var_mean: variance of mean flux (WARNING: Not computed accurately yet)
- deriv:  dflux/dlnλ averaged over spectra
- deriv2: d²flux/dlnλ² averaged over spectra
"""
function make_template_spectra( chunk_list_timeseries::ACLT; remove_rv_est::Bool = true, oversample_factor::Real = 1, smooth_factor::Real = 1, alg::Symbol = :TemporalGP ) where { ACLT<:AbstractChunkListTimeseries }
    if remove_rv_est
        @assert all(map(md->haskey(md,:rv_est), chunk_list_timeseries.metadata))
    end
    chunk_λ_grids = map(c->RvSpectML.make_grid_for_chunk(chunk_list_timeseries,c,oversample_factor=oversample_factor, remove_rv_est=remove_rv_est), 1:num_chunks(chunk_list_timeseries) )
    local matrix, f_mean, var_mean, deriv, deriv2
    if remove_rv_est
        if alg == :TemporalGP
            ( matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.pack_shifted_chunk_list_timeseries_to_matrix(chunk_list_timeseries,chunk_λ_grids, alg=alg, smooth_factor=smooth_factor )
        else
            if smooth_factor != one(smooth_factor)   @warn "Interpolation algorithm selected does not currently support smoothing."   end
            ( matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.pack_shifted_chunk_list_timeseries_to_matrix(chunk_list_timeseries,chunk_λ_grids, alg=alg, smooth_factor=smooth_factor )
        end
    else
        if alg == :TemporalGP
            ( matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.pack_chunk_list_timeseries_to_matrix(chunk_list_timeseries,chunk_λ_grids, alg=alg, smooth_factor=smooth_factor )
        else
            if smooth_factor != one(smooth_factor)   @warn "Interpolation algorithm selected does not currently support smoothing."   end
            ( matrix, f_mean, var_mean, deriv, deriv2 )  = RvSpectML.pack_chunk_list_timeseries_to_matrix(chunk_list_timeseries,chunk_λ_grids, alg=alg, smooth_factor=smooth_factor )
        end
    end
    return ( matrix, f_mean, var_mean, deriv, deriv2, chunk_λ_grids )
end


text:
using Documenter, FstFileFormat

makedocs(modules=[FstFileFormat],
        doctest=true,
	    sitename = "FstFileFormat.jl",
	    pages = [
		    "The FstFileFormat.jl package" => "index.md"])

deploydocs(
    repo = "github.com/xiaodaigh/FstFileFormat.jl.git"
)


text:
# ------------------------------------------------------------------------------------------
# # Julia for Data Science
#
# * Data
# * **Data processing**
# * Visualization
#
# ### Data processing: Standard machine learning algorithms in Julia
# In what's next, we will see how to use some of the standard machine learning algorithms
# implemented in Julia.
# ------------------------------------------------------------------------------------------

using DataFrames, Statistics

# ------------------------------------------------------------------------------------------
# ### Example 1: Kmeans Clustering
#
# Let's start with some data.
#
# The Sacramento real estate transactions file that we download next is a list of 985 real
# estate transactions in the Sacramento area reported over a five-day period,
# ------------------------------------------------------------------------------------------

download("http://samplecsvs.s3.amazonaws.com/Sacramentorealestatetransactions.csv","houses.csv")
houses = readtable("houses.csv")

# ------------------------------------------------------------------------------------------
# Let's use [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) to plot with the `pyplot`
# backend. (NOTE: this can take a long time the first time you run it, when it's
# initializing the package.)
# ------------------------------------------------------------------------------------------

using Plots
pyplot()
plot(size=(500,500),leg=false)

# ------------------------------------------------------------------------------------------
# Now let's create a scatter plot to show the price of a house vs. its square footage,
# ------------------------------------------------------------------------------------------

x = houses[!, :sq__ft]
# x = houses[7] # equivalent, useful if file has no header
y = houses[!, :price]
# y = houses[10] # equivalent
scatter(x,y,markersize=3)

# ------------------------------------------------------------------------------------------
# *Houses with 0 square feet that cost money?*
#
# The square footage seems to not have been recorded in these cases.
#
# Filtering these houses out is easy to do!
# ------------------------------------------------------------------------------------------

filter_houses = houses[houses[!, :sq__ft] .> 0, :]  # dot broadcasting
x = filter_houses[!, :sq__ft]
y = filter_houses[!, :price]
scatter(x,y)

# ------------------------------------------------------------------------------------------
# This makes sense! The higher the square footage, the higher the price.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can filter a `DataFrame` by feature value too, using the `by` function.
# The `mean()` function comes from the `Statistics` module in the standard library, which we
# get from running `using Statistics` at the top of this file.
# ------------------------------------------------------------------------------------------

by(filter_houses,:type,size)

by(filter_houses,:type,filter_houses->mean(filter_houses[!, :price]))

# ------------------------------------------------------------------------------------------
# Now let's do some kmeans clustering on this data.
#
# First, we can load the `Clustering` package to do this.
# ------------------------------------------------------------------------------------------

#Pkg.add("Clustering")
using Clustering

# ------------------------------------------------------------------------------------------
# Let us see how `Clustering` works with a generic example first.
# ------------------------------------------------------------------------------------------

# make a random dataset with 1000 points
# each point is a 5-dimensional vector
J = rand(5, 1000)
R = kmeans(J, 20; maxiter=200, display=:iter)
# performs K-means over X, trying to group them into 20 clusters
# set maximum number of iterations to 200
# set display to :iter, so it shows progressive info at each iteration

# ------------------------------------------------------------------------------------------
# Now, let's get back to the problem in hand and see how this can be applied over there.
#
# Let's store the features `:latitude` and `:longitude` in an array `X` that we will pass to
# `kmeans`.
#
# First we add data for `:latitude` and `:longitude` to a new `DataFrame` called `X`.
# ------------------------------------------------------------------------------------------

X = filter_houses[!, [:latitude,:longitude]]

# ------------------------------------------------------------------------------------------
# and then we convert `X` to an `Array` via
#
# ```julia
# X = Array(X)
# ```
# or
# ```julia
# X = convert(Array, X)
# ```
#
# Since we know this array has no missing values, we can also change the output type of the
# array to just Float64s, which we'll need for Clustering below:
#
# ```julia
# X = Array{Float64}(X)
# ```
# or
# ```julia
# X = convert(Array{Float64}, X)
# ```
# to turn `X` into an `Array` that stores `Float64`s.
# ------------------------------------------------------------------------------------------

X = Array{Float64}(X)

# ------------------------------------------------------------------------------------------
# We now take the transpose of `X` using the `transpose()` function. A transpose is required
# since `kmeans()` function takes each row as a `feature`, and each column a `data point`.
# ------------------------------------------------------------------------------------------

X = transpose(X)
#X = X'  # (conjugate transposition) also does the same thing (but only for real-valued arrays).
X

# ------------------------------------------------------------------------------------------
# As a first pass at guessing how many clusters we might need, let's use the number of zip
# codes in our data.
#
# (Try changing this to see how it impacts results!)
# ------------------------------------------------------------------------------------------

k = length(unique(filter_houses[!, :zip]))
# there should be atleast 2 distinct features (k>=2) to group the data points
println("unique zip codes are ",k)

# ------------------------------------------------------------------------------------------
# We can use the `kmeans` function to do kmeans clustering!
# ------------------------------------------------------------------------------------------

using Clustering
C = kmeans(X, k) # try changing k

# ------------------------------------------------------------------------------------------
# Now let's create a new data frame, `df`, with all the same data as `filter_houses` that
# also includes a column for the cluster to which each house has been assigned.
# ------------------------------------------------------------------------------------------

df = DataFrame(cluster=C.assignments, city=filter_houses[!, :city],
    latitude=filter_houses[!, :latitude], longitude=filter_houses[!, :longitude], zip=filter_houses[!, :zip])

# ------------------------------------------------------------------------------------------
# Let's plot each cluster as a different color.
# ------------------------------------------------------------------------------------------

clusters_figure = plot()
for i = 1:k
    clustered_houses = df[df[!, :cluster].== i,:]
    xvals = clustered_houses[!, :latitude]
    yvals = clustered_houses[!, :longitude]
    scatter!(clusters_figure,xvals,yvals,markersize=4)
end
xlabel!("Latitude")
ylabel!("Longitude")
title!("Houses color-coded by cluster")
display(clusters_figure)

# ------------------------------------------------------------------------------------------
# And now let's try coloring them by zip code.
# ------------------------------------------------------------------------------------------

unique_zips = unique(filter_houses[!, :zip])
zips_figure = plot()
for uzip in unique_zips
    subs = filter_houses[filter_houses[!, :zip].==uzip,:]
    x = subs[!, :latitude]
    y = subs[!, :longitude]
    scatter!(zips_figure,x,y)
end
xlabel!("Latitude")
ylabel!("Longitude")
title!("Houses color-coded by zip code")
display(zips_figure)

# ------------------------------------------------------------------------------------------
# Let's see the two plots side by side.
# ------------------------------------------------------------------------------------------

plot(clusters_figure,zips_figure,layout=(2, 1))

# ------------------------------------------------------------------------------------------
# Not exactly! but almost... Now we know that ZIP codes are not randomly assigned!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Example 2: Nearest Neighbor with a KDTree
#
# For this example, let's start by loading the `NearestNeighbors` package.
# ------------------------------------------------------------------------------------------

using NearestNeighbors

# ------------------------------------------------------------------------------------------
# With this package, we'll look for the `knearest` neighbors of one of the houses, `point`.
# ------------------------------------------------------------------------------------------

knearest = 10
id = 70 # try changing this
point = X[:,id]

# ------------------------------------------------------------------------------------------
# Now we can build a `KDTree` and use `knn` to look for `point`'s nearest neighbors!
# ------------------------------------------------------------------------------------------

kdtree = KDTree(X)
idxs, dists = knn(kdtree, point, knearest, true)

# ------------------------------------------------------------------------------------------
# We'll first generate a plot with all of the houses in the same color,
# ------------------------------------------------------------------------------------------

x = filter_houses[!, :latitude];
y = filter_houses[!, :longitude];
scatter(x,y);

# ------------------------------------------------------------------------------------------
# and then overlay the data corresponding to the nearest neighbors of `point` in a different
# color.
# ------------------------------------------------------------------------------------------

x = filter_houses[idxs,:latitude];
y = filter_houses[idxs,:longitude];
scatter!(x,y)

# ------------------------------------------------------------------------------------------
# There are those nearest neighbors in red!
#
# We can see the cities of the neighboring houses by using the indices, `idxs`, and the
# feature, `:city`, to index into the `DataFrame` `filter_houses`.
# ------------------------------------------------------------------------------------------

cities = filter_houses[idxs,:city]

# ------------------------------------------------------------------------------------------
# ### Example 3: PCA for dimensionality reduction
#
# Let us try to reduce the dimensions of the price/area data from the houses dataset.
#
# We can start by grabbing the square footage and prices of the houses and storing them in
# an `Array`.
# ------------------------------------------------------------------------------------------

F = filter_houses[!, [:sq__ft,:price]]
F = convert(Array{Float64,2},F)'

# ------------------------------------------------------------------------------------------
# Recall how the data looks when we plot housing prices against square footage.
# ------------------------------------------------------------------------------------------

scatter(F[1,:],F[2,:])
xlabel!("Square footage")
ylabel!("Housing prices")

# ------------------------------------------------------------------------------------------
# We can use the `MultivariateStats` package to run PCA
# ------------------------------------------------------------------------------------------

# Pkg.add("MultivariateStats")
using MultivariateStats

# ------------------------------------------------------------------------------------------
# Use `fit` to fit the model
# ------------------------------------------------------------------------------------------

M = fit(PCA, F)

# ------------------------------------------------------------------------------------------
# Note that you can choose the maximum dimension of the new space by setting `maxoutdim`,
# and you can change the method to, for example, `:svd` with the following syntax.
#
# ```julia
# fit(PCA, F; maxoutdim = 1,method=:svd)
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# It seems like we only get one dimension with PCA! Let's use `transform` to map all of our
# 2D data in `F` to `1D` data with our model, `M`.
# ------------------------------------------------------------------------------------------

y = transform(M, F)

# ------------------------------------------------------------------------------------------
# Let's use `reconstruct` to put our now 1D data, `y`, in a form that we can easily overlay
# (`Xr`) with our 2D data in `F` along the principle direction/component.
# ------------------------------------------------------------------------------------------

Xr = reconstruct(M, y)

# ------------------------------------------------------------------------------------------
# And now we create that overlay, where we can see points along the principle component in
# red.
#
# (Each blue point maps uniquely to some red point!)
# ------------------------------------------------------------------------------------------

scatter(F[1,:],F[2,:])
scatter!(Xr[1,:],Xr[2,:])

# ------------------------------------------------------------------------------------------
# ### Example 4: Learn how to build a simple multi-layer-perceptron on the MNIST dataset
#
# MNIST from: https://github.com/FluxML/model-zoo/blob/master/mnist/mlp.jl
#
# Let's start by loading `Flux`, importing a few things from `Flux` explicitly, and bringing
# the `repeated` function into our scope.
# ------------------------------------------------------------------------------------------

using Flux, Flux.Data.MNIST
using Flux: onehotbatch, argmax, crossentropy, throttle
using Base.Iterators: repeated

# ------------------------------------------------------------------------------------------
# We can now store all the MNIST images in `imgs` and take a peak into this vector to see
# what the data looks like
# ------------------------------------------------------------------------------------------

imgs = MNIST.images()
imgs[3]

# ------------------------------------------------------------------------------------------
# Let's look at the type of an individual image.
# ------------------------------------------------------------------------------------------

typeof(imgs[3])

# ------------------------------------------------------------------------------------------
# #### Reorganizing our array of images
#
# We see this is a 2D array that stores `ColorTypes`. To work more easily with this data,
# let's convert all `ColorTypes` to floating point numbers.
# ------------------------------------------------------------------------------------------

fpt_imgs = float.(imgs)

# ------------------------------------------------------------------------------------------
# Now we can see what `imgs[3]` looks like as an array of floats, rather than as an array of
# colors!
# ------------------------------------------------------------------------------------------

fpt_imgs[3]

# ------------------------------------------------------------------------------------------
# **Let's stack the images to create one large 2D array, `X`, that stores the data for each
# image as a column.**
#
# To do this, we can **first** use `reshape` to unravel each image, creating a 1D array
# (`Vector`) of floats from a 2D array (`Matrix`) of floats.
# ------------------------------------------------------------------------------------------

unraveled_fpt_imgs = reshape.(fpt_imgs, :);
typeof(unraveled_fpt_imgs)

# ------------------------------------------------------------------------------------------
# (Note that `Vector` is an alias for a 1D `Array`.)
# ------------------------------------------------------------------------------------------

Vector

# ------------------------------------------------------------------------------------------
# This makes `unraveled_fpt_imgs` a `Vector` of `Vector`s where `imgs[3]` is now
# ------------------------------------------------------------------------------------------

unraveled_fpt_imgs[3]

# ------------------------------------------------------------------------------------------
# After using `reshape` to get a `Vector` of `Vector`s, we can use `hcat` to build a
# `Matrix`, `X`, from `unraveled_fpt_imgs` where the `Vector`s stored in
# `unraveled_fpt_imgs` will become the columns of `X`.
#
# Note that we're using the "splat" command below, `...`, which allows you to pass all the
# elements of an object to a function, rather than just passing the object itself.
# ------------------------------------------------------------------------------------------

X = hcat(unraveled_fpt_imgs...)

# ------------------------------------------------------------------------------------------
# #### How to go back to images from this 2D `Array`
#
# So now each column in X is an image reshaped to a vector of floating points. Let's pick
# one column and see what the digit is.
#
# Let's try to view the second image in the original array, `imgs`, by taking the second
# column of `X`
# ------------------------------------------------------------------------------------------

onefigure = X[:,2]

# ------------------------------------------------------------------------------------------
# We'll `reshape` this array to a 2D, 28x28 array,
# ------------------------------------------------------------------------------------------

t1 = reshape(onefigure,28,28)

# ------------------------------------------------------------------------------------------
# and finally use `colorview` from the `Images` package to view the handwritten digit.
# ------------------------------------------------------------------------------------------

using Images

colorview(Gray, t1)

# ------------------------------------------------------------------------------------------
# *Our data is in working order!*
#
# For our machine to learn the digit with which each image is associated, we'll need to
# train it using correct answers. Therefore we'll make use of the `labels` associated with
# these images from MNIST.
# ------------------------------------------------------------------------------------------

labels = MNIST.labels() # the true labels

# ------------------------------------------------------------------------------------------
# One-hot-encode the labels with `onehotbatch`
# ------------------------------------------------------------------------------------------

Y = onehotbatch(labels, 0:9)

# ------------------------------------------------------------------------------------------
# which gives a binary indicator vector for each figure
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Build the network
# ------------------------------------------------------------------------------------------

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax)

# ------------------------------------------------------------------------------------------
# Define the loss functions and accuracy
# ------------------------------------------------------------------------------------------

loss(x, y) = Flux.crossentropy(m(x), y)
accuracy(x, y) = mean(argmax(m(x)) .== argmax(y))

methodswith(typeof(ps))

# ------------------------------------------------------------------------------------------
# Use `X` to create our training data and then declare our evaluation function:
# ------------------------------------------------------------------------------------------

datasetx = repeated((X, Y), 200)
evalcb = () -> @show(loss(X, Y))
ps = Flux.params(m)

# ------------------------------------------------------------------------------------------
# So far, we have defined our training data and our evaluation functions.
#
# Let's take a look at the function signature of Flux.train!
# ------------------------------------------------------------------------------------------

?Flux.train!

# ------------------------------------------------------------------------------------------
# **Now we can train our model and look at the accuracy thereafter.**
# ------------------------------------------------------------------------------------------

opt = ADAM()
Flux.train!(loss, ps, datasetx, opt, cb = throttle(evalcb, 10))

accuracy(X, Y)

# ------------------------------------------------------------------------------------------
# Now that we've trained our model, let's create test data, `tX`,
# ------------------------------------------------------------------------------------------

tX = hcat(float.(reshape.(MNIST.images(:test), :))...)

# ------------------------------------------------------------------------------------------
# and run our model on one of the images from `tX`
# ------------------------------------------------------------------------------------------

test_image = m(tX[:,1])

argmax(test_image) - 1

# ------------------------------------------------------------------------------------------
# The largest element of `test_image` is the 8th element, so our model says that test_image
# is a "7".
#
# Now we can look at the original image.
# ------------------------------------------------------------------------------------------

using Images
t1 = reshape(tX[:,1],28,28)
colorview(Gray, t1)

# ------------------------------------------------------------------------------------------
# and there we have it!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Example 5: Linear regression in Julia (we will write our own Julia code and Python
# code)
#
# Let's try to find the best line fit of the following data:
# ------------------------------------------------------------------------------------------

xvals = repeat(1:0.5:10, inner=2)
yvals = 3 .+ xvals .+ 2 .* rand(length(xvals)) .-1
scatter(xvals, yvals, color=:black, leg=false)

# ------------------------------------------------------------------------------------------
# We want to fit a line through this data.
#
# Let's write a Julia function to do this.
# ------------------------------------------------------------------------------------------

function find_best_fit(xvals,yvals)
    meanx = mean(xvals)
    meany = mean(yvals)
    stdx = std(xvals)
    stdy = std(yvals)
    r = cor(xvals,yvals)
    a = r*stdy/stdx
    b = meany - a*meanx
    return a,b
end

# ------------------------------------------------------------------------------------------
# To fit the line, we just need to find the slope and the y-intercept (a and b).
#
# Then add this fit to the existing plot!
# ------------------------------------------------------------------------------------------

a,b = find_best_fit(xvals,yvals)
ynew = a .* xvals .+ b

plot!(xvals,ynew)

# ------------------------------------------------------------------------------------------
# Let's generate a much bigger dataset,
# ------------------------------------------------------------------------------------------

xvals = 1:100000;
xvals = repeat(xvals,inner=3);
yvals = 3 .+ xvals .+ 2 .* rand(length(xvals)) .- 1;

@show size(xvals)
@show size(yvals)

# ------------------------------------------------------------------------------------------
# and now we can time how long it takes to find a fit to this data.
# ------------------------------------------------------------------------------------------

@time a,b = find_best_fit(xvals,yvals)

# ------------------------------------------------------------------------------------------
# Now we will write the same code using Python
# ------------------------------------------------------------------------------------------

using PyCall
using Conda

py"""
import numpy
def find_best_fit_python(xvals,yvals):
    meanx = numpy.mean(xvals)
    meany = numpy.mean(yvals)
    stdx = numpy.std(xvals)
    stdy = numpy.std(yvals)
    r = numpy.corrcoef(xvals,yvals)[0][1]
    a = r*stdy/stdx
    b = meany - a*meanx
    return a,b
"""

find_best_fit_python = py"find_best_fit_python"

xpy = PyObject(xvals)
ypy = PyObject(yvals)
@time a,b = find_best_fit_python(xpy,ypy)

# ------------------------------------------------------------------------------------------
# **Let's use the benchmarking package to time these two.**
# ------------------------------------------------------------------------------------------

using BenchmarkTools

@btime a,b = find_best_fit_python(xvals,yvals)

@btime a,b = find_best_fit(xvals,yvals)

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-julia-for-data-science/2. Julia for Data Science - Algorithms.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials