text:
"""
    daughter = mother(this::CWT{W, T, <:ContWaveClass, N},
                        s::Real, nInOctave::Int, ω::AbstractArray{<:Real,1}) where {W, T, N}

given a CWT object, return a rescaled version of the mother wavelet, in the
fourier domain. ω is the frequency, which is fftshift-ed. s is the scale
variable.
"""
function mother(this::CWT{W,T,Morlet,N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W,T,N}
    constant = this.σ[3] * (π)^(1 / 4)
    gauss = exp.(-(this.σ[1] .- ω / s).^2 / (sWidth^2))
    shift = this.σ[2] * exp.(-1 / 2 * (ω / s).^2)
    daughter = constant .* (gauss .- shift)
    return normalize(daughter, s, this.p)
end

function mother(this::CWT{W,T,<:Paul,N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W,T,N}
    daughter = zeros(length(ω))
    constant = (2^this.α) / sqrt((this.α) * gamma(2 * (this.α)))
    polynomial = (ω[ω .>= 0] / s).^(this.α)
    expDecay = exp.(-(ω[ω .>= 0] / s))
    daughter[ω .>= 0] = constant .* polynomial .* expDecay
    return normalize(daughter, s, this.p)
end

function mother(this::CWT{W,T,<:Dog,N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W,T,N}
    constant = im^(this.α) / sqrt(gamma((this.α) + 1 / 2))
    polynomial = (ω / s).^(this.α)
    gauss = exp.(-(ω / s).^2 / 2)
    daughter =  constant .* polynomial .* gauss
    return normalize(daughter, s, this.p)
end

function mother(this::CWT{W,T,Morse,N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W,T,N}

    ga = this.waveType.ga;
    be = this.waveType.be;
    cf = this.waveType.cf;
    p = this.p;

    fo = morsefreq(this) * 2 * pi;
    fact = cf/fo;

    om = (ω / s);

    if be == 0
        daughter = 2 * exp.(-om.^ga);
    else
        daughter = 2 * exp.(-be.*log.(fo) .+ fo.^ga .+ be.*log.(om) .- om.^ga);
    end

    daughter[1] = 1/2*daughter[1]; # Due to unit step function
    # Ensure nice lowpass filters for beta=0;
    # Otherwise, doesn't matter since wavelets vanishes at zero frequency

    if any(daughter .!= daughter)
        @warn "the given values of gamma and beta are numerically unstable"
        daughter[daughter .!= daughter] .= 0;
    end

    return ContinuousWavelets.normalize(daughter, s, p)
end

function mother(this::CWT{W,T,<:ContOrtho,N}, s::Real, itpψ,
                ω::AbstractArray{<:Real,1}, n, n1) where {W,T,N}
    daughter = itpψ(range(1, stop=length(itpψ), step=length(itpψ) / n1 * s))
    daughter = padTo(daughter, n)
    daughter = circshift(daughter, -round(Int, n1 / s / 2))
    return daughter
end
function normalize(daughter, s, p)
    if p == Inf
        normTerm = maximum(abs.(daughter))
    else
        normTerm = s^(1 / p)
    end
    return daughter ./ normTerm
end

"""
    father(c::CWT, ω, averagingType::aT) where {aT<:Averaging}

this creates the averaging function, which covers the low frequency
information, and is emphatically not analytic. aT determines whether it has the
same form as the wavelets (`Father`), or just a bandpass `Dirac`.
`c.averagingLength` gives the number of octaves (base 2) that are covered by
the averaging function. The width is then derived so that it matches the next
wavelet at 1σ.

For the Morlet wavelet, the distribution is just a Gaussian, so it has variance
1/s^2 and mean σ[1]*s set the variance so that the averaging function has σ/2 at
the central frequency of the last scale.

For the Paul wavelets, it's a easy calculation to see that the mean of a paul
wavelet of order m is (m+1)/s, while σ=sqrt(m+1)/s. So we set the variance so
that the averaging function has 1σ at the central frequency of the last scale.

the derivative of a Gaussian (Dog) has a pretty nasty form for the mean and
variance; eventually, if you set 1/2 * σ_{averaging}=⟨ω⟩_{highest scale wavelet}, you
will get the scale of the averaging function to be
`s*gamma((c.α+2)/2)/sqrt(gamma((c.α+1)/2)*(gamma((c.α+3)/2)-gamma((c.α+2)/2)))`

"""
function father(c::CWT, ω, averagingType::Father, sWidth)
    s = 2^(min(c.averagingLength + getMinScaling(c) - 1, 0))
    s0, ω_shift = locationShift(c, s, ω, sWidth)
    averaging = adjust(c) .* mother(c, s0, 1, ω_shift)
end

function father(c::CWT{B,T,W}, ω, averagingType::Father,
                       fullVersion, s, N, n1) where {W <: ContOrtho,B,T}
    itp = genInterp(fullVersion)

    φ = itp(range(1, stop=length(fullVersion),
                  step=length(fullVersion) / n1 * s))
    φ = circshift(padTo(φ, N), -round(Int, n1 / s / 2))
end

# dirac version (that is, just a window around zero)
function father(c::CWT{<:WaveletBoundary,T}, ω,
                       averagingType::Dirac, sWidth) where {T}
    s = 2^(c.averagingLength) * sWidth
    averaging = zeros(T, size(ω))
    upperBound = getUpperBound(c, s)
    averaging[abs.(ω) .<= upperBound] .= 1
    return averaging
end

function father(c::CWT{W, T, <:Morse}, ω, averagingType::ContinuousWavelets.Father, sWidth) where {W,T}
    s = 2^(min(c.averagingLength + getMinScaling(c) - 1, 0))
    s0, ω_shift = locationShift(c, s, ω, sWidth)
    averaging = adjust(c) .* mother(c, s0, sWidth, ω_shift)
    return averaging
end

@doc """
      (daughters, ω) = computeWavelets(n1::Integer, c::CWT{W}; T=Float64, J1::Int64=-1, dt::S=NaN, s0::V=NaN) where {S<:Real,
                                                                   W<:WaveletBoundary, V}
just precomputes the wavelets used by transform c::CWT{W}. For details, see cwt
"""
function computeWavelets(n1::Integer, c::CWT{B,CT,W}; T=Float64, J1::Int64=-1, dt::S=NaN, s0::V=NaN, space=false) where {S <: Real,B <: WaveletBoundary,V,W,CT}
    # don't alter scaling with sampling information if it doesn't exists
    fλ = (4 * π) / (c.σ[1] + sqrt(2 + c.σ[1]^2))
    if isnan(dt) || (dt < 0)
        dt = 1
    end
    # smallest resolvable scale
    if isnan(s0) || (s0 < 0)
        s0 = 2 * dt / fλ
    end

    nOctaves, totalWavelets, sRange, sWidth = getNWavelets(n1, c)
    # padding determines the actual number of elements
    n, nSpace = setn(n1, c)
    # indicates whether we should keep a spot for the father wavelet
    isAve = !(typeof(c.averagingType) <: NoAve)
    # I guess matlab did occasionally do something useful

    ω = computeOmega(n1, nSpace, n)
    daughters = analyticOrNot(c, n, totalWavelets)

    # if the nOctaves is small enough there are none not covered by the
    # averaging, just use that
    if totalWavelets == 1
        @warn "an averaging length of $(c.averagingLength) results in fewer octaves than the minimum frequency $(getMinScaling(c)) for this wavelet type. Either decreasing the averagingLength or increase extraOctaves"
        onlyFather = father(c, ω, c.averagingType, sWidth[1])
        return onlyFather, ω
    end

    for (curWave, s) in enumerate(sRange)
        daughters[:,curWave + isAve] = mother(c, s, sWidth[curWave], ω)
    end
    if isAve # should we include the father?
        @debug "c = $(c), $(c.averagingType), size(ω)= $(size(ω))"
        daughters[:, 1] = father(c, ω, c.averagingType, sWidth[1])# [1:(n1+1)]
    end

    # adjust by the frame bound
    if c.frameBound > 0
        daughters = daughters .* (c.frameBound / norm(daughters, 2))
    end
    testFourierDomainProperties(daughters, isAve)
    if space
        x = zeros(T, n1); x[ceil(Int, n1 / 2)] = 1
        return cwt(x, c, daughters), ω
    else
        return (daughters, ω)
    end
end

function computeWavelets(n1::Integer, c::CWT{B,CT,W}; T=Float64, space=false) where {S <: Real,B <: WaveletBoundary,V,W <: ContOrtho,CT}

    # padding determines the actual number of elements
    n, nSpace = setn(n1, c)
    nOctaves, totalWavelets, sRange, sWidth = getNWavelets(n1, c)
    # indicates whether we should keep a spot for the father wavelet
    isAve = !(typeof(c.averagingType) <: NoAve)
    # I guess matlab did occasionally do something useful

    ω = computeOmega(n1, nSpace, n)
    daughters = zeros(nSpace, totalWavelets)
    φ, ψ, ψLen = getContWaveFromOrtho(c, nSpace)
    itpψ = genInterp(ψ)
    # if the nOctaves is small enough there are none not covered by the
    # averaging, so just use that
    if totalWavelets == 1
        onlyFather = father(c, ω, c.averagingType, φ, c.averagingLength, nSpace, n1)
        return rfft(onlyFather, 1), ω
    end

    for (curWave, s) in enumerate(sRange)
        daughters[:,curWave + isAve] = mother(c, s, itpψ, ω, nSpace, n1)
    end
    if isAve
        daughters[:, 1] = father(c, ω, c.averagingType, φ, sRange[1],
                                        nSpace, n1)
    end
    # switch to fourier domain and normalize appropriately
    daughters = rfft(daughters, 1)
    repeatSRange = [sRange[1:isAve]...; sRange...]
    for i in 1:size(daughters, 2)
        daughters[:, i] = normalize(daughters[:, i],
                                            repeatSRange[i],
                                            c.p)
    end

    # adjust by the frame bound
    if c.frameBound > 0
        daughters = daughters .* (c.frameBound / norm(daughters, 2))
    end
    testFourierDomainProperties(daughters, isAve)
    if space
        x = zeros(T, n1); x[ceil(Int, n1 / 2)] = 1
        return cwt(x, c, daughters), ω
    else
        return (daughters, ω)
    end
end

# not
function analyticOrNot(c::CWT{W,T,<:Union{Dog,ContOrtho},N}, n, totalWavelets) where {W,T,N}
    # Odd derivatives (and any orthogonal wavelets) introduce an imaginary
    # term, so we need a complex representation
    if abs(c.α % 2) == 1
        daughters = zeros(Complex{T}, n, totalWavelets)
    else
        daughters = zeros(T, n, totalWavelets)
    end
    return daughters
end

function analyticOrNot(c::CWT{W,T,<:Union{Morlet,Paul,Morse},N}, n, totalWavelets) where {W,T,N}
    daughters = zeros(T, n, totalWavelets)
    return daughters
end

function computeOmega(nOriginal, nSpace, nFreq)
    range(0, nOriginal >> 1 + 1, length=nFreq) # max size is the last frequency in the rfft of the original data size
end
# convenience methods

# it's ok to just hand the total size, even if we're not transforming across
# all dimensions
function computeWavelets(Y::Tuple, c::CWT{W}; T=Float64) where {S <: Real,W <: WaveletBoundary}
    return computeWavelets(Y[1], c; T=T)
end
function computeWavelets(Y::AbstractArray{<:Integer}, c::CWT{W}; T=Float64) where {S <: Real,W <: WaveletBoundary}
    return computeWavelets(Y[1], c, T=T)
end

# also ok to just hand the whole thing being transformed
function computeWavelets(Y::AbstractArray{<:Number}, c::CWT{W}; T=Float64) where {S <: Real,W <: WaveletBoundary}
    return computeWavelets(size(Y)[1], c, T=T)
end


text:
using OpenQuantumBase, Random

bath = EnsembleFluctuator([1, 1], [1, 2])
interaction = InteractionSet(Interaction(ConstantCouplings(["Z"], unit=:ħ), bath))

fluct = OpenQuantumBase.fluctuator_from_interactions(interaction)[1]
cache_exp = -1.0im * sum(fluct.b0 .* [σz, σz])

cache = zeros(ComplexF64, 2, 2)
p = ODEParams(nothing, 5.0, (tf, t) -> t / tf)
update_cache!(cache, fluct, p, 2.5)
@test cache ≈ cache_exp

cache = zeros(ComplexF64, 4, 4)
update_vectorized_cache!(cache, fluct, p, 2.5)
@test cache ≈ one(cache_exp)⊗cache_exp - transpose(cache_exp)⊗one(cache_exp)

Random.seed!(1234)
random_1 = rand([-1, 1], 2, 1)
τ, idx = findmin(rand(fluct.dist, 1))
Random.seed!(1234)
OpenQuantumBase.reset!(fluct, (x, y) -> rand([-1, 1], x, y))

random_1[idx] *= -1
@test fluct.b0[1] == random_1[1] && fluct.b0[2] == random_1[2]
@test fluct.next_τ ≈ τ
@test fluct.next_idx == idx

b0 = copy(fluct.b0)
Random.seed!(1234)
random_1 = rand(fluct.dist, 1)
τ, idx = findmin(random_1)
b0[idx] *= -1
Random.seed!(1234)
OpenQuantumBase.next_state!(fluct)
@test fluct.b0[1] == b0[1] && fluct.b0[2] == b0[2]
@test fluct.next_τ ≈ τ
@test fluct.next_idx == idx

ext: jl


text:
abstract F{f}
abstract C{c}

type Cursor{S, Co<:Tuple, Fi<:Tuple} <: IO
    ncols::Int
    fields::Type{Fi}
    coltypes::Type{Co}
    io::S
    status::Int
    row::Int
end

# primitives

ncols(curs::Cursor) = curs.ncols
Base.read(curs::Cursor) = nothing
step!(curs::Cursor) = nothing
reset!(curs::Cursor) = nothing
Base.close(curs::Cursor) = nothing

# out utilities

function materialize!(curs::Cursor, restype::DataType)
    Co = typeof(curs).parameters[2]
    Fi = typeof(curs).parameters[3]
    n = ncols(curs)
    columns, nrows = _materialize(curs, restype)
    return columns, Fi, nrows
end

# show

function Base.summary(curs::Cursor)
    fieldnames = gen_fieldnames(collect(curs.fields.parameters))
    res = "Cursor with columns $fieldnames"
end

function gen_fieldnames(fields)
    # fieldnames = fields.parameters
    n = length(fields)
    res = "("
    for i in eachindex(fields)
        i == n ? res = res * "$(fields[i]))" :
                 res = res * "$(fields[i]), "
    end
    return res
end


text:
#! basic implementation

"""
basic cobb douglas
"""
function cobbdouglas(x, w)
    z = one(eltype(x))
    for i in eachindex(x)
        z *= x[i]^w[i]
    end
    return z
end

"""
basic ces
"""
function _ces(x, w, r)
    z = zero(eltype(x))
    for i in eachindex(x)
        z += w[i] * x[i]^r
    end
    return z^inv(r)
end

"""
constant elasticity of substitution function
"""
function ces(x, w, r)
    n = length(x)
    (n == length(w)) || throw(DimensionMismatch())
    r == zero(eltype(x)) ? cobbdouglas(x, w) : _ces(x, w, r)
end


text:
test_gage_min_distance(
    Array[[1., 0., 0.], [0., 0., 1.]],
    [0., 0., 1.],
    0.,
    "one simple array and one array of arrays"
)


text:
#*****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Experiment 4
test combined time series and binary logical operations
"""
module NNGrammarExpt4

export circuit4, restarts

import Compat.ASCIIString
using TFTools
using Datasets
using TensorFlow
using TensorFlow.CoreTypes #DT_FLOAT32, DT_BOOL, PyVectorType
using TensorFlow.API: l2_loss, AdamOptimizer, GradientDescentOptimizer, cast,
    round_, reshape_, l2_normalize, RMSPropOptimizer,
    reduce_max, reduce_min, minimum_, maximum_, transpose_, less_, greater,
    expand_dims, tile, shape, mul, reduce_sum, moments, conv2d, truncated_normal,
    concat, sub_, matmul, clip_by_value
using TensorFlow.API.TfNn: sigmoid
using StatsBase
using RLESUtils, Confusion

using PyCall
@pyimport tensorflow as tf

function restarts(f::Function, N::Int64; kwargs...)
   [f(; kwargs...) for i = 1:N]
end

function circuit4(;
    dataname::AbstractString="vhdist3",
    labelfield::AbstractString="nmac", #"F_x1_lt_100_and_x2_lt_500",
    learning_rate::Float64=0.001,
    max_training_epochs::Int64=500,
    target_cost::Float64=0.001,
    batch_size::Int64=4820,
    mux_hidden_units::Vector{Int64}=Int64[50,30],
    display_step::Int64=1,
    b_debug::Bool=false,
    nshow::Int64=1)

    Dl = dataset(dataname, :nmac; transform=x->Float32(x)) #DFSetLabeled
    val_inputs = constant(Float32[0, 50, 100, 250, 500, 1000, 3000, 5000])
    data_set = TFDataset(Dl)

    # Construct model
    (n_examples, n_steps, n_feats) = size(Dl)
    n_featsflat = n_steps * n_feats
    n_vals = get_shape(val_inputs)[1]

    # inputs
    feats = Placeholder(DT_FLOAT32, [-1, n_steps, n_feats])
    inputs = Tensor(feats)

    normalizer = Normalizer(data_set; converttype=Float32)
    normed_input = normalize01(normalizer, inputs)
    feats_flat = reshape_(normed_input, constant(Int32[-1, n_featsflat]))

    #softness parameter
    softness = collect(linspace(0.1, 1.0, max_training_epochs))
    softness_pl = Placeholder(DT_FLOAT32, [1])

    #toggle hard or soft output
    harden_pl = Placeholder(DT_BOOL, [1])

    # common (embedding) layer
    #embed_in = constant(rand(Float32, 1, 5))
    #embed_hidden_units = [5]
    #embed_blk = ReluStack(embed_in, embed_hidden_units)
    #embed_out = out(embed_blk)

    # mux select input
    #muxselect = feats_flat #simple flatten
    #muxselect = constant(ones(Float32, 1, 1)) #constant 1
    #muxselect = embed_out #relustack embedding

    overrides = zeros(Int64, 8)
    #overrides = [1, 2, 3, 5, 1, 1, 1, 1]
    #overrides = [0, 2, 0, 5, 0, 1, 1, 0]

    #convolution stuff
    HT = 3
    WT = 2
    STRIDES = PyVectorType([1, 1, 1, 1])
    filt_weights = Variable(truncated_normal(constant([HT, WT, 1, 1]), constant(0.0), constant(5e-2)))
    conv1 = conv2d(expand_dims(Tensor(normed_input), constant(3)), Tensor(filt_weights), STRIDES, "SAME")
    muxselect = reshape_(conv1, constant(Int32[-1, n_featsflat]))

    # f1 feat select
    f1_in = inputs
    f1_mux = SoftMux(n_feats, mux_hidden_units, f1_in, muxselect, Tensor(harden_pl), Tensor(softness_pl); override=overrides[1])
    f1_out = out(f1_mux)

    # f2 feat select
    f2_in = inputs
    f2_mux = SoftMux(n_feats, mux_hidden_units, f2_in, muxselect, Tensor(harden_pl), Tensor(softness_pl); override=overrides[2])
    f2_out = hardout(f2_mux)

    # v1 value select
    #v1_in = val_inputs
    n_vals=1
    v1_in = constant(Float32[100.0])
    v1_mux = SoftMux(n_vals, mux_hidden_units, v1_in, muxselect, Tensor(harden_pl); override=overrides[3])
    v1_out = out(v1_mux)

    # v2 value select
    #v2_in = val_inputs
    n_vals = 1
    v2_in = constant(Float32[500.0])
    v2_mux = SoftMux(n_vals, mux_hidden_units, v2_in, muxselect, Tensor(harden_pl); override=overrides[4])
    v2_out = hardout(v2_mux)

    compare_ops = [op_lt, op_gt]
    logical_ops = [op_and, op_or]
    temporal_ops = [op_F, op_G]

    # a1 float op block
    a1_in = (f1_out, v1_out)
    a1_blk = SoftOpsMux(a1_in, compare_ops, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl); opargs=Any[Tensor(softness_pl)], override=overrides[5])
    a1_out = out(a1_blk)

    # a2 float op block
    a2_in = (f2_out, v2_out)
    a2_blk = SoftOpsMux(a2_in, compare_ops, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl); opargs=Any[Tensor(softness_pl)], override=overrides[6])
    a2_out = hardout(a2_blk)

    # l1 logical op block
    l1_in = (a1_out, a2_out)
    l1_blk = SoftOpsMux(l1_in, logical_ops, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl); override=overrides[7])
    l1_out = hardout(l1_blk)

    # t1 temporal op block
    t1_in = (l1_out,)
    #t1_in = (constant(rand(Float32, shape(l1_out))), )
    t1_blk = SoftOpsMux(t1_in, temporal_ops, mux_hidden_units, muxselect, Tensor(harden_pl), Tensor(softness_pl); override=overrides[8])
    t1_out = out(t1_blk)

    # outputs
    pred = t1_out
    labels = Placeholder(DT_FLOAT32, [-1])

    #take nnout of batch, compute moments, take variance component and sum
    f1_var = reduce_sum(moments(f1_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    v1_var = reduce_sum(moments(v1_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    f2_var = reduce_sum(moments(f2_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    v2_var = reduce_sum(moments(v2_mux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    a1_var = reduce_sum(moments(a1_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    a2_var = reduce_sum(moments(a2_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    l1_var = reduce_sum(moments(l1_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    t1_var = reduce_sum(moments(t1_blk.softmux.nnout, Tensor(Int32[0]))[2], Tensor(0))
    sum_var = f1_var + v1_var + f2_var + v2_var + constant(300.0) .* a1_var + constant(300.0) .* a2_var + constant(300.0) .* l1_var + constant(300.0) .* t1_var

    #packed_nnouts = concat(Tensor(1), Tensor([f1_mux.nnout, v1_mux.nnout, f2_mux.nnout, v2_mux.nnout, a1_blk.softmux.nnout, a2_blk.softmux.nnout, l1_blk.softmux.nnout, t1_blk.softmux.nnout]))
    #m = matmul(packed_nnouts, packed_nnouts, false, true)
    #allpairs_cossim = reduce_sum(sub_(constant(8.0), m))

    # Define loss and optimizer
    #cost = l2_loss(pred - labels) # Squared loss
    cost = l2_loss(pred - labels) + constant(1.0) .* sum_var
    #cost = l2_loss(pred - labels) + constant(1.0) .* allpairs_cossim

    #optimizer
    #optimizer = minimize(AdamOptimizer(learning_rate), cost)
    opt = Optimizer(tf.train[:GradientDescentOptimizer](learning_rate))
    gvs = opt.x[:compute_gradients](cost.x)
    capped_gvs = [(tf.nn[:l2_normalize](tf.clip_by_value(grad, -1.0, 1.0), 0), var) for (grad, var) in gvs]
    optimizer = Operation(opt.x[:apply_gradients](capped_gvs))

    #compiled hardselect
    ckt = Circuit([
        f1_mux, v1_mux,
        f2_mux, v2_mux,
        a1_blk, a2_blk,
        l1_blk,
        t1_blk,
        ],
        Vector{ASCIIString}[
        ["alt", "range"],
        ["0", "50", "100", "250", "500", "1000", "3000", "5000"],
        ["alt", "range"],
        ["0", "50", "100", "250", "500", "1000", "3000", "5000"],
        ["<", ">"],
        ["<", ">"],
        ["and", "or"],
        ["F", "G"]])

    #debug
    f1_grad = gradient_tensor(cost, f1_mux.weight)
    v1_grad = gradient_tensor(cost, v1_mux.weight)
    f2_grad = gradient_tensor(cost, f2_mux.weight)
    v2_grad = gradient_tensor(cost, v2_mux.weight)
    a1_grad = gradient_tensor(cost, a1_blk.softmux.weight)
    a2_grad = gradient_tensor(cost, a2_blk.softmux.weight)
    l1_grad = gradient_tensor(cost, l1_blk.softmux.weight)
    t1_grad = gradient_tensor(cost, t1_blk.softmux.weight)
    #/debug

    # Iniuializing the variables
    init = initialize_all_variables()

    # Rock and roll
    println("Optimization Start: $(now())")
    sess = Session()
    run(sess, init)

    #debug
    #fd = FeedDict(feats => data_set.X, labels => data_set.Y)
    #run(sess, normed_input, fd)
    #/debug

    # Training cycle
    for epoch in 1:max_training_epochs
        avg_cost = 0.0
        total_batch = div(num_examples(data_set), batch_size)

        #Loop over all batches
        for i in 1:total_batch
            batch_xs, batch_ys = next_batch(data_set, batch_size)
            fd = FeedDict(feats => batch_xs, labels => batch_ys, softness_pl => [softness[epoch]], harden_pl => [false])
            # Fit training using batch data
            run(sess, optimizer, fd)

            #debug
            #softsel = softselect_by_example(sess, ckt, fd)
            #grads = run(sess, Tensor([f1_grad, v1_grad, f2_grad, v2_grad, a1_grad, a2_grad, l1_grad, t1_grad]), fd)
            #conv1_out = run(sess, conv1, fd)
            #@show run(sess, pred, fd)
            #@show run(sess, pred - labels, fd)
            #@show run(sess, muxselect, fd)
            #@show run(sess, f1_grad, fd)
            #@show run(sess, v1_grad, fd)
            #/debug

            # Compute average loss
            batch_average_cost = run(sess, cost, fd)
            avg_cost += batch_average_cost / (total_batch * batch_size)
        end

        # Display logs per epoch step
        if epoch % display_step == 0
            #softsel = softselect_by_example(sess, ckt, fd)
            #grads = run(sess, Tensor([f1_grad, v1_grad, f2_grad, v2_grad, a1_grad, a2_grad, l1_grad, t1_grad]), fd)
            #println(softsel)
            #println(capped_gvs)
            #grads = Tensor([grad for (grad, var) in capped_gvs])
            println("Epoch $(epoch)  cost=$(avg_cost)")
            #println("Norm=$(norm(grads))")
            if avg_cost < Float32(target_cost)
                break;
            end
        end
    end
    println("Optimization Finished: $(now())")

    # Test model
    correct_prediction = (round_(pred) == labels)
    # Calculate accuracy
    accuracy = mean(cast(correct_prediction, DT_FLOAT32))

    #reload data_set to recover original order (next_batch will internally shuffle)
    data_set = TFDataset(Dl)
    X = data_set.X
    Y = data_set.Y

    #soft metrics
    fd = FeedDict(feats => data_set.X, labels => data_set.Y, softness_pl => [softness[end]], harden_pl => [false])
    Ypred_soft = run(sess, round_(pred), fd)
    acc_soft = run(sess, accuracy, fd)
    stringout = simplestring(sess, ckt, fd; order=[8,1,5,2,7,3,6,4])
    top5 = topstrings(stringout, 5)
    softsel = softselect_by_example(sess, ckt, fd)
    grads = run(sess, Tensor([f1_grad, v1_grad, f2_grad, v2_grad, a1_grad, a2_grad, l1_grad, t1_grad]), fd)
    conf = confusion(Ypred_soft.==1.0, Y.==1.0)
    conf_indices = confusion_indices(Ypred_soft.==1.0, Y.==1.0)

    #hard metrics
    fd = FeedDict(feats => data_set.X, labels => data_set.Y, softness_pl => [softness[end]], harden_pl => [true])
    Ypred_hard = run(sess, round_(pred), fd)
    acc_hard = run(sess, accuracy, fd)

    println("Soft Accuracy:", acc_soft)
    println("Hard Accuracy:", acc_hard)
    println(top5)
    d=Dict{ASCIIString,Any}()

    top5, acc_hard
end

op_F(x::Tensor) = reduce_max(x, Tensor(1))
op_G(x::Tensor) = reduce_min(x, Tensor(1))
op_and(x::Tensor, y::Tensor) = minimum_(x, y) #element-wise
op_or(x::Tensor, y::Tensor) = maximum_(x, y) #element-wise

#TODO: move these to a central location
function op_gt(x::Tensor, y::Tensor, W::Tensor=constant(1.0/100.0))
    tmp = expand_dims(y, Tensor(1))
    ytiled = tile(tmp, Tensor([1, get_shape(x)[2]]))
    result = sigmoid(mul(W, x - ytiled))
    result
end
function op_lt(x::Tensor, y::Tensor, W::Tensor=constant(1.0/100.0))
    tmp = expand_dims(y, Tensor(1))
    ytiled = tile(tmp, Tensor([1, get_shape(x)[2]]))
    result = sigmoid(mul(W, ytiled - x))
end

end #module


text:
struct MercuryBarycenter <: abstractBarycenter end

naif_id(::Type{MercuryBarycenter}) = 1
parent(::Type{MercuryBarycenter}) = SSB


text:
## Design Decisions
#
# - should we have extra code for where hte number of smaples is 1 or just have
# generic routines which return first element of an array.
#
# - Second majo design question is who what kind of random variable to pass to the solver
# initially we tried to make the preimage smaplers agnostic to what type of random variable
# it was passed.  The randvar need only implement call(X, A).  However
# -- If the algorith is parallelised, it may need generate the C typed rand Varvar
# -- sol 1. Make these preimage samplers expect the RandVar type and they can do their own conversion
# -- sol 2. Have some kind of buffer between the different cases

# How might that work, well we say

## Unconditional Sampling
## ======================
rand{T}(X::ExecutableRandVar{T}) = call(X,LazyRandomVector(Float64))

"Generate `n` unconditioned random samples from distribution of X"
rand{T}(X::ExecutableRandVar{T}, n::Integer) =
  T[call(X, LazyRandomVector(Float64)) for i = 1:n]

rand{T}(X::SymbolicRandVar{T}, n::Integer) =
  rand(convert(ExecutableRandVar{T},X),n)

##  Rand Arrays
"Generate `n` unconditioned random array samples from distribution of X"
rand{T}(X::ExecutableRandArray{T}, n::Integer) =
  Array{T}[call(X, LazyRandomVector(Float64)) for i = 1:n]

"Generate `n` unconditioned random array samples from distribution of X"
rand{T,N}(Xs::RandArray{T,N}, n::Integer) = rand(convert(ExecutableRandArray{T},Xs),n)

## RandVar{Bool} Preimage Samples
## ==============================

# Note args named x_sampler sample *from* x, e.g.
# partition_sampler samples set from partition (not a partition itself)

"`n` abstract samples from preimage: Y^-1({true})"
function abstract_sample_partition(
    Y::SymbolicRandVar{Bool},
    n::Integer;
    partition_alg::Type{BFSPartition} = BFSPartition,
    args...)

  partition = pre_partition(Y, partition_alg; args...)
  rand(preiamge, n)
end

"`n` point Sample from preimage: Y^-1({true})"
function point_sample_partition{T<:PartitionAlgorithm}(
    Y::SymbolicRandVar{Bool},
    n::Integer;
    partition_alg::Type{T} = BFSPartition,
    partition_sampler::Function = point_sample,
    args...)
  # FIXME: Float64 too specific
  p = pre_partition(Y, partition_alg; args...)
  s_p = SampleablePartition(p)
  partition_sampler(s_p, n)
end

## Markokv Chain Conditional Sampling
## ==================================
"`n` approximate point sample from preimage: Y^-1({true})"
function point_sample_mc{T<:MCMCAlgorithm}(
    Y::SymbolicRandVar{Bool},
    n::Integer;
    ChainAlg::Type{T} = AMS,    # Generate Markov Chain of samples
    chain_sampler::Function = point_sample, # Sample from Markov Chain
    args...)

  # FIXME: Float64 too specific
  chain = pre_mc(Y, n, ChainAlg; args...)
  chain_sampler(chain)
end

## Samples from X given Y
## ======================

# FIXME, this is not the best way to do it.  This method, finds the preiamge samples
# Then just runs them using interva arithmetic.  In order to do this properly you
# need to pave both ways

"`n` conditional samples from `X` given `Y` is true"
function rand{T}(
    X::SymbolicRandVar{T},
    Y::SymbolicRandVar{Bool},
    n::Integer;
    preimage_sampler::Function = point_sample_mc,
    args...)

  executable_X = convert_psuedo(ExecutableRandVar{T}, X)
  preimage_samples = preimage_sampler(Y, n; args...)
  T[call(executable_X, sample) for sample in preimage_samples]
end

"When `X` is a rand var"
function rand{T}(
    X::RandArray{T},
    Y::SymbolicRandVar{Bool},
    n::Integer;
    preimage_sampler::Function = point_sample_mc,
    args...)

  executable_X = convert_psuedo(ExecutableRandArray{T}, X)
  preimage_samples = preimage_sampler(Y, n; args...)
  Array{T}[call(executable_X, sample) for sample in preimage_samples]
end

"Sample from a tuple of values `(X_1, X_2, ..., X_m) conditioned on `Y`"
function rand(
    X::Tuple,
    Y::SymbolicRandVar{Bool},
    n::Integer;
    preimage_sampler::Function = point_sample_mc,
    args...)

  preimage_samples = preimage_sampler(Y, n; args...)

  # There are two natural ways to return the tuples
  # 1. tuple of m (num in tuple) vectors, each n samples long
  # 2. vector of `n` tuples of length `m` <-- we do this oen

  # types = map(x->Vector{rangetype(x)}, X)
  samples = Array[]
  for x in X
    RT = rangetype(x)
    executable_X = executionalize(x)
    xsamples = RT[call(executable_X, sample) for sample in preimage_samples]
    push!(samples, xsamples)
  end
  samples
  map(tuple, samples...)
  # tuple(samples...)
end

## One Sample
## ==========

"Generate a sample from a rand array `Xs` conditioned on `Y`"
rand(Xs::RandArray, Y::SymbolicRandVar{Bool}; args...) =  rand(Xs,Y,1;args...)[1]

"Generate a sample from a rand array `Xs` conditioned on `Y`"
rand(Xs::Ex, Y::SymbolicRandVar{Bool}; args...) =  rand(Xs,Y,1;args...)[1]

"Generate a sample from a randvar `X` conditioned on `Y`"
rand(X::SymbolicRandVar, Y::SymbolicRandVar{Bool}; args...) = rand(X,Y,1; args...)[1]

"Generate an unconditioned random sample from X"
rand{T}(X::SymbolicRandVar{T}) = rand(X,1)[1]

"Generate single conditionam sample of tuple `X` of RandVar/Arrays given `Y`"
rand(X::Tuple, Y::SymbolicRandVar{Bool}; args...) = rand(X,Y,1;args...)[1]

"Generate single sample from random array"
rand(X::RandArray; args...) = rand(X, 1)[1]


text:
function testAssembly(diffusion_system::System)
    u = addVariable!(diffusion_system, "u")

    diffusion_kernel = Diffusion(u)

    addKernel!(diffusion_system, diffusion_kernel)

    initialize!(diffusion_system)

    solver = JuliaDenseImplicitSolver(diffusion_system)

    initialize!(solver)

    MOOSE.assembleResidualAndJacobian(solver, diffusion_system)

    # From MOOSE
    jac = [0.666666666667 -0.166666666667 -0.333333333333 -0.166666666667 0 0 0 0 0 ;
           -0.166666666667 1.33333333333 -0.333333333333 -0.333333333333 -0.166666666667 -0.333333333333 0 0 0 ;
           -0.333333333333 -0.333333333333 2.66666666667 -0.333333333333 -0.333333333333 -0.333333333333 -0.333333333333 -0.333333333333 -0.333333333333 ;
           -0.166666666667 -0.333333333333 -0.333333333333 1.33333333333 0 0 -0.333333333333 -0.166666666667 0 ;
           0 -0.166666666667 -0.333333333333 0 0.666666666667 -0.166666666667 0 0 0 ;
           0 -0.333333333333 -0.333333333333 0 -0.166666666667 1.33333333333 -0.333333333333 0 -0.166666666667 ;
           0 0 -0.333333333333 -0.333333333333 0 -0.333333333333 1.33333333333 -0.166666666667 -0.166666666667 ;
           0 0 -0.333333333333 -0.166666666667 0 0 -0.166666666667 0.666666666667 0 ;
           0 0 -0.333333333333 0 0 -0.166666666667 -0.166666666667 0 0.666666666667 ]

    for i in 1:length(jac)
        @test abs(solver.mat[i] - jac[i]) < 1e-5
    end
end

@testset "Assembly" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Float64}(mesh)

    testAssembly(diffusion_system)
end

@testset "AssemblyFD" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Dual{4, Float64}}(mesh)

    testAssembly(diffusion_system)
end


text:
# These tests do not require an FT device which supports D2XX to be connected
#
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.

module TestLibFTD2XX

using Test
using LibFTD2XX
import LibFTD2XX.Wrapper

@testset "high level" begin

  # libversion
  if Sys.iswindows()
    ver = libversion()
    @test ver isa VersionNumber
  else
    @test_throws UndefVarError libversion()
  end

  # createdeviceinfolist
  numdevs = LibFTD2XX.createdeviceinfolist()
  @test numdevs == 0

  # LibFTD2XX.getdeviceinfodetail
  @test_throws D2XXException LibFTD2XX.getdeviceinfodetail(0)

  # FT_HANDLE functions...
  @testset "FT_HANDLE" begin

    # open by description
    @test_throws Wrapper.FT_DEVICE_NOT_FOUND open("", OPEN_BY_DESCRIPTION)

    # open by serialnumber
    @test_throws Wrapper.FT_DEVICE_NOT_FOUND open("", OPEN_BY_SERIAL_NUMBER)

    handle = FT_HANDLE() # create invalid handle...

    # bytesavailable
    @test_throws D2XXException bytesavailable(handle)

    # read
    @test_throws D2XXException read(handle, 0)
    @test_throws D2XXException read(handle, 1)
    @test_throws ErrorException read(handle, -1) # exception type set by Base/io.jl

    # write
    txbuf = ones(UInt8, 10)
    @test_throws D2XXException write(handle, txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws ErrorException write(handle, Int.(txbuf)) # No byte I/O...

    # readavailable
    @test_throws D2XXException readavailable(handle)

    # baudrate
    @test_throws D2XXException baudrate(handle, 9600)
    @test_throws DomainError baudrate(handle, 0)
    @test_throws DomainError baudrate(handle, -1)

    # flush and eof
    @test_throws D2XXException flush(handle)
    @test_throws D2XXException eof(handle)

    # driverversion
    if Sys.iswindows()
      @test_throws D2XXException driverversion(handle)
    else
      @test_throws UndefVarError driverversion(handle)
    end

    # datacharacteristics
    @test_throws D2XXException datacharacteristics(handle, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)

    # timeouts tests...
    timeout_read, timeout_wr = 200, 100 # milliseconds
    @test_throws D2XXException timeouts(handle, timeout_read, timeout_wr)
    @test_throws DomainError timeouts(handle, timeout_read, -1)
    @test_throws DomainError timeouts(handle, -1, timeout_wr)

    # status
    @test_throws D2XXException status(handle)

    # close and isopen
    retval = close(handle)
    @test retval == nothing
    @test !isopen(handle)
    @test LibFTD2XX.Wrapper._ptr(handle) == C_NULL
    retval = close(handle) # check can close more than once without issue...
    @test !isopen(handle)
  end

  # D2XXDevice
  @testset "D2XXDevice" begin

    # Constructor
    @test_throws DomainError D2XXDevice(-1)
    @test_throws D2XXException D2XXDevice(0)

    # D2XXDevices
    devices = D2XXDevices()
    @test length(devices) == numdevs == 0

    device = D2XXDevice(0, 0, 0, 0, 0, "", "", FT_HANDLE()) # blank device...

    # isopen
    @test !isopen(device)

    # open
    @test_throws Wrapper.FT_DEVICE_NOT_FOUND open(device)

    # bytesavailable
    @test_throws D2XXException bytesavailable(device)

    nb = 1

    # read
    @test_throws D2XXException read(device, nb)
    @test_throws ErrorException read(device, -1) # exception type set by Base/io.jl

    # write
    txbuf = ones(UInt8, 10)
    @test_throws D2XXException write(device, txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws ErrorException write(device, Int.(txbuf)) # No byte I/O...

    # readavailable
    @test_throws D2XXException readavailable(device)

    # baudrate
    @test_throws D2XXException baudrate(device, 9600)
    @test_throws DomainError baudrate(device, 0)
    @test_throws DomainError baudrate(device, -1)

    # flush and eof
    @test_throws D2XXException flush(device)
    @test_throws D2XXException eof(device)

    # driverversion
    if Sys.iswindows()
      @test_throws D2XXException driverversion(device)
    else
      @test_throws UndefVarError driverversion(device)
    end

    # datacharacteristics
    @test_throws D2XXException datacharacteristics(device, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)

    # timeouts tests...
    timeout_read, timeout_wr = 200, 100 # milliseconds
    @test_throws D2XXException timeouts(device, timeout_read, timeout_wr)
    @test_throws DomainError timeouts(device, timeout_read, -1)
    @test_throws DomainError timeouts(device, -1, timeout_wr)

    # status
    @test_throws D2XXException status(device)

    # close and isopen (all devices)
    retval = close.(devices)
    @test all(retval .== nothing)
    @test all(.!isopen.(devices))
    @test all(LibFTD2XX.Wrapper._ptr.(fthandle.(devices)) .== C_NULL)
    close.(devices) # check can close more than once without issue...
    @test all(.!isopen.(devices))
  end

end

end # module TestLibFTD2XX


text:
include("plotter_setup.jl")
data_folder_name = "propagation_qsw"

## CTQW vs GQSW vs LQSW
cla()
figure(figsize=[5,2.5])
prob_global = npzread("$data_folder_name/global.npz")
prob_local = npzread("$data_folder_name/local.npz")
prob_ctqw = npzread("$data_folder_name/ctqw.npz")
@assert length(prob_global) == length(prob_local) == length(prob_ctqw)
n = length(prob_global)
divn=div(n, 2)
positions = -divn:divn
plot(positions, prob_ctqw, "-r", label="CTQW", linewidth=1)
plot(positions, prob_global, "--k", label="GQSW")
plot(positions, prob_local, "-.b", label="LQSW")
legend(loc=1)
xlabel("position")
ylabel("measurement probability")
ylim(0., 0.055)
xlim(-divn-5, divn+5)
savefig("plots/nonmoralizing_qsw/global_vs_local_prob_dist.pdf", bbox_inches = "tight", metadata = Dict("CreationDate" => nothing))

## scaling LQSW
linestyles = [":g", "-.b", "-k"]
ω_select = [1,3,4,6]
ωs = (0.:0.2:1.)[ω_select]
max_n = 300
t_scales = 0.02:0.02:1
data = npzread("$data_folder_name/sec_moments_local.npz")[ω_select,:]
data_t = collect(t_scales .* max_n)[1:size(data,2)]

cla()
figure(figsize=[2.8, 2])
batch = 5

linetyles = ["--", ":", "-.", "-", "-", "--"]
for (i, ω) = enumerate(ωs)
    _, data_y = scaling_exponent_generator(log.(data_t), data[i,:], batch)
    plot(data_t[div(batch, 2)+1:end-div(batch, 2)], data_y, linestyle=linetyles[i], label="\$\\omega = $ω\$", linewidth=1.3)
end

legend(loc=4, fontsize=7, bbox_to_anchor=(1., 0.15))
ylim(0.9, 2.05)
#hlines([1,2], minimum(data_t), maximum(data_t), linestyles="--", linewidth=.8)

xlabel("time \$t\$")
ylabel("scaling exponent \$\\alpha \$")
savefig("plots/nonmoralizing_qsw/scaling_exponent_local.pdf", bbox_inches = "tight", metadata = Dict("CreationDate" => nothing))

## scaling LQSW
ω_select = [1,3,4,6]
ωs = (0.:0.2:1.)[ω_select]
max_n = 300
t_scales = 0.02:0.02:1
data = npzread("$data_folder_name/sec_moments_global.npz")[ω_select,:]
data_t = collect(t_scales .* max_n)[1:size(data,2)]

cla()
figure(figsize=[2.8, 2])
batch = 5

linetyles = ["--", ":", "-.", "-", "-", "--"]
for (i, ω) = enumerate(ωs)
    _, data_y = scaling_exponent_generator(log.(data_t), data[i,:], batch)
    plot(data_t[div(batch, 2)+1:end-div(batch, 2)], data_y, linestyle=linetyles[i], label="\$\\omega = $ω\$", linewidth=1.3)
end

legend(loc=4, fontsize=7, bbox_to_anchor=(1., 0.15))
ylim(0.9, 2.05)
#hlines([1,2], minimum(data_t), maximum(data_t), linestyles="--", linewidth=.8)
yticks(visible=false)

xlabel("time \$t\$")
savefig("plots/nonmoralizing_qsw/scaling_exponent_global.pdf", bbox_inches = "tight", metadata = Dict("CreationDate" => nothing))


text:
using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"

indir = projectdir("notebooks")
outdir = projectdir("scripts")
!isdir(outdir) && mkdir(outdir)

function copy_file(
  fromfile::AbstractString, fromdir::AbstractString, todir::AbstractString
)

  nblines = readlines(joinpath(fromdir, fromfile))
  outfilename = joinpath(todir, fromfile)

  isfile(outfilename) && rm(outfilename)
  outfile = open(outfilename, "w")

  for i in nblines
    #println(i)
    if length(i) > 0 && i[1] == '#'
      continue
    else
      write(outfile, i)
      write(outfile, "\n")
    end
  end

  close(outfile)

end

cd(indir) do
  nbsubdirs = readdir()
  for nbsubdir in nbsubdirs
    if isdir(nbsubdir) && nbsubdir !== "scripts"
      if !isdir(joinpath(outdir, nbsubdir))
        mkdir(joinpath(outdir, nbsubdir))
      end
      # Find all notebooks in nbsubdir, skip subdirs "05" and  "intros"
      if !(nbsubdir in ["05", "intros"])
        nbs = readdir(nbsubdir)
        println("$(nbsubdir): $(nbs)\n")
        # Copy the notebooks to the scripts dir
        for nb in nbs
          nb == ".DS_Store" && continue
          copy_file(nb, joinpath(indir, nbsubdir), joinpath(outdir, nbsubdir))
        end
      else # Handle subdirs in "05" and "intros"
        nb_intro_dirs = readdir(nbsubdir)
        println("$(nbsubdir): $(nb_intro_dirs)\n")
        # Copy the notebooks to the scripts dir
        for nb_intro_dir in nb_intro_dirs
          nb_intro_dir == ".DS_Store" && continue
          intro_nbs = readdir(joinpath(indir, nbsubdir, nb_intro_dir))
          println("$(nbsubdir)/$(nb_intro_dir): $(intro_nbs)\n")
          for nb in intro_nbs
            nb == ".DS_Store" && continue
            #println("Checking for $(joinpath(outdir, nbsubdir, nb_intro_dir))")
            !isdir(joinpath(outdir, nbsubdir, nb_intro_dir)) &&
              mkdir(joinpath(outdir, nbsubdir, nb_intro_dir))
            copy_file(nb, joinpath(indir, nbsubdir, nb_intro_dir),
              joinpath(outdir, nbsubdir, nb_intro_dir))
          end
        end
      end
    else
      println(("$(nbsubdir) ignored"))
    end
  end

end


text:
struct BiPHSplineData <: AbstractData
    id::Vector{Int}
    entry::Vector{Float64}
    time::Vector{Float64}
    status::Vector{Bool}
    b::Vector{Int}
    z::Vector{Float64}
end

BiPHSplineData() = BiPHSplineData(Vector{Int}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0),
Vector{Bool}(undef,0), Vector{Int}(undef,0), Vector{Float64}(undef,0))

mutable struct BiPHSplineParam <: AbstractParam
    spline::NTuple{2,Spline1D}
    beta::NTuple{2,Float64}
    prev::Float64 ## biomarker-positive subgroup prevalence
end

struct BiPHSplineFixedParam <: AbstractFixedParam
    n::Int
    tmax::Float64
    max_entry::Float64
end

struct BiPHSplinePrior <: AbstractPrior
    spline::NTuple{2,SplineHazard.Prior}
    beta::MultivariateDistribution
    prev::UnivariateDistribution
end

BiPHSplinePrior() = BiPHSplinePrior((SplineHazard.Prior(), SplineHazard.Prior()), MvNormal(10*diagm(0 => ones(Float64, 2))), Beta(2, 2))

BiPHSplineModel() = Model(BiPHSplineParam((SplineHazard.Spline.const_spline(1.0, (0.0, 36.0)), SplineHazard.Spline.const_spline(1.0, (0.0, 36.0))),
                                (log(0.95), log(0.83)), 0.4), FixedParam(400, 36.0, 24.0), generate_data)

BiPHSplineModel(param, fparam, gen_data=generate_data) = Model(param, fparam, gen_data)

struct BiPHSplineFit <: AbstractFit
    s::NTuple{2,SplineHazard.Sampler}
    beta::Array{Float64,2} ## hazard ratios
    prev::Vector{Float64} ## biomarker-positive subgroup prevalence
    accept::Float64 ## acceptance rates
    prior::Prior
end

nevents(data::BiPHSplineData) = sum(data.status)

function slice(i::Int, x::BiPHSplineFit)
    BiPHSplineParam((SplineHazard.extract_spline(x.s[1], i), SplineHazard.extract_spline(x.s[2], i)), (x.beta[i,1], x.beta[i,2]), x.prev[i])
end

function sample_prior(B::Int, prior::BiPHSplinePrior)
    beta = rand(prior.beta, B)
    prev = rand(prior.prev, B)

    [BiPHSplineParam((SplineHazard.sample_prior(first(prior.spline)), SplineHazard.sample_prior(last(prior.spline))),
           (beta[1,b], beta[2,b]), prev[b]) for b in 1:B]
end

## return prior distribution (before sampling) or posterior sample (after sampling)
## of treatment effect parameter
function get_outcome(x::BiPHSplineFit)
    if size(x.beta)[1] == 0
        x.prior.beta
    else
        x.beta
    end
end

function inv_cumhaz(U::Float64, beta, group, s::Spline1D, tmax::Float64)
    f(x) = SplineHazard.Spline.cumhaz(x, s)*exp(group*beta) - U

    if f(tmax) < 0
        Inf
    else
        fzero(f, 0, tmax)
    end
end

function generate_data(par::BiPHSplineParam, fpar::BiPHSplineFixedParam)
    ##predict_data(par, fpar, Data(), 0.0)

    n = fpar.n
    J = length(par.beta)

    z = rand(Binomial(1, 0.5), n) .- 0.5
    b = 1 .+ rand(Binomial(1, par.prev), n)

    ## 25 patients per month
    max_entry = fpar.tmax ##/ 25
    entry = rand(Uniform(0, fpar.max_entry), n)
    ##U = -log.(1 .-rand(Uniform(0, 1), n))
    U = 1 .- rand(Uniform(0, 1), n)

    T = exp.([quantile(Logistic(log(12), 1/6), 1 .- exp.(log(U[i]).*exp.(-par.beta[b[i]]*z[i]))) for i in 1:n])

    ##C = min.(rand(Exponential(fpar.tmax), n), fpar.tmax)
    C = min.(rand(Uniform(0, 3*36), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    BiPHSplineData(id, entry, time, status, b, z)
end

function predict_data(par::BiPHSplineParam, fpar::BiPHSplineFixedParam, data0::BiPHSplineData, caltime::Float64, nevents=NaN)
    if any(isnan.(nevents))
        n = fpar.n
        N = n.*[1-par.prev, par.prev]
    else
        N = Int.(floor.(nevents ./ [1-par.prev, par.prev]))
        n = sum(N)
    end

    tmax = fpar.tmax
    max_entry = n / 15 ##fpar.max_entry

    beta = par.beta
    n0 = length(data0.time)

    z = rand(Binomial(1, 0.5), n) .- 0.5
    b = 1 .+ rand(Binomial(1, par.prev), n)

    if n < n0
        @error "requested sample size must be at least as large as sample size of old data"
    end

    entry = Vector{Float64}(undef,n)
    chaz = zeros(Float64, n)

    if caltime >= max_entry
        if n > n0
            @error "max_entry has passed, cannot recruit more patients"
        end
    else
        copyto!(entry, rand(Uniform(caltime, max_entry), n))
    end

    left = caltime .- entry

    ## keep old data
    if n0 > 0
        z[1:n0] = data0.z
        b[1:n0] = data0.b

        chaz[1:n0] = [SplineHazard.Spline.cumhaz(data0.time[i], par.spline[b[i]]).*exp.(z[i] * beta[b[i]]) for i in 1:n0]
        entry[1:n0] = data0.entry
        left[1:n0] = caltime .- data0.entry
        keep = data0.time .< left[1:n0]
        chaz[1:n0][keep] = zeros(Float64, sum(keep))
    end

    U = .-log.(1 .- rand(Uniform(0, 1), n)) .+ chaz
    T = [inv_cumhaz(U[i], beta[b[i]], z[i], par.spline[b[i]], 500.0) for i in 1:n]

    time = min.(T, tmax)
    status = T .< tmax

    max_id = 0
    id = zeros(Int, n)

    ## keep old data
    if n0 > 0
        keep = data0.time .< left[1:n0]

        if sum(keep) > 0
            keep2 = zeros(Bool, n)
            keep2[1:n0] = keep

            time[keep2] = data0.time[keep]
            status[keep2] = data0.status[keep]
        end

        id[1:n0] = data0.id
        max_id = maximum(data0.id)
    end

    if n > n0
        id[(n0+1):end] = (max_id+1):(max_id+n-n0)
    end

    BiPHSplineData(id, entry, time, status, b, z)
end

function select_data(t::Float64, data::BiPHSplineData; entry_stop=t)
    sel = data.entry .<= entry_stop
    entry2 = data.entry[sel]
    time2 = data.time[sel]
    status2 = data.status[sel]

    cutoff = max.(t .- entry2, 0) .+ 1e-14

    status2 = status2 .& (time2 .<= cutoff)
    time2 = min.(time2, cutoff)

    Data(data.id[sel], entry2, time2, status2, data.b[sel], data.z[sel])
end

function loglik(data::BiPHSplineData, param::BiPHSplineParam)
    log_likelihood(param.spline, param.beta, data)
end

function log_likelihood_b(b::Int, s::Spline1D, beta, theta, data::BiPHSplineData,
                          use_precalc=false, precalc_haz=0.0, precalc_cumhaz=0.0)

    selb = data.b .== b
    time = data.time[selb]
    status = data.status[selb]
    zb = data.z[selb] * beta

    if use_precalc
        h = precalc_haz * exp.(theta)
        ch = precalc_cumhaz * exp.(theta)
    else
        h = SplineHazard.Spline.hazard(time[status], s)
        ch = SplineHazard.Spline.cumhaz(time, s)
    end

    sum(log.(h) .+ zb[status]) - sum(ch.*exp.(zb))
end

function log_likelihood(s::NTuple{2,Spline1D}, beta, prev, theta1, theta2, data::BiPHSplineData, use_precalc=false,
                        precalc_haz1=0.0, precalc_cumhaz1=0.0,
                        precalc_haz2=0.0, precalc_cumhaz2=0.0)

    ll0 = log_likelihood_b(1, s[1], beta[1], theta1, data, use_precalc, precalc_haz1, precalc_cumhaz1)
    ll1 = log_likelihood_b(2, s[2], beta[2], theta2, data, use_precalc, precalc_haz2, precalc_cumhaz2)

    ll0 + sum(data.b .== 1)*log(1 - prev) + ll1 + sum(data.b .== 2)*log(prev)
end

function laplace_approx(data::BiPHSplineData, param::BiPHSplineParam, prior::BiPHSplinePrior)
    sel1 = data.b .== 1
    sel2 = data.b .== 2

    precalc_haz1, precalc_cumhaz1 = SplineHazard.Spline.eval_basis(param.spline[1], data.time[data.status .& sel1], data.time[sel1])
    precalc_haz2, precalc_cumhaz2 = SplineHazard.Spline.eval_basis(param.spline[2], data.time[data.status .& sel2], data.time[sel2])

    J, K1, K2 = length(param.beta), size(precalc_haz1)[2], size(precalc_haz2)[2]

    theta10 = log.(param.spline[1].c)
    theta20 = log.(param.spline[2].c)
    beta0 = param.beta

    function f(x)
        ## prev, v1, theta1, beta1, v2, theta2, beta2
        prev = expit(x[1])
        v1 = exp(2*x[2])
        theta1 = x[3:(K1+2)]
        beta1 = x[K1+3]

        v2 = exp(2*x[K1+4])
        theta2 = x[(K1+4+1):(K1+4+K2)]
        beta2 = x[end]

        -(log_likelihood(param.spline, [beta1, beta2], prev, theta1, theta2, data, true,
                         precalc_haz1, precalc_cumhaz1, precalc_haz2, precalc_cumhaz2)
          - dot(theta1,theta1)/(2*v1) - K1*log(v1)/2 - dot(theta2,theta2)/(2*v2) - K2*log(v2)/2) + logpdf(prior.prev, prev) + log(prev) + log(1-prev) + x[2] + x[K1+4]
    end

    res = optimize(f, [logit(0.5), 0.0, theta10..., beta0[1], 0.0, theta20..., beta0[2]])

    @show res

    est = Optim.minimizer(res)
    H = Symmetric(ForwardDiff.hessian(f, est))

    beta1 = est[K1+3]
    beta2 = est[end]

    H_chol = cholesky!(H; check=false)

    if issuccess(H_chol)
        L = inv(H_chol.U')
        Sigma = (L'*L)[[K1+3,end],[K1+3,end]]
        [Normal(beta1, sqrt(Sigma[1,1])), Normal(beta2, sqrt(Sigma[2,2]))]
        ##MvNormal([beta1, beta2], (L'*L)[[4,7],[4,7]])
    else
        [Dirac.DiracPM(beta1), Dirac.DiracPM(beta2)]
    end
end

function sample_param1(s::NTuple{2,Spline1D}, beta, prev, data, prior, V, t)
    function log_target(x)
        p = expit(x[1])
        log_likelihood(s, beta, p, [0.0], [0.0], data) + logpdf(prior.prev, p) + log(p) + log(1-p)
    end

    AM_RR_single(log_target, [logit(prev)], V, t)
end

function sample_param2(s::NTuple{2,Spline1D}, beta, prev, data, prior, V, t)
    function log_target(x)
        p = expit(x[3])
        log_likelihood(s, x[1:2], p, [0.0], [0.0], data) + logpdf(prior.beta, x[1:2]) + logpdf(prior.prev, p) + log(p) + log(1-p)
    end

    AM_RR_single(log_target, vcat(beta, logit(prev)), V, t)
end

function get_approx(param::BiPHSplineParam, new_data::BiPHSplineData, prior::BiPHSplinePrior)
    ##laplace_approx(new_data, param, prior)

    x = fit(5000, new_data; warmup=2000, prior=prior)
    [Normal(mean(x.beta[:,1]), sqrt(var(x.beta[:,1]))), Normal(mean(x.beta[:,2]), sqrt(var(x.beta[:,2])))]

    # S = cov(x.beta)
    # S_chol = cholesky!(S; check=false)

    # if !issuccess(S_chol)
    #     Dirac.MultivariateDiracPM(mean(x.beta, dims=1)[1,:])
    # else
    #     MvNormal(mean(x.beta, dims=1)[1,:], Distributions.PDMats.PDMat(S_chol))
    # end
end

logit(p) = log(p) - log(1-p)
expit(x) = exp(x) / (1+exp(x))

function fit(M::Int, data::BiPHSplineData; prior::BiPHSplinePrior=BiPHSplinePrior(), warmup::Int=Int(floor(M/2)))
    n = length(data.time)

    if M == 0 || n == 0
        s1 = SplineHazard.create_sampler((sp, x) -> log_likelihood_b(1, sp, x, [0.0], data))
        s2 = SplineHazard.create_sampler((sp, x) -> log_likelihood_b(2, sp, x, [0.0], data))
        return Fit((s1, s2), Matrix{Float64}(undef,0,0), Vector{Float64}(undef,0), 0.0, prior)
    end

    s1 = SplineHazard.setup_sampler(M, data.time[(data.b .== 1) .& data.status], (sp, x) -> log_likelihood_b(1, sp, x, [0.0], data))
    s2 = SplineHazard.setup_sampler(M, data.time[(data.b .== 2) .& data.status], (sp, x) -> log_likelihood_b(2, sp, x, [0.0], data))

    sp1 = SplineHazard.extract_spline(s1, 1)
    sp2 = SplineHazard.extract_spline(s2, 1)

    J = length(prior.beta)
    beta = Array{Float64,2}(undef,M+1,J)
    beta[1,:] = zeros(Float64, J)

    prev = Vector{Float64}(undef, M+1)
    prev[1] = 0.5

    if prior.beta isa MultivariateDiracPM
        V = diagm(0 => ones(Float64, 1))
        m = [logit(prev[1])]
    else
        V = diagm(0 => ones(Float64, J+1))
        m = vcat(beta[1,:], logit(prev[1]))
    end

    eps = 0.001*V

    attempt = 0
    accept = 0
    ac = 0

    for t in 1:M
        attempt += 1

        if prior.beta isa MultivariateDiracPM
            par, ac = sample_param1((sp1, sp2), beta[t,:], prev[t], data, prior, V, t)
            beta[t+1,:] = prior.beta.x
        else
            par, ac = sample_param2((sp1, sp2), beta[t,:], prev[t], data, prior, V, t)
            beta[t+1,:] = par[1:2]
        end

        prev[t+1] = expit(par[end])
        new_m = (par + m*t) / (t+1)
        V = updateVar(t+1, par, new_m, m, V, eps)
        m = new_m
        accept += ac

        SplineHazard.update!(s1, beta[t+1,1])
        sp1 = SplineHazard.extract_spline(s1, t)

        SplineHazard.update!(s2, beta[t+1,2])
        sp2 = SplineHazard.extract_spline(s2, t)
    end

    index = (warmup+1):M
    BiPHSplineFit((s1, s2), beta[index,:], prev[index], accept/attempt, prior)
end

function summary(x::BiPHSplineFit)
    B, J = size(x.beta)

    ##SplineHazard.summary(x, time)

    DataFrame(variable=[repeat(["beta"], inner=J)..., "prev"],
              mean=[mean(x.beta,dims=1)[1,:]..., mean(x.prev)],
              stderr=sqrt.([var(x.beta,dims=1)[1,:]..., var(x.prev)]),
              neff=[[effective_sample_size(x.beta[:,j]) for j in 1:J]..., effective_sample_size(x.prev)])
end

function fit(chains::Int, M::Int, data::BiPHSplineData; prior::BiPHSplinePrior=BiPHSplinePrior(), warmup::Int=Int(floor(M/2)))
    @distributed vcat for k in 1:chains
        fit(M, data; warmup=warmup, prior=prior)
    end
end

function merge_chains(a::BiPHSplineFit, b::BiPHSplineFit)
    Fit(a.s, vcat(a.beta, b.beta), vcat(a.prev, b.prev), 0.5*(a.accept + b.accept), a.prior)
end

merge_chains(a) = a
merge_chains(a...) = foldr(merge_chains, a)

function summary(x::Vector{BiPHSplineFit})
    chains = length(x)

    z = merge_chains(x...)
    df = summary(z)

    if chains > 1
        tmp = hcat(z.beta, z.prev)
        len = Int(floor(size(x[1].beta, 1)/2))
        df.Rhat = [potential_scale_reduction([tmp[(1+(k-1)*len):(k*len),j] for k in 1:(2*chains)]...) for j in 1:size(tmp,2)]
    end

    df
end

function sim(B::Int, M::Int; warmup=Int(floor(M/2)),
             m=BiPHSplineModel(), prior=BiPHSplinePrior())
    res = Array{Float64, 2}(undef,B, 4)

    J = length(m.param.beta)
    beta = SharedArray{Float64,2}((B,J))
    var_beta = SharedArray{Float64,2}((B,J))
    cpb = SharedArray{Bool,2}((B,J))

    prev = SharedArray{Float64,1}(B)
    var_prev = SharedArray{Float64,1}(B)
    cpp = SharedArray{Bool,1}(B)

    @sync @distributed for b in 1:B
        data = generate_data(m.param, m.fparam)
        data = select_data(m.fparam.tmax, data)
        x = fit(M, data, warmup=warmup, prior=prior)

        beta[b,:] = mean(x.beta, dims=1)
        var_beta[b,:] = var(x.beta, dims=1)
        cpb[b,:] = [(quantile(x.beta[:,j], 0.025) < m.param.beta[j]) .& (quantile(x.beta[:,j], 0.975) > m.param.beta[j]) for j in 1:length(m.param.beta)]

        prev[b,:] = mean(x.prev, dims=1)
        var_prev[b,:] = var(x.prev, dims=1)
        cpp[b,:] = [(quantile(x.prev, 0.025) < m.param.prev) .& (quantile(x.prev, 0.975) > m.param.prev) for j in 1:length(m.param.prev)]
    end

    DataFrame(variable=[repeat(["beta"], inner=J)..., "prev"],
              mean=[mean(beta,dims=1)[1,:]..., mean(prev)],
              sse=sqrt.([var(beta,dims=1)[1,:]..., var(prev)]),
              ese=sqrt.([mean(var_beta,dims=1)[1,:]..., mean(var_prev)]),
              cp=[mean(cpb,dims=1)[1,:]..., mean(cpp)])
end

function biphspline_test_chains(M, chains)
    model = BiPHSplineModel()
    data = model.generate_data(model.param, model.fparam)
    x = fit(chains, M, data)
    summary(x)
end

function biphspline_test_laplace_approx()
    m = BiPHSplineModel()
    data = m.generate_data(m.param, m.fparam)

    la = laplace_approx(data, m.param, BiPHSplinePrior())
    x = fit(5000, data; warmup=2000)

    df = summary(x)
    la, df
end


text:
VERSION ≤ v"0.7-" && __precompile__()
module ForceImport

#   This file is part of ForceImport.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

export @force, @merge

"""

        @force using Module

Forces imports of exported methods from `Module`, even if there are conflicts.
"""
macro force(use)
    use.head ∉ [:using,:toplevel] && throw(error("$use not a using statement"))
    pkgs = use.head ≠ :using ? use.args : [use]
    out = []
    for p ∈ pkgs
        w = typeof(p.args[end]) == Symbol ? [p.args[end]] : p.args[end].args
        m = VERSION > v"0.7-" ? p.args[end].args[end] : p.args[end]
        s = if VERSION > v"0.7-"
            :([Expr(:import,Expr(:.,Symbol.($(string.(w)))...,j)) for j ∈ names($(esc(m)))])
        else
            :([Expr(:import,Symbol.($(string.(w)))...,j) for j ∈ names($(esc(m)))])
        end
        push!(out,Expr(:import,p.args...),:($(esc(:eval))(Expr(:toplevel,$s...))))
    end
    return Expr(:block,out...)
end

"""

        ForceImport.@port using Module

Force exports methods from `Module` in current module.
"""
macro port(use)
    use.head ∉ [:using,:toplevel] && throw(error("$use not a using statement"))
    pkgs = use.head ≠ :using ? use.args : [use]
    out = []
    for p ∈ pkgs
        w = typeof(p.args[end]) == Symbol ? [p.args[end]] : p.args[end].args
        m = VERSION > v"0.7-" ? p.args[end].args[end] : p.args[end]
        s = if VERSION > v"0.7-"
            :([Expr(:export,Expr(:.,Symbol.($(string.(w)))...,j)) for j ∈ names($(esc(m)))])
        else
            :([Expr(:export,Symbol.($(string.(w)))...,j) for j ∈ names($(esc(m)))])
        end
        push!(out,Expr(:export,p.args...),:($(esc(:eval))(Expr(:toplevel,$s...))))
    end
    return Expr(:block,out...)
end

"""

        @merge fun(x) = x

Automatically imports `fun` from its module if necessary for function definition.
"""
macro merge(expr)
    if !( (expr.head == :function) | ( (expr.head == :(=)) &&
            (typeof(expr.args[1]) == Expr) && (expr.args[1].head == :call) ) )
        throw(error("ForceImport: $expr is not a function definition"))
    end
    fun = expr.args[1].args[1]
    return Expr(:quote,quote
        for name in names(current_module())
            try
                eval(ForceImport.imp($(string(fun)),name))
            catch
            end
        end
        eval($(Expr(:quote,expr)))
    end)
end

function imp(fun::String,name::Symbol)
    :(Symbol($fun) ∈ names($name) && (import $name.$(Symbol(fun))))
end

#__init__() = nothing

end # module


text:
module Agg

using Cxx
import SFrames.SFrameMod: SFrame
import SFrames.Util: cstring

for f in  [:SUM, :MAX, :MIN, :MEAN, :AVG, :VAR, :VARIANCE, :STD, :STDV, :SELECT_ONE, :CONCAT]
    cxxstr = "aggregate::$f(\$(cstring(col)));"
    macro_e = Expr(:macrocall, Symbol("@icxx_str"), cxxstr)
    @eval function $f(col)
        $macro_e
    end
end

function COUNT()
    icxx"aggregate::COUNT();"
end

function QUANTILE(col, quantile)
    icxx"aggregate::QUANTILE($(cstring(col)), $(Float64(quantile)));"
end

function QUANTILE(col, quantiles::AbstractVector)
    c = icxx"vector<double>();"
    for _ in quantiles
        icxx"$c.push_back($(Float64(_)));"
    end
    icxx"aggregate::QUANTILE($(cstring(col)), $c);"
end

# type ParityAgg
#     parity::Int
#     sum::Int
# end
#
# emit(p::ParityAgg ) = p.parity
# function add_element_simple(p::ParityAgg, elem)
#     p.sum = p.sum + Int(elem)
# end
#
# function partial_finalize(p::ParityAgg)
#     p.parity = p.sum%2
# end

end


text:
struct TobogganGeology
    map::Array{Char, 2}
end

Base.getindex(a::TobogganGeology, i::Int, j::Int) = let
    new_j = j % size(a.map, 2)
    new_j = new_j == 0 ? size(a.map, 2) : new_j
    a.map[i, new_j]
end

text = reduce(vcat, permutedims.(collect.(readlines("day3/input.txt"))))

t = TobogganGeology(text)

function run_slope_guy(right::Int, down::Int)
    global t
    start = (1, 1)
    trees = 0
    while start[1] <= size(t.map, 1)
        trees += t[start...] == '#'
        start = (start[1] + down, start[2] + right)
    end
    return trees
end

println(run_slope_guy(3, 1))

println(run_slope_guy(1,1)*run_slope_guy(3,1)*run_slope_guy(5,1)*run_slope_guy(7,1)*run_slope_guy(1,2))


text:
#####
##### These functions return the transforms required to solve Poisson's equation with
##### periodic boundary conditions or staggered Neumann boundary conditions.
#####
##### Fast Fourier transforms (FFTs) are used in the periodic dimensions and
##### real-to-real discrete cosine transforms are used in the wall-bounded dimensions.
##### Note that the DCT-II is used for the DCT and the DCT-III for the IDCT
##### which correspond to REDFT10 and REDFT01 in FFTW.
#####
##### They operate on an array with the shape of `A`, which is needed to plan
##### efficient transforms. `A` will be mutated.
#####

function plan_forward_transform(A::Array, ::Periodic, dims, planner_flag=FFTW.PATIENT)
    length(dims) == 0 && return nothing
    return FFTW.plan_fft!(A, dims, flags=planner_flag)
end

function plan_forward_transform(A::Array, ::Bounded, dims, planner_flag=FFTW.PATIENT)
    length(dims) == 0 && return nothing
    return FFTW.plan_r2r!(A, FFTW.REDFT10, dims, flags=planner_flag)
end

function plan_backward_transform(A::Array, ::Periodic, dims, planner_flag=FFTW.PATIENT)
    length(dims) == 0 && return nothing
    return FFTW.plan_ifft!(A, dims, flags=planner_flag)
end

function plan_backward_transform(A::Array, ::Bounded, dims, planner_flag=FFTW.PATIENT)
    length(dims) == 0 && return nothing
    return FFTW.plan_r2r!(A, FFTW.REDFT01, dims, flags=planner_flag)
end

non_batched_topologies = ((Periodic, Bounded, Periodic),
                          (Periodic, Bounded, Bounded),
                          (Bounded, Periodic, Bounded),
                          (Bounded, Bounded, Periodic),
                          (Bounded, Bounded, Bounded))

function plan_forward_transform(A::CuArray, topo, dims, planner_flag)
    length(dims) == 0 && return nothing
    return CUDA.CUFFT.plan_fft!(A, dims)
end

function plan_backward_transform(A::CuArray, topo, dims, planner_flag)
    length(dims) == 0 && return nothing
    return CUDA.CUFFT.plan_ifft!(A, dims)
end

forward_orders(::Type{Periodic}, ::Type{Bounded} , ::Type{Bounded})  = (3, 2, 1)
forward_orders(::Type{Periodic}, ::Type{Bounded} , ::Type{Periodic}) = (2, 1, 3)
forward_orders(::Type{Bounded} , ::Type{Periodic}, ::Type{Bounded})  = (1, 3, 2)
forward_orders(::Type{Bounded} , ::Type{Bounded} , ::Type{Periodic}) = (1, 2, 3)
forward_orders(::Type{Bounded} , ::Type{Bounded} , ::Type{Bounded})  = (1, 2, 3)

backward_orders(::Type{Periodic}, ::Type{Bounded} , ::Type{Bounded})  = (1, 2, 3)
backward_orders(::Type{Periodic}, ::Type{Bounded} , ::Type{Periodic}) = (3, 1, 2)
backward_orders(::Type{Bounded} , ::Type{Periodic}, ::Type{Bounded})  = (2, 1, 3)
backward_orders(::Type{Bounded} , ::Type{Bounded} , ::Type{Periodic}) = (3, 1, 2)
backward_orders(::Type{Bounded} , ::Type{Bounded} , ::Type{Bounded})  = (1, 2, 3)

" Used by FFTBasedPoissonSolver "
function plan_transforms(arch, grid::RegularRectilinearGrid, storage, planner_flag)
    Nx, Ny, Nz = size(grid)
    topo = topology(grid)
    periodic_dims = findall(t -> t == Periodic, topo)
    bounded_dims = findall(t -> t == Bounded, topo)

    if arch isa GPU && topo in non_batched_topologies

        rs_storage = reshape(storage, (Ny, Nx, Nz))
        forward_plan_x = plan_forward_transform(storage   , topo[1](), [1], planner_flag)
        forward_plan_y = plan_forward_transform(rs_storage, topo[2](), [1], planner_flag)
        forward_plan_z = plan_forward_transform(storage   , topo[3](), [3], planner_flag)

        backward_plan_x = plan_backward_transform(storage   , topo[1](), [1], planner_flag)
        backward_plan_y = plan_backward_transform(rs_storage, topo[2](), [1], planner_flag)
        backward_plan_z = plan_backward_transform(storage   , topo[3](), [3], planner_flag)

        forward_plans = (forward_plan_x, forward_plan_y, forward_plan_z)
        backward_plans = (backward_plan_x, backward_plan_y, backward_plan_z)
        f_order = forward_orders(topo...)
        b_order = backward_orders(topo...)

        forward_transforms = (
            DiscreteTransform(forward_plans[f_order[1]], Forward(), arch, grid, [f_order[1]]),
            DiscreteTransform(forward_plans[f_order[2]], Forward(), arch, grid, [f_order[2]]),
            DiscreteTransform(forward_plans[f_order[3]], Forward(), arch, grid, [f_order[3]])
        )

        backward_transforms = (
            DiscreteTransform(backward_plans[b_order[1]], Backward(), arch, grid, [b_order[1]]),
            DiscreteTransform(backward_plans[b_order[2]], Backward(), arch, grid, [b_order[2]]),
            DiscreteTransform(backward_plans[b_order[3]], Backward(), arch, grid, [b_order[3]])
        )

    else
        # This is the case where batching transforms is possible. It's always possible on the CPU
        # since FFTW is awesome so it includes all topologies on the CPU.
        #
        # On the GPU batching is possible when the topology is not one of non_batched_topologies
        # (where an FFT is needed along dimension 2), so it includes (Periodic, Periodic, Periodic),
        # (Periodic, Periodic, Bounded), and (Bounded, Periodic, Periodic).

        forward_periodic_plan = plan_forward_transform(storage, Periodic(), periodic_dims, planner_flag)
        forward_bounded_plan = plan_forward_transform(storage, Bounded(), bounded_dims, planner_flag)

        forward_transforms = (
            DiscreteTransform(forward_bounded_plan, Forward(), arch, grid, bounded_dims),
            DiscreteTransform(forward_periodic_plan, Forward(), arch, grid, periodic_dims)
        )

        backward_periodic_plan = plan_backward_transform(storage, Periodic(), periodic_dims, planner_flag)
        backward_bounded_plan = plan_backward_transform(storage, Bounded(), bounded_dims, planner_flag)

        backward_transforms = (
            DiscreteTransform(backward_periodic_plan, Backward(), arch, grid, periodic_dims),
            DiscreteTransform(backward_bounded_plan, Backward(), arch, grid, bounded_dims)
        )
    end

    transforms = (forward = forward_transforms, backward = backward_transforms)

    return transforms
end

" Used by FourierTridiagonalPoissonSolver "
function plan_transforms(arch, grid::VerticallyStretchedRectilinearGrid, storage, planner_flag)
    Nx, Ny, Nz = size(grid)
    TX, TY, TZ = topo = topology(grid)

    # Limit ourselves to x, y transforms (z uses a tridiagonal solve)
    periodic_dims = findall(t -> t == Periodic, (TX, TY))
    bounded_dims = findall(t -> t == Bounded, (TX, TY))

    if arch isa GPU && topo in non_batched_topologies
        if (TX, TY) == (Periodic, Bounded)
            forward_plan_x = plan_forward_transform(storage, Periodic(), [1], planner_flag)
            forward_plan_y = plan_forward_transform(reshape(storage, (Ny, Nx, Nz)), Bounded(), [1], planner_flag)

            backward_plan_x = plan_backward_transform(storage, Periodic(), [1], planner_flag)
            backward_plan_y = plan_backward_transform(reshape(storage, (Ny, Nx, Nz)), Bounded(),  [1], planner_flag)

            forward_transforms = (
                DiscreteTransform(forward_plan_y, Forward(), arch, grid, [2]),
                DiscreteTransform(forward_plan_x, Forward(), arch, grid, [1])
            )

            backward_transforms = (
                DiscreteTransform(backward_plan_x, Backward(), arch, grid, [1]),
                DiscreteTransform(backward_plan_y, Backward(), arch, grid, [2])
            )

        elseif (TX, TY) == (Bounded, Periodic)
            forward_plan_x = plan_forward_transform(storage, Bounded(), [1], planner_flag)
            forward_plan_y = plan_forward_transform(reshape(storage, (Ny, Nx, Nz)), Periodic(), [1], planner_flag)

            backward_plan_x = plan_backward_transform(storage, Bounded(), [1], planner_flag)
            backward_plan_y = plan_backward_transform(reshape(storage, (Ny, Nx, Nz)), Periodic(),  [1], planner_flag)

            forward_transforms = (
                DiscreteTransform(forward_plan_x, Forward(), arch, grid, [1]),
                DiscreteTransform(forward_plan_y, Forward(), arch, grid, [2])
            )

            backward_transforms = (
                DiscreteTransform(backward_plan_y, Backward(), arch, grid, [2]),
                DiscreteTransform(backward_plan_x, Backward(), arch, grid, [1])
            )

        elseif (TX, TY) == (Bounded, Bounded)
            forward_plan_x = plan_forward_transform(storage, Bounded(), [1], planner_flag)
            forward_plan_y = plan_forward_transform(reshape(storage, (Ny, Nx, Nz)), Bounded(), [1], planner_flag)

            backward_plan_x = plan_backward_transform(storage, Bounded(), [1], planner_flag)
            backward_plan_y = plan_backward_transform(reshape(storage, (Ny, Nx, Nz)), Bounded(),  [1], planner_flag)

            forward_transforms = (
                DiscreteTransform(forward_plan_x, Forward(), arch, grid, [1]),
                DiscreteTransform(forward_plan_y, Forward(), arch, grid, [2])
            )

            backward_transforms = (
                DiscreteTransform(backward_plan_x, Backward(), arch, grid, [1]),
                DiscreteTransform(backward_plan_y, Backward(), arch, grid, [2])
            )
        end
    else
        # This is the case where batching transforms is possible. It's always possible on the CPU
        # since FFTW is awesome so it includes all topologies on the CPU.
        #
        # On the GPU batching is possible when the topology is not one of non_batched_topologies
        # (where an FFT is needed along dimension 2), so it includes (Periodic, Periodic, Periodic),
        # (Periodic, Periodic, Bounded), and (Bounded, Periodic, Periodic).

        forward_periodic_plan = plan_forward_transform(storage, Periodic(), periodic_dims, planner_flag)
        forward_bounded_plan = plan_forward_transform(storage, Bounded(), bounded_dims, planner_flag)

        forward_transforms = (
            DiscreteTransform(forward_bounded_plan, Forward(), arch, grid, bounded_dims),
            DiscreteTransform(forward_periodic_plan, Forward(), arch, grid, periodic_dims)
        )

        backward_periodic_plan = plan_backward_transform(storage, Periodic(), periodic_dims, planner_flag)
        backward_bounded_plan = plan_backward_transform(storage, Bounded(), bounded_dims, planner_flag)

        backward_transforms = (
            DiscreteTransform(backward_periodic_plan, Backward(), arch, grid, periodic_dims),
            DiscreteTransform(backward_bounded_plan, Backward(), arch, grid, bounded_dims)
        )
    end

    transforms = (forward = forward_transforms, backward = backward_transforms)

    return transforms
end


text:
using Printf
using NCDatasets
using PyCall
using PyPlot

using Oceananigans
using Oceananigans.Utils
using JULES

cmocean = pyimport("cmocean")

const plt = PyPlot
plt.ioff()

km = kilometers
Tvar = Energy

ds = NCDataset("three_gas_dry_rising_thermal_bubble_$(Tvar).nc")

for n in 1:length(ds["time"])
    fig, axes = plt.subplots(nrows=2, ncols=3, sharex=true, sharey=true, figsize=(16, 9), dpi=200)

    xC, xF = ds["xC"] / km, ds["xF"] / km
    zC, zF = ds["zC"] / km, ds["zF"] / km

    ρ₀ = ds["ρ₀"][:, 1, :]'
    ρ₁₀ = ds["ρ₁₀"][:, 1, :]'
    ρ₂₀ = ds["ρ₂₀"][:, 1, :]'
    ρ₃₀ = ds["ρ₃₀"][:, 1, :]'
    ρe₀ = ds["ρe₀"][:, 1, :]'

    ρ = ds["ρ"][:, 1, :, n]'
    ρ₁ = ds["ρ₁"][:, 1, :, n]'
    ρ₂ = ds["ρ₂"][:, 1, :, n]'
    ρ₃ = ds["ρ₃"][:, 1, :, n]'
    ρw = ds["ρw"][:, 1, 2:end, n]'
    ρu = ds["ρu"][:, 1, :, n]'
    ρe = ds["ρe"][:, 1, :, n]'

    fig.suptitle("Three gas dry rising thermal bubble: time = $(prettytime(ds["time"][n]))")

    ax_ρ₁ = axes[1, 1]
    im = ax_ρ₁.pcolormesh(xC, zC, ρ₁, cmap=cmocean.cm.dense, vmin=0, vmax=1)
    ax_ρ₁.set_title("ρ₁(x, z)")
    ax_ρ₁.set_ylabel("z (km)")
    fig.colorbar(im, ax=ax_ρ₁, label="kg/m³")

    ax_ρ₂ = axes[1, 2]
    im = ax_ρ₂.pcolormesh(xC, zC, ρ₂, cmap=cmocean.cm.dense, vmin=0, vmax=1)
    ax_ρ₂.set_title("ρ₂(x, z)")
    fig.colorbar(im, ax=ax_ρ₂, label="kg/m³")

    ax_ρ₃ = axes[1, 3]
    im = ax_ρ₃.pcolormesh(xC, zC, ρ₃, cmap=cmocean.cm.dense, vmin=0, vmax=1)
    ax_ρ₃.set_title("ρ₃(x, z)")
    fig.colorbar(im, ax=ax_ρ₃, label="kg/m³")

    ax_u = axes[2, 1]
    im = ax_u.pcolormesh(xF, zC, ρu ./ ρ, cmap=cmocean.cm.balance, vmin=-10, vmax=10)
    ax_u.set_title("u(x, z)")
    ax_u.set_xlabel("x (km)")
    ax_u.set_ylabel("z (km)")
    fig.colorbar(im, ax=ax_u, label="m/s")

    ax_w = axes[2, 2]
    im = ax_w.pcolormesh(xC, zC, ρw ./ ρ, cmap=cmocean.cm.balance, vmin=-10, vmax=10)
    ax_w.set_title("w(x, z)")
    ax_w.set_xlabel("x (km)")
    fig.colorbar(im, ax=ax_w, label="m/s")

    ax_e = axes[2, 3]
    im = ax_e.pcolormesh(xC, zC, (ρe .- ρe₀) ./ ρ, cmap=cmocean.cm.thermal, vmin=0, vmax=600)
    ax_e.set_title("e′(x, z)")
    ax_e.set_xlabel("x (km)")
    fig.colorbar(im, ax=ax_e, extend="max", label="J/kg")

    plt.xlim(-10, 10)
    plt.ylim(0, 10)

    filename = @sprintf("three_gas_dry_rising_thermal_bubble_%s_%03d.png", Tvar, n)
    @info "Saving $filename..."
    plt.savefig(filename)
    plt.close(fig)
end

close(ds)

run(`ffmpeg -r 15 -f image2 -i three_gas_dry_rising_thermal_bubble_$(Tvar)_%03d.png -vcodec libx264 -crf 22 -pix_fmt yuv420p three_gas_dry_rising_thermal_bubble_$Tvar.mp4`)


text:
module Example

using Glimmer, Glimmer.FlexUI
import Luxor

println("Start [$(splitext(basename(@__FILE__))[1])]")

#---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
#---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Luxor")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

#---------------------------------------------------------------
# Define Variables
#---------------------------------------------------------------
angle = addVariable!(app, Variable(name="angle", type="flota64",value=8))
count = addVariable!(app, Variable(name="count", type="flota64",value=16))

image = addVariable!(app, Variable(name="image", type="image", value="", ))

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------
ui = VContainer(
    Card(
        title="Controls",
        content=VContainer(
            Slider(
                text="Angle (π over value)",
                trailingText="[\$()]",
                min=2,
                max=16,
                value=8,
                variable="angle"
            ),
            Slider(
                text="Count",
                trailingText="[\$()]",
                min=1,
                max=30,
                value=16,
                variable="count"
            ),
        ),
    ),
    Card(
        title="Result Image",
        content=VContainer(
            Image(
                source="\$(image)",
            ),
        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control
)
# set the controls for the application
controls!(app, ui)

#---------------------------------------------------------------
# the render function - preparing the image
#---------------------------------------------------------------
function render()
    Luxor.Drawing(600, 400, :svg)
    Luxor.origin()
    Luxor.background("white")

    for θ in range(0, step=π/angle[], length=count[])
        Luxor.gsave()
        Luxor.scale(0.2)
        Luxor.rotate(θ)
        Luxor.translate(350, 0)
        Luxor.julialogo(action=:fill, bodycolor=Luxor.randomhue())
        Luxor.grestore()
    end

    Luxor.gsave()
    Luxor.scale(0.3)
    Luxor.juliacircles()
    Luxor.grestore()

    Luxor.translate(150, -150)
    Luxor.scale(0.3)
    Luxor.julialogo()
    Luxor.finish()

    image["svg"] = Luxor.preview()      # updating the UI with the new image
end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module


text:
using Gadfly, Distributions
using Changepoints

srand(1)

# Simulate time series with changepoints
n = 1000          # Sample size
λ = 25            # Freq of changepoints
μ, σ = Normal(0,1), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)

# Find the Changepoints for a specific penalty
pelt_cps, pelt_cost = @PELT data Normal(?, 1.0) 6.0

# For a range of penalties construct a cost function and use the CROPS function
seg_cost = NormalMeanSegment(data)
crops_output = CROPS(seg_cost , n, (4.0,100.0))
# ...or use the PELT macro
@PELT data Normal(?, 1.0) 4.0 100.0

# 1st plot: time series
p1 = plot(y=data, Geom.line, Theme(default_color=colorant"black"))
draw(PNG("example.png", 5inch, 5inch), p1)

# 2nd plot: time series with chpts
p2 = plot(data, pelt_cps)
draw(PNG("example_pelt.png", 5inch, 5inch), p2)

# 3rd: elbow plot
p3 = plot(crops_output)
draw(PNG("elbowplot.png", 5inch, 5inch), p3)


text:
@testset "1562.find-latest-group-of-size-m.jl" begin
    @test find_latest_step([3,5,1,2,4], 1) == 4
    @test find_latest_step([3,5,1,2,4], 2) == -1
    @test find_latest_step([1], 1) == 1
    @test find_latest_step([2, 1], 2) == 2
end


text:
using Hiccup

typealias AString AbstractString

import Juno: Model, view

render(::Inline, m::Model) = m.data

view(x::AString) = x
view(x::Associative) = x

view(n::Node) =
  d(:type     => :dom,
    :tag      => tag(n),
    :attrs    => n.attrs,
    :contents => map(view, n.children))

render(::Inline, n::Node) = view(n)

render(::Inline, x::HTML) = view(x)

render(::Console, x::Node) =
  @msg result(view(x))

render(::Editor, x::Node) = view(x)

import Juno: Tree, SubTree

render(i::Inline, t::Tree) =
  isempty(t.children) ?
    render(i, t.head) :
    d(:type     => :tree,
      :head     => render(i, t.head),
      :children => map(x->render(i, x),
                       t.children))

render(i::Inline, t::SubTree) =
  d(:type  => :subtree,
    :label => render(i, t.label),
    :child => render(i, t.child))

import Juno: Copyable

render(i::Inline, x::Copyable) =
  d(:type => :copy,
    :view => render(i, x.view),
    :text => x.text)

import Juno: Link, link

render(i::Inline, l::Link) =
  d(:type => :link,
    :file => l.file,
    :line => l.line-1,
    :contents => map(x->render(i, x), l.contents))

render(::Clipboard, l::Link) =
  "$(l.file):$(l.line)"

import Juno: icon, fade


text:
function boltzmannresult2DataFrame(pa::Vector, avec::Vector, varargin...)
    #if string representation for a is provided, check its dimensionality and use it
    #otherwise do not include them into the data-frame
    strings_used = true;
    nargin = length(varargin)
    if nargin == 0
        strings_used = false;
    elseif nargin == 1
        #check size
        astring = varargin[1]
        if length(astring) != length(pa)
            error("String representation and value-vector for a have different lengths.")
        end
    else
        error("Wrong number of arguments. Only one optional argument allowed (vector with a string-representation of a).")
    end

    #fill data frame
    if strings_used
        pa_df = DataFrame(p_a=pa, a=avec, a_string=astring)
    else
        pa_df = DataFrame(p_a=pa, a=avec)
    end

    return pa_df
end

function BAmarginal2DataFrame(pa::Vector, avec::Vector, varagin...)
    return boltzmannresult2DataFrame(pa, avec, varagin...)
end

#the function assumes that pago has one row per a-value and one column per w-value
function BAconditional2DataFrame(pagw::Matrix, avec::Vector, wvec::Vector, varargin...)
    #if string representations for a and w are provided, check their dimensionality and use them
    #otherwise do not include them into the data-frames
    strings_used = true;
    nargin = length(varargin)
    if nargin == 0
        strings_used = false;
    elseif nargin == 2
        #check size (size of astring will be checked in boltmannresult2DataFrame)
        astring = varargin[1]
        wstring = varargin[2]
        if length(astring) != size(pagw,1)
            error("String representation for a and value-matrix for pagw mismatch in size.")
        end
        if length(wstring) != size(pagw,2)
            error("String representation for w and value-matrix for pagw mismatch in size.")
        end
    else
        error("Wrong number of arguments. Either provide string-representations for both a AND w or for neither of them.")
    end

    #map matrix onto a vector
    na = size(pagw,1)
    nw = size(pagw,2)
    pagw_v = vec(pagw) #conversion using column-major convention (column-wise)
    avec_v = vec(repmat(avec,1,nw))
    wvec_v = vec(repmat(wvec',na,1))
    if strings_used
        astring_v = vec(repmat(astring,1,nw))
        wstring_v = vec(repmat(wstring',na,1))
    end

    #fill data frames
    if strings_used
        pagw_df = DataFrame(p_agw=pagw_v, a=avec_v, w=wvec_v, a_string=astring_v, w_string=wstring_v)
    else
        pagw_df = DataFrame(p_agw=pagw_v, a=avec_v, w=wvec_v)
    end

    return pagw_df
end

#the function assumes that pago has one row per a-value and one column per o-value
function BAresult2DataFrame(pa::Vector, pagw::Matrix, avec::Vector, wvec::Vector, varargin...)
    nargin = length(varargin)
    if nargin > 0
        pa_df = BAmarginal2DataFrame(pa,avec,varargin[1])
    else
        pa_df = BAmarginal2DataFrame(pa,avec)
    end

    if nargin > 1
        pagw_df = BAconditional2DataFrame(pagw,avec,wvec,varargin...)
    else
        pagw_df = BAconditional2DataFrame(pagw,avec,wvec)
    end

    return pa_df, pagw_df
end

#convert performance measures to DataFrame representation
function performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj)
    return DataFrame(I_aw=I, H_a=Ha, H_agw = Hagw, E_U = EU, RD_obj = RDobj)
end

#convert performance measures to DataFrame - this method is intended for the three-variable general case
function performancemeasures2DataFrame(I_ow, I_ao, I_awgo, I_aw, Ho, Ha, Hogw, Hago, Hagow, Hagw, EU, ThreeVarRDobj)
    return DataFrame(I_ow=I_ow, I_ao=I_ao, I_awgo=I_awgo, I_aw=I_aw, H_o=Ho, H_a=Ha, H_ogw=Hogw,
                     H_ago=Hago, H_agow=Hagow, H_agw=Hagw, E_U=EU, Objective_value=ThreeVarRDobj)
end


text:
"Return a string of the value as a percentage, rounded to one digit."
perc(x::Real) = @sprintf("%.1f%%", x*100)

"Return a string of `count` and percentage of `total` (in parentheses)."
count_and_perc(count, total) = "$(count)  ($(perc(count/total)))"


text:
struct UncertainVectorValue{T1<:Number, T2<:Number, T3<:Number,
        S<:ValueSupport} <: AbstractUncertainValue
    distribution::Distribution{Multivariate, S}
    value::AbstractVector{T1}
    lower::AbstractVector{T2}
    upper::AbstractVector{T3}
end

dimension(usv::UncertainVectorValue) = length(distribution)

export UncertainVectorValue


text:
# D(u)_function.jl

"""
Implementation of FHIP's susceptibility formula D(u) for the Optical Absorption of Polarons.

See FHIP 1962, equation (35c):
https://link.aps.org/doi/10.1103/PhysRev.127.1004.
"""

"""
ℜD(x::Float64, y::Float64, v::Float64, w::Float64, β::Float64)

    Calculate the real part of D(u) (equation (35c) in FHIP) for a complex argument u = x + iy. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant.
"""
function ℜD(x, y, v, w, β)
    R = (v^2 - w^2) / (w^2 * v)
    P = 1 / (exp(v * β) - 1)
    w^2 / v^2 * (
        R * cos(v * x) * (exp(-v * y) + 2 * P * cosh(v * y)) - (x^2 - y^2) / β - y -
        R * (2 * P + 1)
    )
end

"""
ℑD(x::Float64, y::Float64, v::Float64, w::Float64, β::Float64)

    Calculate the imaginary part of D(u) (equation (35c) in FHIP) for a complex argument u = x + iy. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant.
"""
function ℑD(x, y, v, w, β)
    R = (v^2 - w^2) / (w^2 * v)
    P = 1 / (exp(v * β) - 1)
    w^2 / v^2 * (R * sin(v * x) * (exp(-v * y) - 2 * P * sinh(v * y)) - x * (2 * y / β - 1))
end

# Create a complex number from the real and imaginary parts of D(u).

D(x, y, v, w, β) = ℜD(x, y, v, w, β) + 1im * ℑD(x, y, v, w, β)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SOLAR"
version = v"0.3.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/bbopt/solar.git", "d0c2932a3f25d4b0a71010b87b1f929c6fc1e020")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/solar/src
make COMPILATOR="c++" LIBS="-lm"
mkdir -p $bindir
cp "../bin/solar${exeext}" "${bindir}/solar${exeext}"
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("solar", :solar)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9.1.0", julia_compat="1.6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function _plot_fspec_grid(fs::Dict{Symbol,Dict{Int64,JAXTAM.FFTData}},
        obs_row)
    example = _recursive_first(fs)

    plt_1 = plot(fs; norm=:rms,   save=false, save_csv=true)
    plt_2 = plot(fs; norm=:leahy, save=false)

    plt_3 = plot(fs, norm=:leahy, save=false, freq_lims=(0, 1),     rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 0 to 1 Hz")
    plt_4 = plot(fs, norm=:leahy, save=false, freq_lims=(1, :end),  rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 1 to :end Hz")
    plt_4 = plot(fs, norm=:leahy, save=false, freq_lims=(50, :end), rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 50 to :end Hz")

    Plots.plot(plt_1, plt_2, plt_3, plt_4, layout=grid(4,1), size=(1140,600*4))

    savefig(example, obs_row, example.e_range; plot_name="grid")
end

# Sort of works but doesn't
# macro error_wrapper(mission, obs_row, func)
#     local escaped_obs_row = esc(obs_row)

#     return quote
#         local _wrapped_result = nothing
#         try
#             _wrapped_result = $func
#         catch err
#             if typeof(err) == JAXTAMError
#                 _wrapped_result = err
#                 _log_add($mission, $escaped_obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
#                 @warn err
#             else
#                 rethrow(err)
#             end
#         end

#         _wrapped_result
#     end
# end

function _error_wrapper(func, args...; mission=nothing, obs_row=nothing, kwargs...)
    if mission == nothing
        # Assume first arg is mission
        @assert args[1] isa Mission
        mission = args[1]
    end

    if obs_row == nothing
        # Assume second arg is obs_row or obsid
        if args[2] isa String
            obs_row = master_query(mission, :obsid, args[2])
        elseif args[2] isa DataFrameRow
            obs_row = args[2]
        end
    end

    try
        return func(args...; kwargs...)
    catch err
        if typeof(err) == JAXTAMError
            _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
            @warn err
            return nothing
        else
            rethrow(err)
        end
    end
end

"""
    report(mission::Mission, obs_row::DataFrameRow; e_range::Tuple{Float64,Float64}=_mission_good_e_range(mission), overwrite::Bool, nuke::Bool=false, update_masterpage::Bool=true)

Generates a report for the default energy range

Creates plots for:

    * Lightcurve (+ grouped lightcurves)

    * Periodigram (+ grouped periodograms)

    * Power Spectra

        * :rms, full range, log-rebinned, log-log plot

        * :leahy, full range, log-rebinned, log-log plot

        * :leahy, 0 to 1 Hz, no rebin, linear-linear plot

        * :leahy, 1 to end Hz, no rebin, linear-linear plot

        * :leahy, 50 to end Hz, no regin, linear-linear plot

    * Spectrogram

    * Pulsation search plot

Produces HTML report page

Updates the homepage
"""
function report(mission, obs_row; e_range=_mission_good_e_range(mission), overwrite=false, nuke=false, update_masterpage=true)
    path_jaxtam = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM")
    path_web    = abspath(_obs_path_local(mission, obs_row; kind=:web), "JAXTAM")

    if nuke
        GC.gc() # Required due to Feather.jl loading files lazily, meaning they can't be removed from disk
                # until garbage collection runs and un-lazily-loads them
        ispath(path_jaxtam) ? rm(path_jaxtam, recursive=true) : false
        ispath(path_web)    ? rm(path_web,    recursive=true) : false
    end

    if ismissing(_log_query(mission, obs_row, "meta", :downloaded; surpress_warn=true)) || !_log_query(mission, obs_row, "meta", :downloaded)
        _error_wrapper(download, mission, obs_row)
    end

    if !ismissing(JAXTAM._log_query(mission, obs_row, "errors", :read_cl; surpress_warn=true))
        @warn "Error logged at :read_cl stage, no files to be analysed, skipping report gen"
        return nothing
    end

    images = _log_query(mission, obs_row, "images", e_range)

    img_count_groupless = ismissing(images) ? 0 : size(filter(x->ismissing(x[:group]), images), 1)
    # Expect five 'groupless' plots: lightcurve, periodogram, powerspectra, spectrogram, pulsations
    if img_count_groupless < 5 || overwrite
        lc = _error_wrapper(JAXTAM.lcurve, mission, obs_row, 2.0^0; e_range=e_range)
        # If the initial 1s lightcurve failed, everything else will fail, so just exit now
        lc == nothing ? (return nothing) : ""

        _error_wrapper(JAXTAM.plot, lc; save=true, mission=mission, obs_row=obs_row)
        _error_wrapper(JAXTAM.plot_groups, lc; save=true, size_in=(1140,400/2), mission=mission, obs_row=obs_row)

        pg = _error_wrapper(JAXTAM.pgram, lc, mission=mission, obs_row=obs_row)
        if pg != nothing
            # Only try plotting the periodogram if there is one
            _error_wrapper(JAXTAM.plot, pg; save=true, mission=mission, obs_row=obs_row)
            pg = _error_wrapper(JAXTAM.pgram, lc; per_group=true, mission=mission, obs_row=obs_row)
            _error_wrapper(JAXTAM.plot_groups, pg; save=true, size_in=(1140,600/2), mission=mission, obs_row=obs_row)
        end

        lc = nothing; pg = nothing; GC.gc()

        lc = _error_wrapper(JAXTAM.lcurve, mission, obs_row, 2.0^-13; e_range=e_range)
        lc == nothing ? (return nothing) : ""

        gtis = _error_wrapper(JAXTAM.gtis, mission, obs_row, 2.0^-13; lcurve_data=lc, e_range=e_range)
        gtis == nothing ? (return nothing) : ""
        lc = nothing;

        fs = _error_wrapper(JAXTAM.fspec, mission, obs_row, 2.0^-13, 128; gtis_data=gtis, e_range=e_range)
        if fs != nothing # Only try plotting if fspec didn't fail
            @info "Plotting fspec grid"
            _error_wrapper(JAXTAM._plot_fspec_grid, fs, obs_row, mission=mission, obs_row=obs_row)
            @info "Plotting fspec groups"
            _error_wrapper(JAXTAM.plot_groups, fs; save=true, size_in=(1140,600/2), mission=mission, obs_row=obs_row)
            @info "Plotting sgram"
            _error_wrapper(JAXTAM.plot_sgram, fs;  save=true, size_in=(1140,600), mission=mission, obs_row=obs_row)
            @info "Plotting pulses"
            _error_wrapper(JAXTAM.plot_pulses_candle, fs; save=true, size_in=(1140,600/2), mission=mission, obs_row=obs_row)
            @info "Plotting pulses groups"
            _error_wrapper(JAXTAM.plot_pulses_candle_groups, fs; save=true, size_in=(1140,600/2), mission=mission, obs_row=obs_row)
            fs = 0; GC.gc()
        end
    end

    subpage_path = _error_wrapper(_webgen_subpage, mission, obs_row; e_range=e_range)

    if update_masterpage
        webgen_mission(mission)
    end

    return subpage_path
end

"""
    report(mission::Mission, obsid::String; kwargs...)

Multiple dispath to report(mission::Mission, obs_row::DataFrameRow; kwargs...)

Uses `obsid` to select the observation
"""
function report(mission::Mission, obsid::String; e_range=_mission_good_e_range(mission), overwrite=false, nuke=false, update_masterpage=true)
    obs_row = master_query(mission, :obsid, obsid)

    return report(mission, obs_row; e_range=e_range, overwrite=overwrite, nuke=nuke, update_masterpage=update_masterpage)
end

"""
    report_all(::Mission, ::DataFrameRow; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite::Bool=false, nuke::Bool=false, update_masterpage::Bool=true)

Calls `report` with three default energy ranges
"""
function report_all(mission::Mission, obs_row::DataFrames.DataFrameRow; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite=false, nuke=false, update_masterpage=true)
    if nuke
        path_jaxtam = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM")
        path_web    = abspath(_obs_path_local(mission, obs_row; kind=:web), "JAXTAM")

        GC.gc() # Required due to Feather.jl loading files lazily, meaning they can't be removed from disk
                # until garbace collection runs and un-lazily-loads them
        ispath(path_jaxtam) ? rm(path_jaxtam, recursive=true) : false
        ispath(path_web)    ? rm(path_web,    recursive=true) : false
    end

    for e_range in e_ranges
        report(mission, obs_row; e_range=e_range, overwrite=overwrite, update_masterpage=update_masterpage)
    end

    # Run it again to re-generate report pages so that they all have each others links in them
    for e_range in e_ranges
        println("Report path:\n\t$(report(mission, obs_row; e_range=e_range, overwrite=false))")
    end
end

"""
    report_all(mission::Mission, obsid::String; kwargs...)

Multiple dispath to report_all(mission::Mission, obs_row::DataFrameRow; kwargs...)

Uses `obsid` to select the observation
"""
function report_all(mission::Mission, obsid::String; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite=false, nuke=false, update_masterpage=true)
    obs_row = master_query(mission, :obsid, obsid)

    return report_all(mission, obs_row; e_ranges=e_ranges, overwrite=overwrite, nuke=nuke, update_masterpage=update_masterpage)
end