text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

###################################################################################################
# Batch definition

# Assign available batch + blob clients to available pools
function divide_clients_among_pools()
    num_pools = parse(Int, __params__["_POOL_COUNT"])
    num_credentials = length(__credentials__)

    if num_pools == num_credentials
        return __credentials__, __clients__, __resources__
    elseif num_pools < num_credentials
        return __credentials__[1:num_pools], __clients__[1:num_pools], __resources__[1:num_pools]
    else
        credentials = Array{Any}(undef, 0)
        clients = Array{Any}(undef, 0)
        resources = Array{Any}(undef, 0)
        count = 1
        for i=1:num_pools
            push!(credentials, __credentials__[count])
            push!(clients, __clients__[count])
            push!(resources, __resources__[count])
            if count >= num_credentials
                count = 1
            else
                count += 1
            end
        end
        return credentials, clients, resources
    end
end

"""
    `create_pool(; enable_auto_scale=false, auto_scale_formula=nothing,
    auto_scale_evaluation_interval_minutes=nothing, image_resource_id=nothing)`

 Create a batch pool using parameters specified in as a JSON file. E.g. set `ENV["PARAMETERS"] = "parameters.json"` to
 point to the parameter file.

 *Optional arguments*:

 - `enable_auto_scale` (Bool): Enable auto-scaling of the pool. Requires the `auto_scale_formula` and
    `auto_scale_evaluation_interval_minutes` to be set. If the number of VMs in the pool cannot be specified.

 - `auto_scale_formula` (String): Formula for auto-scaling the pool.
    See "https://docs.microsoft.com/en-us/azure/batch/batch-automatic-scaling" for details.

 - `image_resource_id` (String): Image ID of the VM image.

 *Output*:

 - `Nothing`

 See also:  [`create_pool_and_resource_file`](@ref)
 """
function create_pool(; enable_auto_scale=false, auto_scale_formula=nothing,
    auto_scale_evaluation_interval_minutes=nothing, image_resource_id=nothing,
    container_registry=nothing)

    # Autoscaling?
    for client in __clients__
        create_blob_containers(client["blob_client"], [__container__])
    end

    if enable_auto_scale
        num_nodes = nothing
    else
        num_nodes = __params__["_NODE_COUNT_PER_POOL"]
    end

    # Internode?
    enable_inter_node = parse(Bool, __params__["_INTER_NODE_CONNECTION"])

    # Container?
    if __params__["_CONTAINER"] == "None"
        docker_container = nothing
    else
        docker_container = __params__["_CONTAINER"]
    end

    # Divide clients among pools
    if ~isnothing(__credentials__)
        credential_per_pool, clients_per_pool, resources_per_pool = divide_clients_among_pools()
    else
        __verbose__ && @warn "Could not create pool(s). No credentials were found."
        return nothing
    end
    num_pools = parse(Int, __params__["_POOL_COUNT"])

    for i=1:num_pools
        pool_id = join([__params__["_POOL_ID"], "_", i])
        try
            azureclusterlesshpc.create_pool(clients_per_pool[i]["batch_client"], pool_id, __params__["_POOL_VM_SIZE"], num_nodes,
                __params__["_NODE_OS_PUBLISHER"], __params__["_NODE_OS_OFFER"], __params__["_NODE_OS_SKU"],
                enable_inter_node=enable_inter_node, enable_auto_scale=enable_auto_scale, auto_scale_formula=auto_scale_formula,
                auto_scale_evaluation_interval_minutes=auto_scale_evaluation_interval_minutes, image_resource_id=image_resource_id,
                container=docker_container, container_registry=container_registry)

            # Keep track of active pools
            __verbose__ && print(join(["Created pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"], " with ", num_nodes, " nodes.\n"]))
            push!(__active_pools__, Dict("pool_id" => pool_id, "clients" => clients_per_pool[i], "credentials" => credential_per_pool[i],
                "resources" => resources_per_pool[i]))
        catch e
            if typeof(e) == PyCall.PyError && e.val.error.code == "PoolExists"
                __verbose__ && print(join(["Pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"]," already exists.\n"]))
                push!(__active_pools__, Dict("pool_id" => pool_id, "clients" => clients_per_pool[i], "credentials" => credential_per_pool[i],
                    "resources" => resources_per_pool[i]))
            else
                __verbose__ && print(join(["Failed to start pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"],". Verify that you have correct credentials, a batch account with AAD authentication and that pool parameters are correct.\n"]))
            end
        end
    end
end

"""
    `create_pool_and_resource_file(startup_script; enable_auto_scale=false, auto_scale_formula=nothing,
    auto_scale_evaluation_interval_minutes=nothing, image_resource_id=nothing)`

 Create a batch pool using a given startup script that will be executed by nodes joining the pool. The pool parameters are specified
 in a separate JSON file. E.g. set `ENV["PARAMETERS"] = "parameters.json"` to point to the parameter file.

 *Required arguments*:

 - `startup_script` (String): Path to a bash script that is executed by nodes joining the batch pool.

 *Optional arguments*:

 - `enable_auto_scale` (Bool): Enable auto-scaling of the pool. Requires the `auto_scale_formula` and
    `auto_scale_evaluation_interval_minutes` to be set. If the number of VMs in the pool cannot be specified.

 - `auto_scale_formula` (String): Formula for auto-scaling the pool.
    See "https://docs.microsoft.com/en-us/azure/batch/batch-automatic-scaling" for details.

 - `image_resource_id` (String): Image ID of the VM image.

 *Output*:

 - `Nothing`

 See also:  [`create_pool`](@ref)
 """
function create_pool_and_resource_file(startup_script; enable_auto_scale=false, auto_scale_formula=nothing,
    auto_scale_evaluation_interval_minutes=nothing, image_resource_id=nothing)

    # Create container if it doesn't exist
    for client in __clients__
        create_blob_containers(client["blob_client"], [__container__])
    end
    if enable_auto_scale
        num_nodes = nothing
    else
        num_nodes = __params__["_NODE_COUNT_PER_POOL"]
    end

    # Enable inter node connection?
    enable_inter_node = parse(Bool, __params__["_INTER_NODE_CONNECTION"])

    # Container?
    if __params__["_CONTAINER"] == "None"
        docker_container = nothing
    else
        docker_container = __params__["_CONTAINER"]
    end

    # Divide clients among pools
    if ~isnothing(__credentials__)
        credential_per_pool, clients_per_pool, resources_per_pool = divide_clients_among_pools()
    else
        __verbose__ && @warn "Could not create pool(s). No credentials were found."
        return nothing
    end
    num_pools = parse(Int, __params__["_POOL_COUNT"])

    for i=1:num_pools
        pool_id = join([__params__["_POOL_ID"], "_", i])
        try
            azureclusterlesshpc.create_pool_and_resource_file(clients_per_pool[i], pool_id, __params__["_POOL_VM_SIZE"], num_nodes,
                __params__["_NODE_OS_PUBLISHER"], __params__["_NODE_OS_OFFER"], __params__["_NODE_OS_SKU"], startup_script, __container__;
                enable_inter_node=enable_inter_node, enable_auto_scale=enable_auto_scale, auto_scale_formula=auto_scale_formula,
                auto_scale_evaluation_interval_minutes=auto_scale_evaluation_interval_minutes, image_resource_id=image_resource_id,
                container=docker_container)

            # Keep track of active pools
            __verbose__ && print(join(["Created pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"], " with ", num_nodes, " nodes.\n"]))
            push!(__active_pools__, Dict("pool_id" => pool_id, "clients" => clients_per_pool[i], "credentials" => credential_per_pool[i],
                "resources" => resources_per_pool[i]))
        catch e
            if typeof(e) == PyCall.PyError && e.val.error.code == "PoolExists"
                __verbose__ && print(join(["Pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"]," already exists.\n"]))
                push!(__active_pools__, Dict("pool_id" => pool_id, "clients" => clients_per_pool[i], "credentials" => credential_per_pool[i],
                    "resources" => resources_per_pool[i]))
            else
                __verbose__ && print(join(["Failed to start pool ", i ," of ", num_pools, " in ", credential_per_pool[i]["_REGION"],". Verify that you have correct credentials, a batch account with AAD authentication and that pool parameters are correct.\n"]))
            end
        end
    end
end


text:
### Instability Plots
# Generates the code for the intability/suboptimality in the paper
# Presentation uses similar plots but with different paramters, see presentation plots.jl

using Distributions, Random, DelimitedFiles
include("../src/JS_SAA_main.jl")

function genSinglePathCurves(subcase; seed=8675409, usePoisson=true)
	K = 10000
	N = 10

	if subcase == :TradeoffBad
		p0 = .3
		s = .5
		path = string("../Results/TradeoffBad", usePoisson, ".csv")
	elseif subcase == :TradeoffGoodP0
		p0 = .75
		s = .5
		path = string("../Results/TradeoffGoodP0", usePoisson, ".csv")
	elseif subcase == :TradeoffGoodS
		p0 = .3
		s = .2
		path = string("../Results/TradeoffGoodS", usePoisson, ".csv")
	end
	Random.seed!(seed)

	#Gen Pks uniformly over interval for now
	ps = rand(K)
	ps = .6 .+ .3 * rand(K)

	f = open(path, "w")

	if usePoisson
		Nhats = rand(Poisson(N), K)
		Nhats[Nhats .== 0] .= 1
	else
		Nhats = fill(N, K)
	end
	mhats = JS.nv_sim_path(ps, Nhats);

	alpha_grid = range(0, stop=2N, length=100)

	#note, this actually computes N * ZLoo
	#loo_perf = map(a->JS.nv_loo(mhats, p0, a, Nhats, s), alpha_grid)
	alphaOR, jstar, outOR = JS.nv_oracle_alpha(mhats, ps, p0, alpha_grid, Nhats, s)
	alphaLOO, jstar, loo_perf = JS.nv_loo_alpha(mhats, p0, alpha_grid, Nhats, s)
	loo_perf /= N

	saa_subopt = map(a-> JS.nv_saa(mhats, p0, a, Nhats, s), alpha_grid) ./ N
	saa_zero = saa_subopt[1]
	saa_subopt = saa_subopt .- saa_zero
	instab = loo_perf - saa_subopt .- saa_zero;

	#Header
	writedlm(f, ["Alpha" "OR" "LOO" "SAA-SubOpt" "Instability"], ',')
	writedlm(f, [vec(alpha_grid) outOR loo_perf saa_subopt instab], ',')
	close(f)
end

genSinglePathCurves(:TradeoffBad)
genSinglePathCurves(:TradeoffGoodP0)
genSinglePathCurves(:TradeoffGoodS)

genSinglePathCurves(:TradeoffBad, usePoisson=false)
genSinglePathCurves(:TradeoffGoodP0, usePoisson=false)
genSinglePathCurves(:TradeoffGoodS, usePoisson=false)

ext: jl


text:
export grid_sample, ∇grid_sample

@inline in_bounds(h, w, H, W) = 1 ≤ h ≤ H && 1 ≤ w ≤ W
# Borders are considered out-of-bounds for gradient.
@inline clip_coordinate(coordinate, dim_size) = min(dim_size, max(1, coordinate))
@inline function ∇clip_coordinate(coordinate::C, dim_size) where C
    if coordinate ≤ 1
        return C(1), C(0)
    elseif coordinate ≥ dim_size
        return C(dim_size), C(0)
    end
    coordinate, C(1)
end

@inline unnormalize(coordinate, dim_size) = ((coordinate + 1.0) * 0.5) * (dim_size - 1.0) + 1.0
@inline ∇unnormalize(coordinate, dim_size) = unnormalize(coordinate, dim_size), (dim_size - 1.0) * 0.5

@inline compute_source_index(coordinate, dim_size, ::Val{:zeros}) = unnormalize(coordinate, dim_size)
@inline compute_source_index(coordinate, dim_size, ::Val{:border}) = clip_coordinate(unnormalize(coordinate, dim_size), dim_size)

@inline ∇compute_source_index(coordinate, dim_size, ::Val{:zeros}) = ∇unnormalize(coordinate, dim_size)
@inline function ∇compute_source_index(coordinate, dim_size, ::Val{:border})
    source_coordinate, grad_in = ∇unnormalize(coordinate, dim_size)
    source_coordinate, grad_clip = ∇clip_coordinate(source_coordinate, dim_size)
    source_coordinate, grad_in * grad_clip
end

"""
    grid_sample(input::AbstractArray{T, 4}, grid::AbstractArray{T, 4}; padding_mode = :zeros)

Given `input`, compute output by sampling `input` values at pixel
locations from `grid`. Uses bilinear interpolation to calculate output values.

This implementation assumes the extrema (`-1` and `1`) are considered
as referring to the center points of the input’s corner pixels
(i.e. align corners is `true`).

# Arguments

- `input`: Input array in `(W_in, H_in, C, N)` shape.
- `grid`: Input grid in `(2, W_out, H_out, N)` shape.
    Where for each `(W_out, H_out, N)` grid contains `(x, y)`
    coordinates that specify sampling locations normalized by the `input` shape.

    Therefore, `x` and `y` should have values in `[-1, 1]` range.
    For example, `(x = -1, y = -1)` is the left-top pixel of `input`,
    and `(x = 1, y = 1)` is the right-bottom pixel of `input`.

    Out-of-bound values are handled according to the `padding_mode`.
- `padding_mode`: Out-of-bound padding.
    `:zeros` to use `0` for out-of-bound grid locations.
    `:border` to use border values for out-of-bound grid locations.
    Default is `:zeros`.

# Returns

`(W_out, H_out, C, N)` sampled grid from `input`.

# Examples

In the example below, grid contains two out-of-bound sampling locations,
which are handled differently, depending on the `padding_mode`.

```jldoctest
julia> x = reshape(collect(1.0:4.0), (2, 2, 1, 1))
2×2×1×1 Array{Float64, 4}:
[:, :, 1, 1] =
 1.0  3.0
 2.0  4.0

julia> grid = Array{Float64}(undef, 2, 3, 2, 1);

julia> grid[:, 1, 1, 1] .= (-3, -1);

julia> grid[:, 2, 1, 1] .= (0, -1);

julia> grid[:, 3, 1, 1] .= (1, -1);

julia> grid[:, 1, 2, 1] .= (-1, 1);

julia> grid[:, 2, 2, 1] .= (0, 1);

julia> grid[:, 3, 2, 1] .= (3, 1);

julia> grid_sample(x, grid; padding_mode=:zeros)
3×2×1×1 Array{Float64, 4}:
[:, :, 1, 1] =
 0.0  3.0
 1.5  3.5
 2.0  0.0

julia> grid_sample(x, grid; padding_mode=:border)
3×2×1×1 Array{Float64, 4}:
[:, :, 1, 1] =
 1.0  3.0
 1.5  3.5
 2.0  4.0
```
"""
function grid_sample(input::AbstractArray{T, 4}, grid; padding_mode = :zeros) where T
    _, _, iC, iN = size(input)
    _, gW, gH, _ = size(grid)
    output = similar(input, T, (gW, gH, iC, iN))
    grid_sample!(output, input, grid, padding_mode)
end
function grid_sample!(output, input, grid, padding_mode)
    pad = Val(padding_mode)
    iW, iH, iC, iN = size(input)
    _, gW, gH, _ = size(grid)
    # Loop over each output pixel.
    Threads.@threads for n in 1:iN
        for w in 1:gW, h in 1:gH
            _grid_sample_kernel!(output, input, grid, pad, w, h, n, iW, iH, iC)
        end
    end
    output
end
@inline function _grid_sample_kernel!(
    output, input, grid, padding_mode, w, h, n, iW, iH, iC,
)
    # Get the corresponding (x, y) coordinates from the grid.
    @inbounds x, y = grid[1, w, h, n], grid[2, w, h, n]
    ix = compute_source_index(x, iW, padding_mode)
    iy = compute_source_index(y, iH, padding_mode)
    # Get corner pixel values from (ix, iy) in north-east-south-west directions.
    ix_nw, iy_nw = floor(Int, ix), floor(Int, iy)
    ix_ne, iy_ne = ix_nw + 1, iy_nw
    ix_sw, iy_sw = ix_nw, iy_nw + 1
    ix_se, iy_se = ix_ne, iy_sw
    # Get surfaces to each neighbor (a.k.a. interpolation weights).
    nw = (ix_se - ix) * (iy_se - iy)
    ne = (ix - ix_sw) * (iy_sw - iy)
    sw = (ix_ne - ix) * (iy - iy_ne)
    se = (ix - ix_nw) * (iy - iy_nw)
    # ∀ channel: Calculate bilinear weighted pixel value.
    @inbounds for c in 1:iC
        r = 0.0
        if in_bounds(iy_nw, ix_nw, iH, iW)
            r += input[ix_nw, iy_nw, c, n] * nw
        end
        if in_bounds(iy_ne, ix_ne, iH, iW)
            r += input[ix_ne, iy_ne, c, n] * ne
        end
        if in_bounds(iy_sw, ix_sw, iH, iW)
            r += input[ix_sw, iy_sw, c, n] * sw
        end
        if in_bounds(iy_se, ix_se, iH, iW)
            r += input[ix_se, iy_se, c, n] * se
        end
        output[w, h, c, n] = r
    end
end

"""
    ∇grid_sample(Δ::AbstractArray{T, 4}, input::AbstractArray{T, 4}, grid::AbstractArray{T, 4}; padding_mode = :zeros) where T

# Arguments

- `Δ`: Input gradient in `(W_out, H_out, C, N)` shape
    (same as output of the primal computation).
- `input`: Input from primal computation in `(W_in, H_in, C, N)` shape.
- `grid`: Grid from primal computation in `(2, W_out, H_out, N)` shape.
- `padding_mode`: Out-of-bound padding.
    `:zeros` to use `0` for out-of-bound grid locations.
    `:border` to use border values for out-of-bound grid locations.
    Should be the same as in primal computation.
    Default is `:zeros`.

# Returns

`dinput` (same shape as `input`) and `dgrid` (same shape as `grid`) gradients.
"""
function ∇grid_sample(Δ::AbstractArray{T, 4}, input::AbstractArray{T, 4}, grid; padding_mode = :zeros) where T
    dx = zeros(T, size(input))
    dgrid = similar(grid)
    ∇grid_sample!(dx, dgrid, Δ, input, grid, padding_mode)
end
function ∇grid_sample!(dx, dgrid, Δ, input, grid, padding_mode)
    pad = Val(padding_mode)
    iW, iH, iC, iN = size(input)
    gW, gH = size(grid, 2), size(grid, 3)
    # Loop over each output pixel.
    Threads.@threads for n in 1:iN
        for w in 1:gW, h in 1:gH
            _∇grid_sample_kernel!(dx, dgrid, Δ, input, grid, pad, w, h, n, iW, iH, iC)
        end
    end
    dx, dgrid
end
@inline function _∇grid_sample_kernel!(
    dx, dgrid, Δ, input, grid, padding_mode, w, h, n, iW, iH, iC,
)
    # Get corresponding (x, y) from grid.
    @inbounds x, y = grid[1, w, h, n], grid[2, w, h, n]
    # Compute multipliers for gradinets on ix, iy.
    ix, gix_mult = ∇compute_source_index(x, iW, padding_mode)
    iy, giy_mult = ∇compute_source_index(y, iH, padding_mode)
    # Get corner pixel values from (ix, iy) in north-east-south-west directions.
    ix_nw, iy_nw = floor(Int, ix), floor(Int, iy)
    ix_ne, iy_ne = ix_nw + 1, iy_nw
    ix_sw, iy_sw = ix_nw, iy_nw + 1
    ix_se, iy_se = ix_ne, iy_sw
    # Get surfaces to each neighbor (a.k.a. interpolation weights).
    nw = (ix_se - ix) * (iy_se - iy)
    ne = (ix - ix_sw) * (iy_sw - iy)
    sw = (ix_ne - ix) * (iy - iy_ne)
    se = (ix - ix_nw) * (iy - iy_nw)
    # ∀ channel: Calculate billinear weighted pixel value.
    gix, giy = 0.0, 0.0
    @inbounds for c in 1:iC
        g_out = Δ[w, h, c, n]
        # Calculate dx and dgrid partials.
        if in_bounds(iy_nw, ix_nw, iH, iW)
            _safe_add!(dx, g_out * nw, ix_nw, iy_nw, c, n)
            nw_val = input[ix_nw, iy_nw, c, n]
            gix -= nw_val * (iy_se - iy) * g_out
            giy -= nw_val * (ix_se - ix) * g_out
        end
        if in_bounds(iy_ne, ix_ne, iH, iW)
            _safe_add!(dx, g_out * ne, ix_ne, iy_ne, c, n)
            ne_val = input[ix_ne, iy_ne, c, n]
            gix += ne_val * (iy_sw - iy) * g_out
            giy -= ne_val * (ix - ix_sw) * g_out
        end
        if in_bounds(iy_sw, ix_sw, iH, iW)
            _safe_add!(dx, g_out * sw, ix_sw, iy_sw, c, n)
            sw_val = input[ix_sw, iy_sw, c, n]
            gix -= sw_val * (iy - iy_ne) * g_out
            giy += sw_val * (ix_ne - ix) * g_out
        end
        if in_bounds(iy_se, ix_se, iH, iW)
            _safe_add!(dx, g_out * se, ix_se, iy_se, c, n)
            se_val = input[ix_se, iy_se, c, n]
            gix += se_val * (iy - iy_nw) * g_out
            giy += se_val * (ix - ix_nw) * g_out
        end
    end
    @inbounds dgrid[1, w, h, n] = gix_mult * gix
    @inbounds dgrid[2, w, h, n] = giy_mult * giy
end

@inline function _safe_add!(dx, value, ix, iy, c, n)
    @inbounds dx[ix, iy, c, n] += value
end

function rrule(::typeof(grid_sample), x, grid; padding_mode)
    y = grid_sample(x, grid; padding_mode=padding_mode)
    function grid_sample_pullback(Δ)
        ∇x, ∇grid = ∇grid_sample(unthunk(Δ), x, grid; padding_mode=padding_mode)
        NoTangent(), ∇x, ∇grid
    end
    return y, grid_sample_pullback
end


text:
using Measurements

@testset "Measurements" begin
    x = rand(2)
    u = [0.5;0.5]
    H = Stage(ts,measurement.(x,u))
    @test Measurements.uncertainty.(quantity(H)) == [0.5;0.5]
end


text:
# basic figure with julia

@time using Plots  # 6.9s
@time pyplot()  # 5.2s

t = collect(0.0:0.01:2.0)
s = sin.(2*pi*t)
c = cos.(2*pi*t)

@time f = plot(t, s, label="sin") # 11s
plot!(t, c, label="cos")
plot!(xlab="time (s)", ylab="voltage (mV)")
plot!(legend=:bottomleft)

@time savefig("test.png") # 5.8s

display(f)

println("done. press <ENTER> to quit")
readline()

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
module FASTQDemultiplexerKallisto

import FASTQDemultiplexer:
    Output, Protocol, InterpretedRecord,
    mergeoutput
using DataFrames
using Weave

include("kallisto.jl")

abstract OutputKallisto <: Output

type OutputKallistoSingle <: OutputKallisto
    kallisto::Kallisto{1}
    results::DataFrame
end

type OutputKallistoDouble <: OutputKallisto
    kallisto::Kallisto{1}
    kallisto2::Kallisto{2}
    results::DataFrame
end

issingle(::OutputKallistoSingle) = true
issingle(::OutputKallistoDouble) = false

function OutputKallisto(protocol::Protocol;
                        index::String = "",
                        index2::String = "",
                        kwargs...)

    results = DataFrame(cellid = UInt[],
                        umiid = UInt[],
                        groupname = PooledDataArray(String,UInt8,0),
                        alignment = Int[])

    k = Kallisto{1}(index)

    if !isempty(index2)
        k2 = Kallisto{2}(index2)
        return OutputKallistoDouble(
            k,k2,vcat(results,DataFrame(alignment2 = Int[])))
    else
        return OutputKallistoSingle(k,results)
    end
end

function Base.write(o::OutputKallisto,
                    ir::InterpretedRecord)
    if ir.groupid >= 0
        alignment = align(o.kallisto,ir)
        if issingle(o)
            push!(o.results,(ir.cellid,ir.umiid,ir.groupname,alignment))
        else
            alignment2 = align(o.kallisto2,ir)
            push!(o.results,(ir.cellid,ir.umiid,ir.groupname,alignment,alignment2))
        end
    end
end

function mergeoutput{OK<:OutputKallisto}(outputs::Vector{OK};
                                         outputdir::String=".",
                                         writealignment = true,
                                         writereport = true,
                                         index::String = "",
                                         index2::String = "",
                                         name::String = basename(index),
                                         name2::String = basename(index2),
                                         # subsamplecells::Int = 100_000,
                                         # subsamplereads::Int = 100_000,
                                         plotsbackend::String = "plotly",
                                         kwargs...)
    if length(outputs) > 1
        results = vcat((o.results for o in outputs)...)
    else
        results = outputs[1].results
    end

    if !(plotsbackend in ["plotly", "gr"])
        error("Unknown backend, try `plotly` or `gr`")
    end

    if writealignment
        mkpath(outputdir)
        writetable(joinpath(outputdir,"alignment.tsv"),results)
    end

    if writereport
        mkpath(outputdir)
        report(results, outputdir, [name,name2], plotsbackend)
    end
end

function report(results, outputdir, names, backend)
    args = Dict{Symbol,Any}()
    args[:results] = results
    args[:indexnames] = names
    args[:backend] = backend

    weave(Pkg.dir("FASTQDemultiplexerKallisto",
                  "src","weave","alignment.jmd"),
          args = args,
          out_path = outputdir)
end

end


text:
using BitCircuits
using Base.Test

or(left, right) = left | right
and(left, right) = left & right

# Demo using the following equivalency:
# p <- ab + ac
# q <- a (b + c)
# p = q

a = Variable(0)
b = Variable(1)
c = Variable(2)

ab = Operation(and, a, b)
ac = Operation(and, a, c)
p = Operation(or, ab, ac)
boc = Operation(or, b, c)
q = Operation(and, a, boc)

@test equal(p, q)
@test evaluate(ab, 0b000011) == true
@test evaluate(ab, 0b000010) == false
@test evaluate(ab, 0b000001) == false
@test evaluate(ab, 0b000000) == false


text:
import Base: show
import Oceananigans: short_show

location_str(::Type{Face})    = "Face"
location_str(::Type{Cell})    = "Cell"
location_str(::Type{Nothing}) = "⋅"

show_location(X, Y, Z) = "($(location_str(X)), $(location_str(Y)), $(location_str(Z)))"

show_location(field::AbstractField{X, Y, Z}) where {X, Y, Z} = show_location(X, Y, Z)

short_show(m::Missing) = "$m"

short_show(field::Field) = string("Field located at ", show_location(field))

short_show(field::AveragedField) = string("AveragedField over dims=$(field.dims) located at ", show_location(field), " of ", short_show(field.operand))
short_show(field::ComputedField) = string("ComputedField located at ", show_location(field), " of ", short_show(field.operand))

show(io::IO, field::Field) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))\n",
          "└── boundary conditions: $(short_show(field.boundary_conditions))")

show_status(::Nothing) = "nothing"
show_status(status) = "time=$(status.time)"

show(io::IO, field::AveragedField) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))", '\n',
          "├── dims: $(field.dims)", '\n',
          "├── operand: $(short_show(field.operand))", '\n',
          "└── status: ", show_status(field.status), '\n')

show(io::IO, field::ComputedField) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))", '\n',
          "├── operand: $(short_show(field.operand))", '\n',
          "└── status: ", show_status(field.status), '\n')

short_show(array::OffsetArray{T, D, A}) where {T, D, A} = string("OffsetArray{$T, $D, $A}")


text:
#=
Given a multiset of integers, return whether it can be partitioned into two subsets whose sums are the same.

For example, given the multiset {15, 5, 20, 10, 35, 15, 10}, it would return true, since we can split it up into {15, 5, 10, 15, 10} and {20, 35}, which both add up to 55.

Given the multiset {15, 5, 20, 10, 35}, it would return false, since we can't split it up into two subsets that add up to the same sum.
=#

function ispartitionable(multiset::Array{Int,1})
    sorted_multiset = sort(multiset, rev = true)
    left_subset = sorted_multiset[1:1]
    right_subset = sorted_multiset[2:end]
    index = 1
    while sum(left_subset) < sum(right_subset) && index < length(right_subset)
        if sum(left_subset) + right_subset[index] <= sum(right_subset) - right_subset[index]
            push!(left_subset, popfirst!(right_subset))
            index = 1
        else
            index = index + 1
        end
    end

    if sum(left_subset) == sum(right_subset)
        return true
    else
        return false
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Test matrices for regularization methods from Hansen's
# Regularization toolbox

"""
Oscillating Matrix
==================
A matrix `A` is called oscillating if `A` is totally
    nonnegative and if there exists an integer `q > 0` such that
    `A^q` is totally positive.

*Input options:*

+ [type,] Σ: the singular vaule spectrum of the matrix.

+ [type,] dim, mode: `dim` is the dimension of the matrix.
        `mode = 1`: geometrically distributed singular values.
        `mode = 2`: arithmetrically distributed singular values.

+ [type,] dim: `mode = 1`.

*Groups:* ['symmetric','posdef', 'random', 'eigen']

*References:*

**Per Christian Hansen**, Test matrices for
    regularization methods. SIAM J. SCI. COMPUT Vol 16,
    No2, pp 506-512 (1995).
"""
function oscillate(Σ::Vector{T}) where T
    n = length(Σ)
    dv = rand(T, 1, n)[:] .+ eps(T)
    ev = rand(T, 1, n-1)[:] .+ eps(T)
    B = Bidiagonal(dv, ev, :U)
    U, S, V = svd(B)
    return U*Diagonal(Σ)*U'
end
function oscillate(::Type{T}, n::Integer, mode::Integer) where T
    κ = sqrt(1/eps(T))
    if mode == 1
        factor = κ^(-1/(n-1))
        Σ = factor.^[0:n-1;]
    elseif mode == 2
        Σ = ones(T, n) - T[0:n-1;]/(n-1)*(1 - 1/κ)
    else
        throw(ArgumentError("invalid mode value."))
    end
    return oscillate(Σ)
end
oscillate(::Type{T}, n::Integer) where T = oscillate(T, n, 2)
oscillate(args...) = oscillate(Float64, args...)
oscillate(::Type, args...) = throw(MethodError(oscillate, Tuple(args)))

struct RegProb{T}
    A::AbstractMatrix{T}  # matrix of interest
    b::AbstractVector{T}  # right-hand side
    x::AbstractVector{T}  # the solution to Ax = b
end

struct RegProbNoSolution{T}
    A::AbstractMatrix{T}  # matrix of interest
    b::AbstractVector{T}  # right-hand side
end

function show(io::IO, p::RegProb)
    println(io, "Test problems for Regularization Methods")
    println(io, "A:")
    display(p.A)
    println(io, "b:")
    display(p.b)
    println(io, "x:")
    display(p.x)
end

function show(io::IO, p::RegProbNoSolution)
    println(io, "Test problems for Regularization Methods with No Solution")
    println(io, "A:")
    display(p.A)
    println(io, "b:")
    display(p.b)
end

# The following test problems are derived from Per Christian Hansen's
# Regularization tools for MATLAB.
# http://www.imm.dtu.dk/~pcha/Regutools/
#
# BSD License
# Copyright (c) 2015, Per Christian Hansen
# All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in
#       the documentation and/or other materials provided with the distribution
#     * Neither the name of the DTU Compute nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

"""
Computation of the Second Derivative
====================================
A classical test problem for regularization algorithms:
    This is a mildly ill-posed problem. It is a discretization
    of a first kind Fredholm integral equation whose kernel K
    is the Green's function for the second derivative.

*Input options:*

+ [type,] dim, example, [matrixonly]: the dimension of the
        matrix is `dim`.  One choose between between the following right-hand
        g and solution f:

        example = 1 gives g(s) = (s^3 - s)/6, f(t) = t.

        example = 2 gives g(s) = exp(s) + (1 -e)s - 1, f(t) = exp(t)

        example = 3 gives  g(s) = | (4s^3 - 3s)/24,  s < 0.5
                                  | (-4s^3 + 12s^2 - 9s + 1)/24, s>= 0.5
                           f(t) = | t, t < 0.5
                           f(t) = | 1- t, t >= 0.5.

        If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be
        generated (`matrixonly = true` by default).

+ [type,] dim, [matrixonly]: `example = 1`.

*Groups:* ["regprob"]

*References:*

**P. C. Hansen**, Regularization tools: A MATLAB pacakge for
    analysis and solution of discrete ill-posed problems.
    Numerical Algorithms, 6(1994), pp.1-35
"""
function deriv2(::Type{T}, n::Integer, example::Integer, matrixonly::Bool = true) where T
    h = convert(T, one(T)/n); sqh = sqrt(h)
    h32 = h*sqh; h2 = h^2; sqhi = one(T)/sqh
    t = 2/3; A = zeros(T, n, n)

    # compute A
    for i = 1:n
        A[i,i] = h2*((i^2 - i + 0.25)*h - (i - t))
        for j = 1:i-1
            A[i,j] = h2*(j - 0.5)*((i - 0.5)*h - 1)
        end
    end
    A = A + tril(A, -1)'

    if matrixonly
        return A
    else
        b = zeros(T, n)
        x = zeros(T, n)
        if example == 1
            # compute b
            [b[i] = h32*(i - 0.5)*((i^2 + (i-1)^2)*h2/2 - 1)/6 for i = 1:n]

            # compute x
            [x[i] = h32*(i - 0.5) for i = 1:n]
        elseif example == 2
            ee = one(T) - exp(one(T))
            [b[i] = sqhi*(exp(i*h) - exp((i-1)*h) + ee*(i-0.5)*h2 - h) for i = 1:n]
            [x[i] = sqhi*(exp(i*h) - exp((i-1)*h)) for i = 1:n]
        elseif example == 3
            mod(n, 2) == 0 || error("The order n must be even.")
            for i = 1:div(n,2)
                s12 = (i*h)^2; s22 = ((i-1)*h)^2
                b[i] = sqhi*(s12 + s22 - 1.5)*(s12 - s22)/24
            end
            for i = div(n, 2)+1:n
                s1 = i*h; s12 = s1^2; s2 = (i-1)*h; s22 = s2^2
                b[i] = sqhi*(-(s12 + s22)*(s12 - s22) + 4*(s1^3 - s2^3)
                       - 4.5*(s12 - s22) + h)/24
            end
            [x[i] = sqhi*((i*h)^2 - ((i-1)*h)^2)/2 for i = 1:div(n, 2)]
            [x[i] = sqhi*(h - ((i*h)^2 - ((i-1)*h)^2)/2) for i = div(n, 2)+1:n]
        else
            throw(ArgumentError("Illegal value of example."))
        end
        return RegProb(A, b, x)
    end
end
deriv2(::Type{T}, n::Integer, matrixonly::Bool = true) where T = deriv2(T, n, 1, matrixonly)
deriv2(args...) = deriv2(Float64, args...)
deriv2(::Type, args...) = throw(MethodError(deriv2, Tuple(args)))

"""
One-Dimensional Image Restoration Model
=======================================
This test problem uses a first-kind Fredholm integral equation
    to model a one-dimentional image restoration situation.

*Input options:*

+ [type,] dim, [matrixonly]: the dimesion of the matrix `dim` must be even.
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**C. B. Shaw, Jr.**, Improvements of the resolution of
    an instrument by numerical solution of an integral equation.
    J. Math. Ana. Appl. 37 (1972), 83-112.
"""
function shaw(::Type{T}, n::Integer, matrixonly::Bool = true) where T
    mod(n, 2) == 0 || error("The dimension of the matrix must be even.")
    h = π/n; A = zeros(T, n, n)

    # compute A
    co = cos.(-π/2 .+ T[.5:n-.5;] .* h)
    psi = π .* sin.(-π/2 .+ T[.5:n-.5;] .* h)
    for i in 1:div(n,2)
        for j in i:n-i
            ss = psi[i] +psi[j]
            A[i,j] = ((co[i] + co[j])*sin(ss)/ss)^2
            A[n-j+1, n-i+1] = A[i,j]
        end
        A[i, n-i+1] = (2*co[i])^2
    end
    A = A + triu(A, 1)'; A = A*h

    if matrixonly
        return A
    else
        # compute x and b
        a1 = 2
        c1 = 6
        t1 = .8
        a2 = 1
        c2 = 2
        t2 = -.5
        x = a1 .* exp.(-c1 .* (-π/2 .+ T[.5:n-.5;] .* h .- t1).^2) .+
            a2 .* exp.(-c2 .* (-π/2 .+ T[.5:n-.5;] .* h .- t2).^2)
        b = A*x
        return RegProb(A, b, x)
    end
end
shaw(args...) = shaw(Float64, args...)
shaw(::Type, args...) = throw(MethodError(shaw, Tuple(args)))

"""
A Problem with a Discontinuous Solution
=======================================

*Input options:*

+ [type,] dim, t1, t2, [matrixonly]: the dimension of matrix is `dim`.
    `t1` and `t2` are two real scalars such that `0 < t1 < t2 < 1`.
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

+ [type,] n, [matrixonly]: `t1 = 1/3` and `t2 = 2/3`.

*Groups:* ["regprob"]

*References:*

**G. M. Wing**, A Primer on Integral Equations of the
    First Kind, Society for Industrial and Applied Mathematics, 1991, p. 109.
"""
function wing(::Type{T}, n::Integer, t1::Real, t2::Real, matrixonly = true) where T
    t1 < t2 || error("t1 must be smaller than t2")
    A = zeros(T, n, n)
    h = T(1/n)

    # compute A
    sti = (T[1:n;] .- 0.5) * h
    for i in 1:n
        A[i,:] .= h .* sti .* exp.(-sti[i] .* sti.^2)
    end

    if matrixonly
        return A
    else
        # compute b
        b = sqrt(h) .* 0.5 .* (exp.(-sti .* t1^2) .- exp.(-sti .* t2^2))./sti

        # compute x
        indices = [findfirst(t1 .< sti):findlast(t2 .> sti);]
        x = zeros(T,n); x[indices] = sqrt(h)*ones(length(indices))
        return RegProb(A, b, x)
    end
end
wing(::Type{T}, n::Integer, matrixonly = true) where T = wing(T, n, 1/3, 2/3, matrixonly)
wing(args...) = wing(Float64, args...)
wing(::Type, args...) = throw(MethodError(wing, Tuple(args)))

"""
Severely Ill-posed Problem Suggested by Fox & Goodwin
=====================================================
This is a model problem discretized by simple quadrature, which does
not satifiy the discrete Picard condition for the small singular values.

*Input options:*

+ [type,] dim, [matrixonly]: `dim` is the dimension of the matrix.
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**C. T. H. Baker**, The Numerical Treatment of Integral
    Equations, Clarendon Press, Oxford, 1977, p. 665.
"""
function foxgood(::Type{T}, n::Integer, matrixonly = true) where T
    h = T(1/n)
    t = h*(T[1:n;] .- one(T)/2)

    A = h*sqrt.((t.^2)*ones(T,n)' + ones(T, n) * (t.^2)')

    if matrixonly
        return A
    else
        x = t
        b = zeros(T, n)
        for i in 1:n
            b[i] = ((one(T) + t[i]^2)^T(1.5) - t[i]^3)/3
        end
        return RegProb(A, b, x)
    end
end
foxgood(args...) = foxgood(Float64, args...)
foxgood(::Type, args...) = throw(MethodError(foxgood, Tuple(args)))

"""
Inverse Heat Equation
=====================

*Input options:*

+ [type,] dim, κ, [matrixonly]: `dim` is the dimension of the matrix and
    `dim` must be even. `κ` controls the ill-conditioning of the matrix.
    (`κ = 5` gives a well-conditioned problem and `κ = 1`
     gives an ill conditoned problem).
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

+ [type,] n, [matrixonly]: `κ = 1`.

*Groups:* ["regprob"]

*References:*

**A. S. Carasso**, Determining surface temperatures
    from interior observations, SIAM J. Appl. Math. 42 (1982), 558-574.
"""
function heat(::Type{T}, n::Integer, κ::Real, matrixonly::Bool = true) where T
    mod(n, 2) == 0 || error("The dimension of the matrix must be even.")
    h = one(T)/n; t = T[h/2:h:1;]
    c = h/(2*κ*sqrt(π))
    d = one(T)/(4*κ^2)

    # compute the matrix A
    m = length(t); k = zeros(T, m)
    [k[i] = c*t[i]^(-1.5)*exp(-d/t[i]) for i in 1:m]
    r = zeros(T, m); r[1] = k[1]
    A = toeplitz(T, k, r)

    if matrixonly
        return A
    else
        # compute the vectors x and b
        x = zeros(T, n)
        for i = 1:div(n,2)
            ti = i*20/n
            if ti < 2
                x[i] = 0.75*ti^2/4
            elseif ti < 3
                x[i] = 0.75 + (ti - 2)*(3 - ti)
            else
                x[i] = 0.75*exp(-(ti - 3)*2)
            end
        end
        x[div(n,2)+1:n] = zeros(T, div(n, 2))
        b = A*x
        return RegProb(A, b, x)
    end
end
heat(::Type{T}, n::Integer, matrixonly::Bool = true) where T = heat(T, n, 1, matrixonly)
heat(args...) = heat(Float64, args...)
heat(::Type, args...) = throw(MethodError(heat, Tuple(args)))

"""
Fredholm Integral Equation of the First Kind
============================================

*Input options:*

+ [type,] dim, [matrixonly]: the dimenstion of the matrix is `dim`.
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**M. L. Baart**, The use of auto-correlation for
    pesudo-rank determination in noisy ill-conditioned linear-squares
    problems, IMA, J. Numer. Anal. 2 (1982), 241-247.
"""
function baart(::Type{T}, n::Integer, matrixonly::Bool = true) where T
    mod(n, 2) == 0 || error("The dimension of the matrix must be even.")
    hs  = T(π/(2*n))
    ht  = T(π/n)
    c   = one(T)/(3*sqrt(2))
    A   = zeros(T, n, n)
    ihs = T[0:n;]*hs
    n1  = n + 1
    nh  = div(n,2)
    f3  = exp.(ihs[2:n1]) .- exp.(ihs[1:n])

    # compute A
    for j = 1:n
        f1  = f3
        co2 = cos((j - one(T)/2)*ht)
        co3 = cos(j*ht)
        f2 = (exp.(ihs[2:n1].*co2) .- exp.(ihs[1:n].*co2))./co2
        if j == nh
            f3 = hs*ones(T, n)
        else
            f3 = (exp.(ihs[2:n1].*co3) .- exp.(ihs[1:n].*co3))./co3
        end
        A[:,j] .= c.*(f1 .+ 4 .* f2 .+ f3)
    end

    if matrixonly
        return A
    else
        # compute vector b
        si      = T[.5:.5:n;] .* hs
        si      = sinh.(si) ./ si
        b       = zeros(T, n)
        b[1]    = 1 + 4 * si[1] + si[2]
        b[2:n] .= si[2:2:2*n-2] .+ 4 .* si[3:2:2*n-1] .+ si[4:2:2*n]
        b      .= b .* sqrt(hs) ./ 3

        # compute vector x
        x = -diff(cos.(T[0:n;] .* ht)) ./ sqrt(ht)
        return RegProb(A, b, x)
    end
end
baart(args...) = baart(Float64, args...)
baart(::Type, args...) = throw(MethodError(baart, Tuple(args)))

"""
Phillips's \"Famous\" Problem
=============================

*Input options:*

+ [type,] dim, [matrixonly]: the dimenstion of the matrix is `dim`.
    If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**D. L. Phillips**, A technique for the numerical
    solution of certain integral equations of the first kind, J. ACM
    9 (1962), 84-97.
"""
function phillips(::Type{T}, n::Integer, matrixonly::Bool = true) where T
    mod(n, 4) == 0 || error("The dimension of the matrix must be a multiple of 4.")

    # compute A
    h = 12/n
    n4 = div(n, 4)
    r1 = zeros(T,n)
    c = cos.(T[-1:n4; ] * 4 * π/n)
    for i in 1:n4
        r1[i] = h + 9 / (h * π^2) * (2 * c[i + 1] - c[i] - c[i + 2])
    end
    r1[n4 + 1] = h / 2 + 9 / (h * π^2) * (cos(4 * π / n) - 1)
    A = toeplitz(T, r1)

    if matrixonly
        return A
    else
        # compute the vector b
        b = zeros(T, n)
        c = π/3
        for i = (div(n,2) + 1):n
            t1 = -6 + i*h
            t2 = t1 - h
            b[i] = t1*(6-abs(t1)/2) +
                ((3 - abs(t1)/2)*sin(c*t1) - 2/c*(cos(c*t1) - one(T)))/c -
                t2*(6 - abs(t2)/2) -
                ((3 - abs(t2)/2)*sin(c*t2) - 2/c*(cos(c*t2) - one(T)))/c
            b[n - i + 1] = b[i]
        end
        b ./= sqrt(h)

        # compute x
        x = zeros(T, n)
        x[2n4+1:3n4] = (h .+ diff(sin.(T[0:h:(3 + 10 * eps(T));] * c)) / c) / sqrt(h)
        x[n4+1:2*n4] = x[3*n4:-1:2*n4+1]
        return RegProb(A, b, x)
    end
end
phillips(args...) = phillips(Float64, args...)
phillips(::Type, args...) = throw(MethodError(phillips, Tuple(args)))

# replicates the grid vectors xgv and ygv to produce a full grid.
function meshgrid(xgv, ygv)
    X = [i for j in ygv, i in xgv]
    Y = [j for j in ygv, i in xgv]
    return X, Y
end

# MATLAB rounding behavior
# This is equivalent to RoundNearestTiesAway
# and it can be used for both Julia v0.3, v0.4
function round_matlab(x::AbstractFloat)
    y = trunc(x)
    ifelse(x==y,y,trunc(2*x-y))
end
round_matlab(::Type{T}, x::AbstractFloat) where T<:Integer = trunc(T,round_matlab(x))

"""
One-dimensional Gravity Surverying Problem
==========================================
Discretization of a 1-D model problem in gravity surveying, in
    which a mass distribution f(t) is located at depth d, while the
    vertical component of the gravity field g(s) is measured at the
    surface.

*Input options:*

+ [type,] dim, example, a, b, d, [matrixonly]: `dim` is the dimension
        of the matrix. Three examples are implemented.

       (a) example = 1 gives f(t) = sin(π*t) + 0.5*sin(2*π*t).

       (b) example = 2 gives f(t) = piecewise linear function.

       (c) example = 3 gives f(t) = piecewise constant function.

       The t integration interval is fixed to [0, 1], while the s
       integration interval [a, b] can be specified by the user.
       The parameter d is the depth at which the magnetic deposit is
       located. The larger the d, the faster the decay of the singular values.

       If matrixonly = false, the matrix A and vectors b and x in the linear system Ax = b will be generated
       (matrixonly = true by default).

+ [type,] dim, example, [matrixonly]: `a = 0, b = 1, d = 0.25`;

+ [type,] dim, [matrixonly]: `example = 1, a = 0, b = 1, d = 0.25`.
*Groups:* ["regprob"]

*References:*

**G. M. Wing and J. D. Zahrt**, A Primer on Integral Equations of
the First Kind, Society for Industrial and Applied Mathematics, Philadelphia, 1991, p. 17.
"""
function gravity(::Type{T}, n::Integer, example::Integer,
                    a::Number, b::Number, d::Number, matrixonly::Bool = true) where T
    dt = one(T)/n
    a  = T(a)
    b  = T(b)
    d  = T(d)
    ds = (b - a)/n
    tv = dt .* (T[1:n;] .- one(T) ./ 2)
    sv = a .+ ds .* (T[1:n;] .- one(T) ./ 2)
    Tm, Sm = meshgrid(tv, sv)
    A = dt .* d .* ones(T, n, n) ./ (d^2 .+ (Sm .- Tm).^2).^T(3/2)
    if matrixonly
        return A
    else
        x = ones(T, n)
        nt = round_matlab(Int, n/3)
        nn = round_matlab(Int, n*7/8)
        if example == 1
            x .= sin.(π .* tv) .+ sin.(2 .* π .* tv) ./ 2
        elseif example == 2
            x[1:nt]        .= 2 ./ nt .* [1:nt;]
            x[(nt + 1):nn] .= ((2 .* nn .- nt) .- [(nt .+ 1):nn;]) ./ (nn .- nt)
            x[(nn + 1):n]  .= (n .- [(nn .+ 1):n;]) ./ (n .- nn)
        elseif example == 3
            x[1:nt] .= 2
        else
            error("Illegal value of example")
        end
        b = A*x
        return RegProb(A, b, x)
    end
end

gravity(::Type{T}, n::Integer, example::Integer, matrixonly::Bool = true) where T =
    gravity(T, n, example, 0, 1, 0.25, matrixonly)

gravity(::Type{T}, n::Integer, matrixonly::Bool = true) where T =
           gravity(T, n, 1, 0, 1, 0.25, matrixonly)
gravity(args...) = gravity(Float64, args...)
gravity(::Type, args...) = throw(MethodError(gravity, Tuple(args)))

"""
Image Deblurring Test Problem
=============================
The generated matrix A is an `n*n-by-n*n` sparse, symmetric,
           doubly block Toeplitz matrix that models blurring of an n-by-n
           image by a Gaussian point spread function.

*Input options:*

+ [type,] dim, band, σ, [matrixonly]: the dimension of the matrix
          is `dim^2`. `band` is the half-bandwidth, only matrix elements within
          a distance `band-1` from the diagonal are nonzero. `σ` controls the
          width of the Gaussin point spread function. The larger the `σ`, the
          wider the function and the more ill posed the problem.
          If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
              (`matrixonly = true` by default).

+ [type,] dim, [matrixonly]: `band = 3, σ = 0.7`.

*Groups:* ["regprob", "sparse"]
"""
function blur(::Type{T}, n::Integer, band::Integer, σ::Number,
                 matrixonly::Bool = true) where T
    σ = T(σ)
    z = [exp.(-(T[0:band-1;].^2) / (2 * σ^2)); zeros(T, n - band)]
    A = toeplitz(T, z)
    A = sparse(A)
    A = (1 / T(2 * π * σ^2)) * kron(A, A)
    if matrixonly
        return A
    else
        # start with an image of all zeros
        n2  = round_matlab(Int, n/2)
        n3  = round_matlab(Int, n/3)
        n6  = round_matlab(Int, n/6)
        n12 = round_matlab(Int, n/12)

        m = max(n, 2 * n6 + 1 + n2 + n12)
        x = zeros(T, m, m)

        # add a large ellipse
        Te = zeros(T, n6, n3)
        for i = 1:n6
            for j = 1:n3
                if (i / n6)^2 + (j / n3)^2 < 1
                    Te[i,j] = 1
                end
            end
        end
        Te = [reverse(Te, dims=2) Te;]
        Te = [reverse(Te, dims=1); Te]

        x[2 .+ [1:2n6;], n3 - 1 .+ [1:2n3;]] = Te

        # add a smaller ellipse
        Te = zeros(T, n6, n3)
        for i = 1:n6
            for j = 1:n3
                if (i / n6)^2 + (j / n3)^2 < 0.6
                    Te[i,j] = 1
                end
            end
        end
        Te = [reverse(Te, dims=2) Te;]
        Te = [reverse(Te, dims=1); Te]
        x[n6 .+ [1:2n6;], n3 - 1 .+ [1:2n3;]] = x[n6 .+ [1:2n6;], n3 - 1 .+ [1:2n3;]] .+ 2Te
        x[findall((i) -> i==3, x)] .= 2

        # Add a triangle
        Te = triu(ones(T, n3, n3))
        mT, nT = size(Te)
        x[n3 + n12 .+ [1:nT;], 1 .+ [1:mT;]] = 3Te

        # add a cross
        Te = zeros(T, 2 * n6 + 1, 2 * n6 + 1)
        mT, nT = size(Te)
        Te[n6+1, 1:nT] = ones(T, nT)
        Te[1:mT, n6+1] = ones(T, mT)
        x[n2 + n12 .+ [1:mT;], n2 .+ [1:nT;]] = 4Te

        x = reshape(x[1:n, 1:n], n^2)
        b = A * x
        return RegProb(A, b, x)
    end
end
blur(::Type{T}, n::Integer, matrixonly::Bool = true) where T = blur(T, n, 3, 0.7, matrixonly)
blur(args...) = blur(Float64, args...)
blur(::Type, args...) = throw(MethodError(blur, Tuple(args)))

"""
Inverse Laplace Transformation
"""
function ilaplace(::Type{T}, n::Integer) where T

end

"""
Stellar Parallax Problem with 26 Fixed, Real Observations
=========================================================
The underlying problem is a Fredholm integral equation of the first
kind with kernel

   `K(s,t) = (1/(σ√(2π)))*exp(-0.5*((s-t)/σ)^2)`,

with `σ = 0.014234`, and it is dscretized by means of a Galerkin method
with `n` orthonormal basis functions. The right-hand side `b` consists of
a measured distribution function of stellar parallaxes, and its length
is fixed at `26`, i.e, the matrix `A` is `26×n`. The exact solution,
which represents the true distribution of stellar parallaxes, is unknown.

*Input options:*

+ [type,] dim, [matrixonly]: the generated matrix is `26×dim`. If
    `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be generated
    (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**W. M. Smart**, Stellar Dynamics, Cambridge University Press, Cambridge,
(1938), p. 30.
"""
function parallax(::Type{T}, n::Integer, matrixonly::Bool=true) where T
    # Intialization
    a = zero(T)
    b = T(0.1)
    m = 26
    σ = T(0.014234)
    hs = T(0.130/m)
    hx = (b-a)/n
    hsh = hs/2
    hxh = hx/2
    ss = (-T(0.03) .+ T[0:m-1;] * hs) * ones(T, n)'
    xx = ones(T, m) * (a .+ T[0:n-1;]' * hx)

    # compute matrix A
    A = 16exp.(-T(0.5).*((ss .+ hsh .- xx .- hxh)./σ).^2)

    A .+= 4(exp.(-T(0.5).*((ss .+ hsh .- xx)./σ).^2) .+
              exp.(-T(0.5).*((ss .+ hsh .- xx .- hx)./σ).^2) .+
              exp.(-T(0.5).*((ss .- xx .- hxh)./σ).^2) .+
              exp.(-T(0.5).*((ss .+ hs .- xx .- hxh)./σ).^2))

    A .+= (exp.(-T(0.5).*((ss .- xx)./σ).^2) .+
           exp.(-T(0.5).*((ss .+ hs .- xx)./σ).^2) .+
           exp.(-T(0.5).*((ss .- xx .- hx)./σ).^2) .+
           exp.(-T(0.5).*((ss .+ hs .- xx .- hx)./σ).^2))

    A = T(sqrt(hs*hx)/(36*σ*sqrt(2*π)))*A

    if matrixonly
        return A
    else
        # compute b
        b = T[3;7;7;17;27;39;46;51;56;50;43;45;43;32;33;29;
             21;12;17;13;15;12;6;6;5;5]/T(sqrt(hs)*640)
        return RegProbNoSolution(A,b)
    end
end
parallax(args...) = parallax(Float64, args...)
parallax(::Type, args...) = throw(MethodError(parallax, Tuple(args)))

"""
Test Problem with \"Spike\" Solution
====================================
Artifically generated discrete ill-posed problem.

*Input options:*

+ [type,] dim, t_max, [matrixonly]: `dim` is the dimension of the
              matrix. `t_max` controls the length of the pulse train.
              If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will be
              generated (`matrixonly = true` by default). The solution x
              consists a unit step at t = .5 and a pulse train of spike
              of decreasing magnitude at t = .5, 1.5, 2.5, ....

+ [type,] dim, [matrixonly]: `t_max = 5`.

*Groups:* ["regprob"]
"""
function spikes(::Type{T}, n::Integer, t_max::Integer, matrixonly::Bool = true) where T
    del = convert(T, t_max/n)

    # compute A
    t, sigma = meshgrid(T[del:del:t_max;], T[del:del:t_max;])
    A = sigma ./ (2 .* sqrt.(π .* t.^3)) .* exp.(-(sigma.^2) ./ (4 .* t))

    if matrixonly
        return A
    else
        # compute b and x
        heights = 2*ones(T, t_max); heights[1] = 25
        heights[2] = 9; heights[3] = 5; heights[4] = 4; heights[5] = 3
        x = zeros(T, n); n_h = one(Integer)
        peak = convert(T, 0.5/t_max); peak_dist = one(T)/ t_max
        if peak < 1
            n_peak = round_matlab(Integer, peak*n)
            x[n_peak] = heights[n_h]
            x[n_peak+1:n] = ones(T, n-n_peak)
            peak = peak + peak_dist; n_h = n_h + one(Integer)
        end
        while peak < 1
            n_peak = round_matlab(Integer, peak*n)
            x[n_peak] = heights[n_h]
            peak = peak + peak_dist; n_h = n_h + 1
        end
        b = A*x
        return RegProb(A, b, x)
    end
end
spikes(::Type{T}, n::Integer, matrixonly::Bool = true) where T =
         spikes(T, n, 5, matrixonly)
spikes(args...) = spikes(Float64, args...)
spikes(::Type, args...) = throw(MethodError(spikes, Tuple(args)))

#
# Two-dimensional tomography problem with sparse matrix
#
function tomo(::Type{T}, n::Integer) where T

end

"""
Integral Equation with No square Integrable Solution
====================================================
Discretization of a first kind Fredholm integral equation with
kernel `K` and right-hand side `g` given by
                    `K(s,t) = 1/(s+t+1), g(s) = 1`,
where both integration intervals are `[0, 1]`. The matrix `A`
is a Hankel matrix.

*Input options:*

+ [type,] dim, [matrixonly]: `dim` is the dimension of the matrix.
              If `matrixonly = false`, the matrix A and vectors b and x in the linear system Ax = b will also
              be generated (`matrixonly = true` by default).

*Groups:* ["regprob"]

*References:*

**F. Ursell**, Introduction to the theory of linear
              integral equations., Chapter 1 in L. M. Delves & J. Walsh,
              Numerical Solution of Integral Equations, Clarendon Press,
              1974.
"""
function ursell(::Type{T}, n::Integer, matrixonly::Bool = true) where T
    r = zeros(T, n); c = copy(r)
    for k = 1:n
        d1 = one(T) + (one(T) + k)/n
        d2 = one(T) + k/n
        d3 = one(T) + (k - one(T))/n
        c[k] = n*(d1*log(d1) + d3*log(d3) - 2*d2*log(d2))
        e1 = one(T) + (n+k)/n
        e2 = one(T) + (n+k-one(T))/n
        e3 = one(T) + (n+k-2)/n
        r[k] = n*(e1*log(e1) + e3*log(e3) - 2*e2*log(e2))
    end
    A = hankel(T, c, r)
    if matrixonly
        return A
    else
        # compute the right-hand side b
        b = ones(T,n)/convert(T, sqrt(n))
        return RegProbNoSolution(A, b)
    end
end
ursell(args...) = ursell(Float64, args...)
ursell(::Type, args...) = throw(MethodError(ursell, Tuple(args)))


text:
function check_kwargs(input_kwargs, valid_set::Array{Symbol}, function_name::String)
    if isempty(input_kwargs)
        return
    else
        for (key, value) in input_kwargs
            if !(key in valid_set)
                throw(ArgumentError("keyword argument $(key) is not a valid input for $(function_name)"))
            end
        end
    end
    return
end


text:
#=

This file specifies how to pretty-print Scattering module types

=#

function Base.show(io::IO, ::MIME"text/plain", x::Aerosol)

    println(io, "\n------------------")
    println(io, "Aerosol Struct")
    println(io, "------------------")
    println(io, "Size distribution: $(x.size_distribution)")
    println(io, "Real refractive Index nᵣ: $(x.nᵣ)")
    println(io, "Imaginary refractive Index nᵣ: $(x.nᵢ)")
    # start,stop = quantile(x.size_distribution,[0.0025,0.9975])
    # x_axis = collect(range(start,stop,20))
    # y = pdf.(x.size_distribution, x_axis)
    # print(io, UnicodePlots.lineplot(x_axis, y,title="Aerosol Size Distribution", xlabel="Radius (μm)" , ylabel="Frequency"))
end

# Overload the show method for ScatteringMatrix
function Base.show(io::IO, x::ScatteringMatrix)

    println(io, "  f₁₁| $(length(x.f₁₁))-element vector")
    println(io, "  f₁₂| $(length(x.f₁₂))-element vector")
    println(io, "  f₂₂| $(length(x.f₂₂))-element vector")
    println(io, "  f₃₃| $(length(x.f₃₃))-element vector")
    println(io, "  f₃₄| $(length(x.f₃₄))-element vector")
    print(io, "  f₄₄| $(length(x.f₄₄))-element vector")
end


text:
using Combinatorics, StatsBase

namesofline(n) = counts([x[1] for x in integer_partitions(n)])

function centerjustpyramid(n)
    maxwidth = length(string(namesofline(n)))
    for i in 1:n
        s = string(namesofline(i))
        println(" " ^ div(maxwidth - length(s), 2), s)
    end
end

centerjustpyramid(25)

const cachecountpartitions = Dict{BigInt,BigInt}()
function countpartitions(n::BigInt)
    if n < 0
        0
    elseif n < 2
        1
    elseif (np = get(cachecountpartitions, n, 0)) > 0
        np
    else
        np = 0
        sgn = 1
        for k = 1:n
            np += sgn * (countpartitions(n - (k*(3k-1)) >> 1) + countpartitions(n - (k*(3k+1)) >> 1))
            sgn = -sgn
        end
        cachecountpartitions[n] = np
    end
end

G(n) = countpartitions(BigInt(n))

for g in [23, 123, 1234, 12345]
    @time println("\nG($g) is $(G(g))")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@inline function _sampleSignedBernoulli(rng)
  return 2*(rand(rng) < 0.5) - 1
end

## A modified version of the algorithm (attributed to S. Vempala) in
## Wästlund, J., 1999. Functions arising by coin flipping. Technical Report,
## KTH, Stockholm.
##
## Specifically, I modified the algorithm so that the calls to f() and the
## sampling of the Bernoulli r.v.s occurs in tandem, so as to avoid the
## infinite expected running time of the initial pass of the original algorithm.
##
## This is also more efficient in terms of calls to f() than the version
## described in
## Flajolet, P., Pelletier, M. and Soria, M., 2011, January. On Buffon machines
## and numbers. In Proc. 22nd ACM-SIAM symposium on Discrete algorithms.
function _wastlundSqrt(f::F, rng::RNG) where
  {F<:Function, RNG<:AbstractRNG}
  flips::Int64 = 0
  x::Bool = f(); flips += 1
  x && (return true, flips)
  s::Int64 = _sampleSignedBernoulli(rng)
  while s > 0
    x = f(); flips += 1
    x && (return true, flips)
    s += _sampleSignedBernoulli(rng) + _sampleSignedBernoulli(rng)
  end
  return false, flips
end


text:
@params mutable struct MeanCompliance{T, TC <: Compliance{T}, TM, TS} <: AbstractFunction{T}
    compliance::TC
    method::TM
    F::TS
    grad_temp::AbstractVector{T}
end
function MeanCompliance(problem::MultiLoad, solver::AbstractDisplacementSolver, args...; method = :exact_svd, sample_once = true, nv = nothing, V = nothing, sample_method = :hutch, kwargs...)
    if method == :exact
        method = ExactMean(problem.F)
    elseif method == :exact_svd
        method = ExactSVDMean(problem.F)
    elseif method == :trace || method == :approx
        if sample_method isa Symbol
            sample_method = sample_method == :hadamard ? hadamard! : hutch_rand!
        end
        if V === nothing
            nv = nv === nothing ? 1 : nv
            method = TraceEstimationMean(problem.F, nv, sample_once, sample_method)
        else
            nv = nv === nothing ? size(V, 2) : nv
            method = TraceEstimationMean(problem.F, view(V, :, 1:nv), sample_once, sample_method)
        end
    else
        if sample_method isa Symbol
            sample_method = sample_method == :hadamard ? hadamard! : hutch_rand!
        end
        if V === nothing
            nv = nv === nothing ? 1 : nv
            method = TraceEstimationSVDMean(problem.F, nv, sample_once, sample_method)
        else
            nv = nv === nothing ? size(V, 2) : nv
            method = TraceEstimationSVDMean(problem.F, view(V, :, 1:nv), sample_once, sample_method)
        end
    end
    comp = Compliance(whichdevice(solver), problem.problem, solver, args...; kwargs...)
    return MeanCompliance(comp, method, problem.F, similar(comp.grad))
end

function (ec::MeanCompliance{T})(x, grad = ec.grad) where {T}
    @unpack compliance, F, method = ec
    @unpack cell_comp, solver, tracing, topopt_trace = compliance
    @unpack elementinfo, u, xmin = solver
    @unpack metadata, Kes, black, white, varind = elementinfo
    @unpack cell_dofs = metadata

    @timeit to "Eval obj and grad" begin
        #if compliance.solver.vars ≈ x && getpenalty(ec).p ≈ getprevpenalty(ec).p
        #    grad .= o.grad
        #    return compliance.comp
        #end

        penalty = getpenalty(ec)
        copyto!(solver.vars, x)
        compliance.fevals += 1
        setpenalty!(solver, penalty.p)
        obj = compute_mean_compliance(ec, method, solver.vars, grad)

        if compliance.logarithm
            compliance.comp = log(obj)
            grad ./= obj
        else
            compliance.comp = obj
        end
        if compliance.grad !== grad
            copyto!(compliance.grad, grad)
        end

        if compliance.tracing
            if ec.reuse
                ec.reuse = false
            else
                push!(topopt_trace.c_hist, obj)
                push!(topopt_trace.x_hist, copy(x))
                if length(topopt_trace.x_hist) == 1
                    push!(topopt_trace.add_hist, 0)
                    push!(topopt_trace.rem_hist, 0)
                else
                    push!(topopt_trace.add_hist, sum(topopt_trace.x_hist[end] .> topopt_trace.x_hist[end-1]))
                    push!(topopt_trace.rem_hist, sum(topopt_trace.x_hist[end] .< topopt_trace.x_hist[end-1]))
                end
            end
        end
    end
    return compliance.comp::T
end

function compute_mean_compliance(ec::MeanCompliance, ::ExactMean, x, grad)
    compute_exact_ec(ec, x, grad, ec.F, size(ec.F, 2))
end
function compute_exact_ec(ec, x, grad, F, n)
    @unpack compliance, grad_temp = ec
    @unpack cell_comp, solver = compliance
    @unpack elementinfo, u, xmin = solver
    @unpack metadata, Kes, black, white, varind = elementinfo
    @unpack cell_dofs = metadata
    penalty = getpenalty(compliance)
    T = eltype(grad)
    obj = zero(T)
    grad .= 0
    for i in 1:size(F, 2)
        @views solver.rhs .= F[:,i]
        solver(assemble_f = false, reuse_chol = (i > 1))
        u = solver.lhs
        obj += compute_compliance(cell_comp, grad_temp, cell_dofs, Kes, u,
                                    black, white, varind, x, penalty, xmin)
        grad .+= grad_temp
    end
    obj /= n
    grad ./= n
    return obj
end

function compute_mean_compliance(ec::MeanCompliance, ap::TraceEstimationMean, x, grad)
    compute_approx_ec(ec, x, grad, ap.F, ap.V, size(ap.F, 2))
end
function compute_approx_ec(ec, x, grad, F, V, n)
    nv = size(ec.method.V, 2)
    @unpack compliance, grad_temp = ec
    @unpack cell_comp, solver = compliance
    @unpack elementinfo, u, xmin = solver
    @unpack metadata, Kes, black, white, varind = elementinfo
    @unpack cell_dofs = metadata
    penalty = getpenalty(compliance)
    T = eltype(grad)
    obj = zero(T)
    grad .= 0
    ec.method.sample_once || ec.method.sample_method(V)
    for i in 1:nv
        @views mul!(solver.rhs, F, V[:,i])
        solver(assemble_f = false, reuse_chol = (i > 1))
        invKFv = solver.lhs
        obj += compute_compliance(cell_comp, grad_temp, cell_dofs, Kes, invKFv,
                                    black, white, varind, x, penalty, xmin)
        grad .+= grad_temp
    end
    obj /= nv * n
    grad ./= nv * n
    return obj
end

function compute_mean_compliance(ec::MeanCompliance, ex::ExactSVDMean, x, grad)
    compute_exact_ec(ec, x, grad, ex.US, ex.n)
end

function compute_mean_compliance(ec::MeanCompliance, ax::TraceEstimationSVDMean, x, grad)
    compute_approx_ec(ec, x, grad, ax.US, ax.V, ax.n)
end

Utilities.getpenalty(c::MeanCompliance) = c.compliance |> getsolver |> getpenalty
@forward_property MeanCompliance compliance

hutch_rand!(x::Array) = x .= rand.(Ref(-1.0:2.0:1.0))
function hadamard3!(V)
	n, nv = size(V)
	H = ones(Int, 1, 1)
	while size(H, 1) < nv
		H = [H H ; H -H]
	end
	H = H[:, 1:nv]
	while size(H, 1) < n
		n1 = nv ÷ 2
		H1 = H[:, 1:n1]
		H2 = H[:, (n1+1):nv]
		H = [H1 H2; H1 -H2]
	end
	V .= H[1:n, :]
	return V
end
function hadamard2!(V)
	n, nv = size(V)
	H = ones(Int, 1, 1)
	while size(H, 1) < nv
		H = [H H ; H -H]
	end
	H = H[:, 1:nv]
	while size(H, 1) < n
		H = [H; -H]
	end
	V .= H[1:n, :]
	return V
end
function hadamard!(V)
	n, nv = size(V)
	H = ones(Int, 1, 1)
	while size(H, 1) < nv
		H = [H H ; H -H]
	end
	H = H[:, 1:nv]
	H = repeat(H, ceil(Int, n/nv), 1)
	V .= H[1:n, :]
	return V
end

function generate_scenarios(dof::Int, size::Tuple{Int, Int}, f, perturb = ()->(rand()-0.5))
    ndofs, nscenarios = size
    I = Int[]; J = Int[]; V = Float64[];
    V = [f * (1 + perturb()) for s in 1:nscenarios]
    I = [dof for s in 1:nscenarios]
    J = 1:nscenarios
    return sparse(I, J, V, ndofs, nscenarios)
end
Nonconvex.getdim(f::MeanCompliance) = 1


text:
"""
    box
        -fully actuated
"""
struct Box{T} <: Model{T}
	# dimensions
	nq::Int # generalized coordinates
	nu::Int # controls
	nw::Int # parameters
	nc::Int # contact points

    mass::T
    inertia::T
    gravity::T

    friction_body_world::Vector{T}
    contact_corner_offset::Vector{Vector{T}}
    friction_joint::Vector{T}
end

# Methods
function lagrangian(model::Box, q, q̇)
    L = 0.0
    L += 0.5 * model.mass * dot(q̇[1:2], q̇[1:2])
    L -= model.mass * model.gravity * q[2]

    L += 0.5 * model.inertia * q̇[3]^2

    return L
end

function kinematics(model::Box, q)
    pos = q[1:2]
    θ = q[3]
    R = rotation_matrix(θ)
    [
        pos + R * model.contact_corner_offset[1];
        pos + R * model.contact_corner_offset[2];
        pos + R * model.contact_corner_offset[3];
        pos + R * model.contact_corner_offset[4];
    ]
end

function signed_distance(model::Box, q)
    kinematics(model, q)[collect([2, 4, 6, 8])]
end

function input_jacobian(model::Box, q)
	[
	 1.0 0.0 0.0;
	 0.0 1.0 0.0;
     0.0 0.0 1.0;
    ]
end

function contact_jacobian(model::Box, q)
	j = kinematics(model, q)
	Symbolics.jacobian(j, q)
end

# nominal configuration
function nominal_configuration(model::Box)
	[0.0; 0.1; 0.0]
end

# friction coefficients
friction_coefficients(model::Box) = model.friction_body_world

# Dimensions
nq = 3 # configuration dimension
nu = 3 # control dimension
nw = 0 # parameters
nc = 4 # number of contact points

# Parameters
r_dim = 0.5
gravity = 10.0#9.81
mass = 1.0 # mass
inertia = 1.0 / 12.0 * mass * ((2.0 * r_dim)^2 + (2.0 * r_dim)^2)

# Contact corner
cc1 = [r_dim, r_dim]
cc2 = [-r_dim, r_dim]
cc3 = [r_dim, -r_dim]
cc4 = [-r_dim, -r_dim]

contact_corner_offset = [cc1, cc2, cc3, cc4]
friction_body_world = [1.0, 1.0, 1.0, 1.0]  # coefficient of friction

# Model
box = Box(nq, nu, nw, nc,
		  mass, inertia, gravity,
		  friction_body_world, contact_corner_offset, zeros(0))

box_contact_kinematics = [q -> kinematics(box, q)[1:2],
                          q -> kinematics(box, q)[3:4],
                          q -> kinematics(box, q)[5:6],
                          q -> kinematics(box, q)[7:8]]

box_contact_kinematics_jacobians = [q -> contact_jacobian(box, q)[1:2, :],
                                    q -> contact_jacobian(box, q)[3:4, :],
                                    q -> contact_jacobian(box, q)[5:6, :],
                                    q -> contact_jacobian(box, q)[7:8, :]]

name(::Box) = :box


text:
using Documenter, BayesianIntegral

makedocs(
    format = Documenter.HTML(),
    sitename = "BayesianIntegral",
    modules = [BayesianIntegral],
    pages = ["index.md",
            "1_bayesian_integration.md",
            "2_training_hyperparameters.md",
            "99_refs.md"]
)

deploydocs(
    repo   = "github.com/s-baumann/BayesianIntegral.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
"""
Extensions for JSON for additional data types
"""
module JSONify

using Requires
import JSON: _writejson, State

"""
    subtypesearch(T::DataType)

Returns: Array of all non-abstract and non-parametrised subtypes of the type `T`.

This function searches recursively through the type hierarchy, whereas `Base.subtypes` searches just 1 level down.

`subtypesearch(concrete_type)` returns an array containing only `concrete_type`.
"""
function subtypesearch(T::DataType)
    function st(T::DataType)
        if isleaftype(T) && isempty(T.parameters)
            produce(T)
        end
        for S in sort(subtypes(T), by=(T) -> (T.size, string(T.name)))
            st(S)
        end
    end
    collect(Task(() -> st(T)))
end

# JSON-serialisable types
json_types = append!([String, Char, Date, DateTime], subtypesearch(Number))
type_map = Dict{String, DataType}(zip([lowercase(string(t.name)) for t in json_types], json_types))

function _writejson(io::IO, state::State, a::DataType)
    if !in(a, json_types)
        throw(ArgumentError("Cannot serialise type: $a"))
    end
    _writejson(io, state, lowercase(string(a.name)))
end

@require DataFrames begin

    """
        DataFrame(data::Dict{String, Vector{Any}})

    Convert dictionary of Arrays to a `DataFrame`.

    The dict is typically constructed by `JSON.parse(JSON.json(my_dataframe))`.

    Assumes a dictionary in the form:
    ```julia
    Dict(
        "names" => Vector{Union{String, Symbol}}(...),
        "types" => Vector{String}(...),
        "columns" => Vector{Any}(...)
    )
    ```
    where

    * `names`: the names of the columns.
    * `types`: the types of the columns.
    * `columns`: the data in the columns.

    and each array is of the same length.
    """
    function DataFrames.DataFrame(data::Dict{String, Any})
        names, types, columns = try
            data["names"], data["types"], data["columns"]
        catch e
            throw(InputError("Missing value: $(e.key)"))
        end

        df = DataFrames.DataFrame()

        for (name, type_name, column) in zip(names, types, columns)
            n = length(column)
            # construct empty DataArray (initially all values are `NA`)
            array = df[Symbol(name)] = DataArrays.DataArray(get(type_map, type_name, Any), n)
            for i = 1:n
                value = column[i]
                # populate if data provided, otherwise leaves value as `NA`
                if value != nothing
                    array[i] = value
                end
            end
        end
        df
    end

    # Custom JSON renderers

    """
    JSONify a DataFrame.

    Converts the DataFrame to a Dict first.
    """
    function _writejson(io::IO, state::State, a::DataFrames.DataFrame)
        # write data out in format of DataFrame converter
        _writejson(io, state, Dict{Symbol, Vector{Any}}(
            :names => DataFrames.names(a),
            :types => DataFrames.eltypes(a),
            :columns => DataFrames.columns(a)
        ))
    end

end

@require DataArrays begin

    """
    JSONify `NAtype` to `null`.
    """
    function _writejson(io::IO, state::State, a::DataArrays.NAtype)
        Base.print(io, "null")
    end

end

# more (de)serialisers added as needed

end


text:
@json struct Milestone
    id::Int
    iid::Int
    project_id::Int
    title::String
    description::String
    state::String
    created_at::DateTime
    updated_at::DateTime
    due_date::Date
    start_date::Date
    web_url::String
end

@json struct TimeStats
    time_estimate::Int
    total_time_spent::Int
    human_time_estimate::String
    human_total_time_spent::String
end

@json struct Pipeline
    id::Int
    sha::String
    ref::String
    status::String
    web_url::String
end

@json struct DiffRefs
    base_sha::String
    head_sha::String
    start_sha::String
end

@json struct MergeRequest
    id::Int
    iid::Int
    project_id::Int
    title::String
    description::String
    state::String
    created_at::DateTime
    updated_at::DateTime
    target_branch::String
    source_branch::String
    upvotes::Int
    downvotes::Int
    author::User
    assignee::User
    source_project_id::Int
    target_project_id::Int
    labels::Vector{String}
    work_in_progress::Bool
    milestone::Milestone
    merge_when_pipeline_succeeds::Bool
    merge_status::String
    merge_error::String
    sha::String
    merge_commit_sha::String
    user_notes_count::Int
    discussion_locked::Bool
    should_remove_source_branch::Bool
    force_remove_source_branch::Bool
    allow_collaboration::Bool
    allow_maintainer_to_push::Bool
    web_url::String
    time_stats::TimeStats
    squash::Bool
    subscribed::Bool
    changes_count::String
    merged_by::User
    merged_at::DateTime
    closed_by::User
    closed_at::DateTime
    latest_build_started_at::DateTime
    latest_build_finished_at::DateTime
    first_deployed_to_production_at::DateTime
    pipeline::Pipeline
    diff_refs::DiffRefs
end

endpoint(::GitLabAPI, ::typeof(get_pull_requests), project::Integer) =
    Endpoint(:GET, "/projects/$project/merge_requests")
into(::GitLabAPI, ::typeof(get_pull_requests)) = Vector{MergeRequest}

endpoint(::GitLabAPI, ::typeof(get_pull_request), project::Integer, number::Integer) =
    Endpoint(:GET, "/projects/$project/merge_requests/$number")
endpoint(
    ::GitLabAPI, ::typeof(get_pull_request),
    owner::AStr, repo::AStr, number::Integer,
) = Endpoint(:GET, "/projects/$(encode(owner, repo))/merge_requests/$number")
into(::GitLabAPI, ::typeof(get_pull_request)) = MergeRequest

endpoint(::GitLabAPI, ::typeof(create_pull_request), project::Integer) =
    Endpoint(:POST, "/projects/$project/merge_requests")
endpoint(::GitLabAPI, ::typeof(create_pull_request), owner::AStr, repo::AStr) =
    Endpoint(:POST, "/projects/$(encode(owner, repo))/merge_requests")
into(::GitLabAPI, ::typeof(create_pull_request)) = MergeRequest

endpoint(::GitLabAPI, ::typeof(update_pull_request), project::Integer, number::Integer) =
    Endpoint(:PUT, "/projects/$project/merge_requests/$number")
endpoint(
    ::GitLabAPI, ::typeof(update_pull_request),
    owner::AStr, repo::AStr, number::Integer,
) = Endpoint(:PUT, "/projects/$(encode(owner, repo))/merge_requests/$number")
into(::GitLabAPI, ::typeof(update_pull_request)) = MergeRequest


text:
mutable struct EBIResponse
   fastq_1_url::String
   fastq_2_url::String
   paired::Bool
   success::Bool
end

function ident_to_fastq_url( ebi_id::String )
   request_url = "http://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=$ebi_id&result=read_run&fields=fastq_ftp&display=txt"
   res = Requests.get(request_url)
   res.status == 200 || error("Couldn't retrieve information for $ebi_id from ebi.ac.uk!")
   buf = IOBuffer(res.data)
   header = chomp(readline(buf))
   if eof(buf)
      error("No fastq data is associated with $ebi_id according to ebi.ac.uk!")
   end
   urls   = split(chomp(readline(buf)), ';', keep=false)
   if length(urls) >= 2
      success = res.status == 200 && urls[1] != "" && urls[2] != "" ? true : false
      ret = EBIResponse(string(urls[1]), string(urls[2]), true, success)
   else
      success = res.status == 200 && urls[1] != "" ? true : false
      ret = EBIResponse(string(urls[1]), "", false, success)
   end
   ret
end


text:
module STMOZOO

# execute your source file and export the module you made
include("sudoku.jl")
export SudokuSolver

end # module


text:
using Documenter, NPDemand

makedocs(sitename = "NPDemand Documentation")
ENV["TRAVIS_PULL_REQUEST"] = false;
ENV["TRAVIS_BRANCH"] = "master";
myTravis = Documenter.Travis();
deploydocs(
    repo = "github.com/jamesbrandecon/NPDemand.jl.git",
    deploy_config = myTravis,
    push_preview = false
)


text:
"""
    ERA5Hourly <: ERA5Module

Specifies that the dataset to be analyzed contains hourly data.  All fields are the same as that specified in the `ERA5Module` docstring.
"""
struct ERA5Hourly{ST<:AbstractString, DT<:TimeType} <: ERA5Module
    modID :: ST
    dtbeg :: Date
    dtend :: Date
    eroot :: ST
end

"""
    ERA5Monthly <: ERA5Module

Specifies that the dataset to be analyzed contains monthly-mean data.  All fields are the same as that specified in the `ERA5Module` docstring.
"""
struct ERA5Monthly{ST<:AbstractString, DT<:TimeType} <: ERA5Module
    modID :: ST
    dtbeg :: Date
    dtend :: Date
    eroot :: ST
end

"""
    ERA5MonthlyHour <: ERA5Module

Specifies that the dataset to be analyzed contains monthly-mean hourly data.  There is one additional field compared to ERA5Module:

Additional fields
=================
- `hours` : A vector of integers containing the hours-of-day to be downloaded/analyzed.
"""
struct ERA5MonthlyHour{ST<:AbstractString, DT<:TimeType} <: ERA5Module
    modID :: ST
    dtbeg :: Date
    dtend :: Date
    hours :: Vector{Int}
    eroot :: ST
end

"""
    ERA5Hourly(;
        dtbeg :: TimeType,
        dtend :: TimeType,
        eroot :: AbstractString = homedir(),
    ) -> ERA5Hourly <: ERA5Module

A function that creates an `ERA5Hourly` module.  All possible hours are downloaded, and data is saved month-by-month.

Keyword Arguments
=================
- `eroot` : The specified directory in which to save the data
- `dtbeg` : The date for which downloads/analysis begins, automatically rounded to the nearest month
- `dtend` : The date for which downloads/analysis finishes, automatically rounded to the nearest month
"""
function ERA5Hourly(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    eroot :: AbstractString = homedir(),
)

    @info "$(modulelog()) - Setting up data structure containing information on the ERA5 Hourly data to be downloaded"
    dtbeg = Date(year(dtbeg),month(dtbeg),1)
	dtend = Date(year(dtend),month(dtend),daysinmonth(dtend))
    checkdates(dtbeg,dtend,)

    return ERA5Hourly{ST,DT}("era5hr",dtbeg,dtend,eroot)

end

"""
    ERA5Monthly(;
        dtbeg :: TimeType,
        dtend :: TimeType,
        eroot :: AbstractString = homedir(),
        hours :: Union{Int,Vector{Int}} = -1,
    ) -> ERA5Monthly <: ERA5Module or ERA5MonthlyHour <: ERA5Module

A function that creates an `ERA5Monthly` or `ERA5MonthlyHour` module depending on the input arguments of `hours`.  Data is saved year-by-year.

Keyword Arguments
=================
- `eroot` : The specified directory in which to save the data
- `dtbeg` : The date for which downloads/analysis begins, automatically rounded to the nearest year
- `dtend` : The date for which downloads/analysis finishes, automatically rounded to the nearest year
- `hours` : If `hours` < 0 or > 24, then it is assumed that we simply want monthly-averages without resolution on hour-of-day and we return an `ERA5Monthly` type.  Otherwise, if an `Int` specifying a single hour-of-day (0-24), or a vector of `Int`s specifying multiple hours of days, then returns an `ERA5MonthlyHour` type.
"""
function ERA5Monthly(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    eroot :: AbstractString = homedir(),
    hours :: Union{Int,Vector{Int}} = -1,
)

    @info "$(modulelog()) - Setting up data structure containing information on the ERA5 Monthly data to be downloaded"
    dtbeg = Date(year(dtbeg),1,1)
    dtend = Date(year(dtend),12,31)
    checkdates(dtbeg,dtend)

    if typeof(hours) <: Int
        if (hours<0) || (hours>24)
            @info "$(modulelog()) - No valid hour of day was specified, and therefore retrieving monthly-averaged data"
            return ERA5Monthly{ST,DT}("era5mo",dtbeg,dtend,eroot)
        else
            if hours == 24; hours = 0 end
            @info "$(modulelog()) - An hour of day ($(hours):00:00) was specified, and therefore retrieving monthly-averaged hourly data"
            return ERA5MonthlyHour{ST,DT}("era5mh",dtbeg,dtend,[hours],eroot)
        end
    else
        hours = sort(unique(mod.(hours,24)))
        @info "$(modulelog()) - Hours of day $(hours) were specified, and therefore retrieving monthly-averaged hourly data for these hours"
        return ERA5MonthlyHour{ST,DT}("era5mh",dtbeg,dtend,hours,eroot)
    end

end

function checkdates(
    dtbeg :: TimeType,
    dtend :: TimeType
)

    if dtend > (now() - Day(5))
        error("$(modulelog()) - You have specified an end date that is in the future of the latest available date of the ERA5 reanalysis dataset")
    end

    if dtend < dtbeg
        error("$(modulelog()) - You have specified an end date that is before your beginning date")
    end

    if (dtbeg < Date(1979,1,1)) || (dtend < Date(1979,1,1))
        error("$(modulelog()) - You have specified a date (be it `dtbeg` or `dtend`) that is before the earliest available date of ERA5 reanalysis data")
    end

end


text:
# Julia wrapper for header: /usr/include/libavfilter/buffersrc.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_buffersrc_get_nb_failed_requests,
    av_buffersrc_parameters_alloc,
    av_buffersrc_parameters_set,
    av_buffersrc_write_frame,
    av_buffersrc_add_frame,
    av_buffersrc_add_frame_flags,
    av_buffersrc_close

function av_buffersrc_get_nb_failed_requests(buffer_src)
    ccall((:av_buffersrc_get_nb_failed_requests, libavfilter), UInt32, (Ptr{AVFilterContext},), buffer_src)
end

function av_buffersrc_parameters_alloc()
    ccall((:av_buffersrc_parameters_alloc, libavfilter), Ptr{AVBufferSrcParameters}, ())
end

function av_buffersrc_parameters_set(ctx, param)
    ccall((:av_buffersrc_parameters_set, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVBufferSrcParameters}), ctx, param)
end

function av_buffersrc_write_frame(ctx, frame)
    ccall((:av_buffersrc_write_frame, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}), ctx, frame)
end

function av_buffersrc_add_frame(ctx, frame)
    ccall((:av_buffersrc_add_frame, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}), ctx, frame)
end

function av_buffersrc_add_frame_flags(buffer_src, frame, flags::Integer)
    ccall((:av_buffersrc_add_frame_flags, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}, Cint), buffer_src, frame, flags)
end

function av_buffersrc_close(ctx, pts::Int64, flags::Integer)
    ccall((:av_buffersrc_close, libavfilter), Cint, (Ptr{AVFilterContext}, Int64, UInt32), ctx, pts, flags)
end


text:
"""
    Dot(;color="", selected=false, col, row)

Generates a dot object containing the following elements:

- `color`: color of dot
- `selected`: dot is selected if true
- `row`: row index
- `col`: col index
"""
mutable struct Dot
    color::String
    selected::Bool
    row::Int
    col::Int
end

function Dot(;color="", selected=false, col, row)
    Dot(color, selected, row, col)
end

"""
    Game(;n_rows=8, n_cols=8, n_rounds=30, round=n_rounds, visible=true)

Generates a game object containing the following elements:

- `dots`: an array of dots representing the board
- `selected_dots`: an array of selected dots
- `n_rounds`: the number of rounds in a game
- `round`: current round
- `score`: the number of connected dots
- `visible`: display GUI if true
"""
mutable struct Game
    dots::Array{Dot,2}
    selected_dots::Vector{Dot}
    n_rounds::Int
    round::Int
    score::Int
    visible::Bool
end

function Game(;n_rows=8, n_cols=8, n_rounds=30, round=n_rounds, visible=true)
    dots = populate(n_rows, n_cols)
    rand_color!.(dots)
    Game(dots, Dot[], n_rounds, round, 0, visible)
end

Broadcast.broadcastable(x::Game) = Ref(x)

mutable struct GUI{G,S}
    gui::G
    style::S
end

function GUI(;width=700, height=600, game = Game())
    gui = GtkWindow("Two Dots", width, height)
    filename = joinpath(@__DIR__, "style.css")
    style = CssProviderLeaf(;filename)
    generate_gui!(game, gui, style)
    GUI(gui, style)
end


text:
# Logistic regression

using SGDOptim

function error_rate(θ::Vector{Float64}, X::Matrix{Float64}, y::Vector{Float64})
    u = X'θ
    countnz(sign(u) .!= sign(y)) / length(y)
end

function logireg_sgd(θ_g::Vector{Float64}, n::Int, σ::Float64)

    # prepare experimental data
    d = length(θ_g)
    X = randn(d, n)
    y = sign(vec(θ_g'X) + σ * randn(n))

    # construct the risk model
    rmodel = riskmodel(LinearPred(d), LogisticLoss())

    # initialize solution
    θ_0 = randn(d)

    # optimize
    θ = sgd(rmodel, θ_0,
        minibatch_seq(X, y, 10),          # configure the way data are supplied
        reg = SqrL2Reg(1.0e-4),           # regularization
        lrate = t->1.0 / (100.0 + t),     # learing rate policy
        cbinterval = 100,                 # how frequently callback is invoked
        callback = simple_trace)          # the callback function

    # compare solution with initial guess
    println()
    @printf("Initial  :  error.rate = %5.2f%%\n", error_rate(θ_0, X, y) * 100.0)
    @printf("Solution :  error.rate = %5.2f%%\n", error_rate(θ,   X, y) * 100.0)
    @printf("gTruth   :  error.rate = %5.2f%%\n", error_rate(θ_g, X, y) * 100.0)
end

logireg_sgd([3.0, 5.0, 2.0], 10000, 0.2)


text:
function calculate_full_rules_vector(data_dictionary::Dict{String,Any}, value_array::BitArray{2})

	# Initialize -
	v = zeros(159)

	function estimate_default(data_dictionary)
		return 1.0
	end
	default() = true

	function calculate_fn_value(x, value_array)
		return value_array[x]
	end
	fn(x) = calculate_fn_value(x, value_array)

	v[1] = default()
	v[2] = ~(fn(158))
	v[3] = ~(fn(159))
	v[4] = max(~(fn(3)), fn(110))
	v[5] = max(~(fn(3)), fn(110))
	v[6] = default()
	v[7] = default()
	v[8] = default()
	v[9] = default()
	v[10] = fn(11)
	v[11] = ~(fn(169))
	v[12] = default()
	v[13] = default()
	v[14] = default()
	v[15] = default()
	v[16] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[17] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[18] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[19] = max(max(~(max(fn(152), fn(43))), fn(111)), fn(110))
	v[20] = default()
	v[21] = default()
	v[22] = default()
	v[23] = default()
	v[24] = max(~(fn(43)), fn(152))
	v[25] = max(~(fn(43)), fn(152))
	v[26] = default()
	v[27] = default()
	v[28] = default()
	v[29] = default()
	v[30] = default()
	v[31] = default()
	v[32] = max(fn(152), min(fn(43), fn(111)))
	v[33] = max(fn(152), min(fn(43), fn(111)))
	v[34] = max(fn(152), min(fn(43), fn(111)))
	v[35] = default()
	v[36] = default()
	v[37] = max(~(fn(48)), ~(fn(2)))
	v[38] = default()
	v[39] = max(fn(163), ~(fn(160)))
	v[40] = max(max(~(fn(168)), ~(min(fn(168), fn(2)))), fn(110))
	v[41] = default()
	v[42] = default()
	v[43] = ~(fn(168))
	v[44] = default()
	v[45] = default()
	v[46] = default()
	v[47] = max(~(fn(163)), min(fn(167), ~(fn(111))))
	v[48] = ~(fn(163))
	v[49] = default()
	v[50] = default()
	v[51] = ~(fn(152))
	v[52] = ~(max(fn(152), fn(43)))
	v[53] = default()
	v[54] = ~(fn(2))
	v[55] = fn(2)
	v[56] = default()
	v[57] = ~(max(fn(67), fn(164)))
	v[58] = default()
	v[59] = default()
	v[60] = default()
	v[61] = max(fn(155), ~(fn(48)))
	v[62] = max(fn(155), ~(fn(48)))
	v[63] = max(fn(155), ~(fn(48)))
	v[64] = default()
	v[65] = default()
	v[66] = default()
	v[67] = fn(157)
	v[68] = default()
	v[69] = default()
	v[70] = default()
	v[71] = ~(max(fn(152), fn(43)))
	v[72] = ~(max(fn(152), fn(43)))
	v[73] = ~(max(fn(152), fn(43)))
	v[74] = ~(max(fn(152), fn(43)))
	v[75] = ~(max(fn(152), fn(43)))
	v[76] = ~(max(fn(152), fn(43)))
	v[77] = ~(max(fn(152), fn(43)))
	v[78] = ~(max(fn(152), fn(43)))
	v[79] = ~(max(fn(152), fn(43)))
	v[80] = ~(max(fn(152), fn(43)))
	v[81] = ~(max(fn(152), fn(43)))
	v[82] = ~(max(fn(152), fn(43)))
	v[83] = ~(max(fn(152), fn(43)))
	v[84] = default()
	v[85] = default()
	v[86] = default()
	v[87] = default()
	v[88] = default()
	v[89] = default()
	v[90] = default()
	v[91] = default()
	v[92] = default()
	v[93] = max(fn(152), min(fn(43), fn(111)))
	v[94] = default()
	v[95] = default()
	v[96] = default()
	v[97] = default()
	v[98] = default()
	v[99] = default()
	v[100] = default()
	v[101] = min(~(fn(152)), fn(103))
	v[102] = min(~(fn(152)), fn(103))
	v[103] = fn(160)
	v[104] = ~(fn(10))
	v[105] = max(fn(111), fn(43))
	v[106] = max(fn(111), fn(43))
	v[107] = ~(min(fn(156), fn(164)))
	v[108] = ~(min(fn(156), fn(164)))
	v[109] = ~(fn(152))
	v[110] = fn(171)
	v[111] = fn(154)
	v[112] = default()
	v[113] = default()
	v[114] = default()
	v[115] = min(min(fn(155), ~(fn(152))), fn(143))
	v[116] = ~(fn(152))
	v[117] = default()
	v[118] = default()
	v[119] = default()
	v[120] = default()
	v[121] = default()
	v[122] = default()
	v[123] = default()
	v[124] = default()
	v[125] = default()
	v[126] = default()
	v[127] = ~(fn(167))
	v[128] = fn(111)
	v[129] = default()
	v[130] = default()
	v[131] = default()
	v[132] = max(fn(152), fn(43))
	v[133] = max(fn(152), fn(43))
	v[134] = default()
	v[135] = default()
	v[136] = min(~(fn(138)), max(~(fn(152)), fn(2)))
	v[137] = min(~(fn(138)), max(~(fn(152)), fn(2)))
	v[138] = fn(39)
	v[139] = default()
	v[140] = default()
	v[141] = default()
	v[142] = max(max(fn(43), fn(111)), fn(143))
	v[143] = fn(144)
	v[144] = max(max(fn(170), fn(162)), fn(166))
	v[145] = max(fn(43), fn(143))
	v[146] = max(fn(43), fn(143))
	v[147] = max(fn(43), fn(143))
	v[148] = max(fn(43), fn(143))
	v[149] = default()
	v[150] = default()
	v[151] = default()
	v[152] = ~(fn(168))
	v[153] = default()
	v[154] = ~(fn(163))
	v[155] = ~(max(max(fn(163), fn(166)), fn(165)))
	v[156] = fn(157)
	v[157] = fn(127)
	v[158] = max(min(~(fn(172)), ~(max(max(fn(182), fn(181)), fn(178)))), fn(161))
	v[159] = min(~(max(fn(176), fn(175))), ~(max(max(max(max(fn(173), fn(179)), fn(177)), fn(174)), fn(180))))

	return v
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8ae2593a-cb48-44c9-b9dc-ab901c358a06
begin
	import Pkg
	Pkg.develop("PoreMatMod")
end

# ╔═╡ 16f0e183-f48d-4dcd-8751-d3c61c875e18
# load required packages (Pluto.jl will automatically install them)
using PoreMatMod, PlutoUI

# ╔═╡ 3cca5289-a829-4717-b796-0834229995d1
using PoreMatMod.ExampleHelpers

# ╔═╡ 8d523993-6e85-443a-9949-12030552b457
md"""
## example: generate a hypothetical MOF structure by decorating its linkers with functional groups
"""

# ╔═╡ 74b45651-21d5-4332-a4a5-866ea1bb02b8
input_file_message()

# ╔═╡ 3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
xtal_folder_message()

# ╔═╡ ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
rc[:paths][:crystals]

# ╔═╡ b66dd0d5-5bac-4b23-a33f-17305b0d4728
moiety_folder_message()

# ╔═╡ a76e79ee-c82c-4771-818d-5380a5fb4c18
rc[:paths][:moieties]

# ╔═╡ be86f07e-0669-46c2-8c79-80e65dfcc6f2
md"""
!!! example \"the task\"
	we have the IRMOF-1 crystal structure, and wish to append acetamido functional groups to six of its (randomly chosen) BDC (1,4-benzodicarboxylate) linkers to give a mixed-linker IRMOF-1 derivative.

**Parent crystal structure**: first, we read in the `.cif` file describing the parent IRMOF-1 crystal structure.
"""

# ╔═╡ b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
begin
	# read in the parent xtal
	parent = Crystal("IRMOF-1.cif")	# load .cif file
	infer_bonds!(parent, true)      # infer bonds
	view_structure(parent)          # view structure
end

# ╔═╡ b402e79a-784b-4d8b-82f1-df4fe1cedad1
fragment_construction_note()

# ╔═╡ 9b2d534a-4f78-4950-add1-9ba645669bb9
md"""
**Query fragment**: next, we construct a query fragment as a _p_-phenylene fragment to match that on the BCD linker of the IRMOF-1 parent structure. we mark one hydrogen atom on the query fragment as "masked" (shown in pink) by tagging its species label with `!` in the input. we need to mask this hydrogen atom because it will eventually be replaced by the acetamido functional group.
"""

# ╔═╡ 4be03110-61ab-4cd6-b3fe-7d51ac5ee771
query = moiety("2-!-p-phenylene.xyz");

# ╔═╡ e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
view_query_or_replacement("2-!-p-phenylene.xyz")

# ╔═╡ 559ef3c3-a176-4d65-8958-810c9b0b32c5
with_terminal() do
	display_query_or_replacement_file("2-!-p-phenylene.xyz")
end

# ╔═╡ d1aa8a19-3702-40f6-b73e-b9ebfc1a7a71
md"""
**Replacement fragment**: next, we construct a replacement fragment as a modified version of the query fragment (with acetamido group in place of one hydrogen atom).
"""

# ╔═╡ 44e7e665-ae68-4cd4-b45f-138a0fb8910e
replacement = moiety("2-acetylamido-p-phenylene.xyz");

# ╔═╡ af606f02-fbd5-4033-a5f3-56a8f740b1a1
view_query_or_replacement("2-acetylamido-p-phenylene.xyz")

# ╔═╡ c8f8dbd3-4191-460d-944f-f5e456ce8b83
with_terminal() do
	display_query_or_replacement_file("2-!-p-phenylene.xyz")
end

# ╔═╡ 127a2bef-3903-4801-bc75-00a6dde2bc6e
md"""
**Find and replace**: finally, we (i) search the parent MOF for the query fragment and (ii) effect the replacements. Ta-da; we have a hypothetical derivatized IRMOF-1 structure. 🎈

n.b.
* the `nb_loc=6` kwarg indicates that we wish to randomly select 6 matches on the IRMOF-1 parent structure to effect the replacement. the `loc` kwarg grants more control over which BDC linkers are functionalized.
* the acetamido functional groups are appended at random positions on each BDC ligand. the `ori` kwarg grants more control over which positions on each linker are functionalized.
* if we omit `nb_loc=6` as a kwarg, a functional group is appended on all BDC linkers of the parent.
"""

# ╔═╡ 74aa19d2-b1a4-4333-9ff9-e6ea74e7d989
begin
	child = replace(parent, query => replacement, nb_loc=6)
	view_structure(child)
end

# ╔═╡ 986ecdc4-455f-457e-a964-f00ddfeb53a2
write_cif_message()

# ╔═╡ 71209370-c445-4ff2-a873-b6e31c46419b
write_cif(child, "acetamido-functionalized_IRMOF-1.cif")

# ╔═╡ Cell order:
# ╟─8ae2593a-cb48-44c9-b9dc-ab901c358a06
# ╟─8d523993-6e85-443a-9949-12030552b457
# ╠═16f0e183-f48d-4dcd-8751-d3c61c875e18
# ╠═3cca5289-a829-4717-b796-0834229995d1
# ╠═74b45651-21d5-4332-a4a5-866ea1bb02b8
# ╟─3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
# ╠═ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
# ╟─b66dd0d5-5bac-4b23-a33f-17305b0d4728
# ╠═a76e79ee-c82c-4771-818d-5380a5fb4c18
# ╟─be86f07e-0669-46c2-8c79-80e65dfcc6f2
# ╠═b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
# ╟─b402e79a-784b-4d8b-82f1-df4fe1cedad1
# ╟─9b2d534a-4f78-4950-add1-9ba645669bb9
# ╠═4be03110-61ab-4cd6-b3fe-7d51ac5ee771
# ╟─e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
# ╟─559ef3c3-a176-4d65-8958-810c9b0b32c5
# ╟─d1aa8a19-3702-40f6-b73e-b9ebfc1a7a71
# ╠═44e7e665-ae68-4cd4-b45f-138a0fb8910e
# ╟─af606f02-fbd5-4033-a5f3-56a8f740b1a1
# ╟─c8f8dbd3-4191-460d-944f-f5e456ce8b83
# ╟─127a2bef-3903-4801-bc75-00a6dde2bc6e
# ╠═74aa19d2-b1a4-4333-9ff9-e6ea74e7d989
# ╟─986ecdc4-455f-457e-a964-f00ddfeb53a2
# ╠═71209370-c445-4ff2-a873-b6e31c46419b


text:
using aI_Test
using Test
using WAV

# OLD EXAMPLE TESTS vvv
@testset "aI_Test.jl" begin
    # 2x + 3y
    @test aI_Test.my_f(2,1) == 7
    @test aI_Test.my_f(2,3) == 13
    @test aI_Test.my_f(1,3) == 11
end

@testset "more testing" begin
    @test derivative_of_my_f(2,1) == 2
end
# OLD EXAMPLE TESTS ^^^

# NEW aI TESTS vvv
# @testset "some WAV testing" begin
#     @test y, fs = wavread("example.wav")
#;pwd()
#file_dir = joinpath(dirname(@__FILE__), "test")

@testset "real aI testing" begin
    #hammerWav = joinpath(dirname(@__FILE__),"Hammerhead.wav")
    example02_ = joinpath(dirname(@__FILE__),"example02.wav")
    pgm_ = joinpath(dirname(@__FILE__),"PGM102T1.wav")

    @test isfile(joinpath(dirname(@__FILE__),"example.wav"))
    @test isfile(joinpath(dirname(@__FILE__),"example02.wav"))
    @test isfile(joinpath(dirname(@__FILE__),"Hammerhead.wav"))
    @test isfile(pgm_)
    #y2, fs2 = WAV.wavread(hammerWav)
    y2, fs2 = WAV.wavread(example02_)
    y_pgm, fs_pgm, nbits_pgm, opt_pgm = WAV.wavread(pgm_)

    #fs2 = 48000.0
    @test fs2 == 16000.0 # versus 48000.0
    @test fs_pgm == 44100
    thisIt_pgm = WAV.getformat(opt_pgm)
    @test thisIt_pgm.nchannels == 2
    #sqrt(sum(y2.^2.) / length(y2))
    #@test sqrt(sum(y2.^2.) / length(y2)) ≈ 0.16209680234518886
    @test sqrt(sum(y2.^2.) / length(y2)) ≈ 0.7070846851353825
    @test sqrt(sum(x -> x*x, y2) / length(y2)) ≈ 0.7070846851353825

    left_pgm = y_pgm[1:end,1]
    right_pgm = y_pgm[1:end,2]
    @test sqrt(sum(x -> x*x, left_pgm) / length(left_pgm)) ≈ 0.06207601152097682
    @test sqrt(sum(x -> x*x, right_pgm) / length(right_pgm)) ≈ 0.06372897603090126

#    isfile("example.wav")
end

#fs2

### COMMENTED vv
# example02_ = joinpath(dirname(@__FILE__),"example02.wav")
# isfile(example02_)
# y2, fs2, nbits, opt = WAV.wavread(example02_)
# fs2
# nbits
# opt
# thisIt = WAV.getformat(opt)
#
# hammerhead_ = joinpath(dirname(@__FILE__),"Hammerhead.wav")
# isfile(hammerhead_)
# y_hammer, fs_hammer, nbits_hammer, opt_hammer = WAV.wavread(hammerhead_)
# fs_hammer
# nbits_hammer
# opt_hammer
# thisIt_hammer = WAV.getformat(opt_hammer)
#
#
# pgm_ = joinpath(dirname(@__FILE__),"PGM102T1.wav")
# isfile(pgm_)
# y_pgm, fs_pgm, nbits_pgm, opt_pgm = WAV.wavread(pgm_)
# fs_pgm
# nbits_pgm
# opt_pgm
# left_pgm = y_pgm[1:end,1]
# right_pgm = y_pgm[1:end,2]
# sqrt(sum(x -> x*x, left_pgm) / length(left_pgm))
# sqrt(sum(x -> x*x, right_pgm) / length(right_pgm))
# thisIt_pgm = WAV.getformat(opt_pgm)
### COMMENTED ^^

# NEW aI TESTS ^^^
#;pwd()
#sqrt(sum(y2.^2.) / length(y2)) == 0.16209680234518886