text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code 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 Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "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 LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES 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.
#
# __END_LICENSE__

include("types.jl")
include("optimize.jl")
include("admm.jl")
include("setup_admm.jl")

"""
    deconvolve(img, psf, γ, σ, reg, method, options)

Deconvolve a focal stack image (img). A 3D psf, camera gain and readnoise variance
have to be known.

The detail of the algorithm is available from
"A convex 3D deconvolution algorithm for low photon count fluorescence imaging"
Scientific Reports 8, Article number: 11489 (2018)
Hayato Ikoma, Michael Broxton, Takamasa Kudo, Gordon Wetzstein
"""
function deconvolve(
    img::Array{T,3}, # 3D captured image [width x height x depth]
    psf::Array{Float32,3}, # 3D psf
    γ::Real, # camera gain
    σ::Real, # camera readnoise std
    reg::Real, # regularization parameter
    method::ADMM = ADMM();
    options::DeconvolutionOptions,
) where {T<:Real}

    @assert all(.!isnan.(img))
    @assert γ > 0
    @assert σ > 0

    I, J, K = size(img)

    # Convert the unit to photoelectron number
    std_img = img ./ γ
    std_σ = σ / γ

    if isnan(method.ρ)
        method.ρ = 6000.0 * reg / maximum(std_img)
    end

    optimizer = setup_optimizer(
        method,
        Float32.(std_img),
        psf,
        Float32(std_σ^2),
        Float32(reg),
        Float32(method.ρ),
    )
    result = optimize(optimizer, options)
    result.x = collect(result.x[1:I, 1:J, 1:K])

    # Release GPU memory
    gc()

    return result
end


text:
# This file was generated, do not modify it. # hide
using Random # hide
rng = Random.GLOBAL_RNG
rand(rng, Float32, Normal())

ext: jl

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
using HDF5

function save_params(file::HDF5File, rbm::RBM, name::AbstractString)
    write(file, "$(name)___weight", rbm.W')
    write(file, "$(name)___vbias", rbm.vbias)
    write(file, "$(name)___bias", rbm.hbias)
end

function load_params(file::HDF5File, rbm::RBM, name::AbstractString)
    rbm.W = read(file, "$(name)___weight")'
    rbm.vbias = read(file, "$(name)___vbias")
    rbm.hbias = read(file, "$(name)___bias")
end

function save_params(file::HDF5File, net::Net)
    for i=1:length(net)
        save_params(file, net[i], getname(net, i))
    end
end
save_params(path::AbstractString, net::Net) = h5open(path, "w") do h5
    save_params(h5, net)
end

function load_params(file, net::Net)
    for i=1:length(net)
        load_params(file, net[i], getname(net, i))
    end
end
load_params(path::AbstractString, net::Net) = h5open(path) do h5
    load_params(h5, net)
end


text:
const default_initialnumvertices = 16
const default_minnumvertices = 4

mutable struct Helper
  numvertices::Int
  minnumvertices::Int
  timetaken::Float64
  function Helper(
      initialnumvertices::Int=default_initialnumvertices,
      minnumvertices::Int=default_minnumvertices,
      timetaken=0.0)
    @assert minnumvertices <= initialnumvertices
    return new(initialnumvertices, minnumvertices, timetaken)
  end
end

function indicesonhypercubesurface(lengths::AbstractVector{Int})
  indices = Set{Vector{Int}}()
  for i in CartesianIndices(Tuple(lengths))
    index = collect(Tuple(i))
    (any(index .== 1) || any(index .== lengths)) || continue
    push!(indices, index)
  end
  return indices
end

function sortgriddedpoints(A::Point, B::Point, n::Int)
  @assert n >= 2
  midpoint = middle(A, B)
  indices = indicesonhypercubesurface(n * ones(Int, 2))
  fractions = ((index .- 1) // (n - 1) for index ∈ indices)
  points = [Point(A + (B - A) * fraction) for fraction ∈ fractions]
  sort!(points, lt=(a, b)->isless(angle(a - midpoint), angle(b - midpoint)))
  return points
end

function distributevertices(f::T, A::Point, B::Point, helper::Helper,
                            vertices=Vector{Vertex}()) where {T<:Function}
  number_points_along_edge = 1 + Int(round(helper.numvertices / 4))
  output = Vector{Vertex}()
  for point ∈ sortgriddedpoints(A, B, number_points_along_edge)
    i = findfirst(v->isequal(v, point), vertices)
    isnothing(i) && push!(output, Vertex(point, Point(f(Vector(point)))))
    isnothing(i) || push!(output, vertices[i])
  end
  return output
end

function refine(A::Point, B::Point)
  M = (A + B) / 2
  return [(A, M), (Point(A.x, M.y), Point(M.x, B.y)),
          (Point(M.x, A.y), Point(B.x, M.y)), (M, B)]
end

function windingangle(vertices, i)
  j = periodicindex(vertices, i+1)
  return angle(vertices[i].value, vertices[j].value)
end
function windingangle(vertices::AbstractVector{Vertex})
  return mapreduce(i -> windingangle(vertices, i), +, eachindex(vertices))
end

"""
The winding number of a complex function around a closed loop
indicates the number of roots of the function encircled by the loop.

### Arguments
 * f: function that maps points in 2D space, ``R^2``,  to to a complex number
 * A: vector of real numbers that represents the lower left corner in ``R^2``
 * B: vector of real numbers that represents the upper right corner in ``R^2``

### Returns
The winding number of the function f: the number roots inside the square A->B
"""
function windingnumber(f::U, A::AbstractVector{T}, B::AbstractVector{T},
                       helper::Helper=Helper()) where {T<:Real, U<:Function}
  return windingnumber(distributevertices(f, Point(A), Point(B), helper))
end
function windingnumber(vertices::AbstractVector{Vertex})
  return windingangle(vertices)/(2π)
end

function boundingboxes!(solutions, f::T, A::Point, B::Point,
    vertices::AbstractVector{Vertex}, helper::Helper,
    xtol_rel, stopval, timelimit) where {T<:Function}
  t = @elapsed nfloat = round(windingnumber(vertices))
  n = isnan(nfloat) ? 0 : Int(nfloat)
  helper.timetaken += t
  ((helper.timetaken > timelimit) || iszero(n)) && return nothing
  isconverged = all(Vector(B - A) .< xtol_rel) ||
    (sum(abs.(f(Vector((A + B) / 2)))) <= stopval)
  if isconverged
    push!(solutions, (Vector(A), Vector(B)))
  else
    helper.numvertices = max(
      Int(round(helper.numvertices / 2)),
      Int(round(default_minnumvertices)))
    for q ∈ refine(A, B)
      newvertices = distributevertices(f, q..., helper, vertices)
      boundingboxes!(solutions, f, q..., newvertices, helper,
                     xtol_rel, stopval, timelimit)
    end
  end
  return nothing
end

function boundingboxes(f::T, A::AbstractVector{U}, B::AbstractVector{U};
                       kwargs...) where {T<:Function, U<:Real}
  @assert length(A) == length(B) == 2
  solutions = Vector{Tuple{Vector{Rational},Vector{Rational}}}()
  kwargs = Dict(kwargs)
  xtol_rel = get(kwargs, :xtol_rel, sqrt(eps()))
  xtol_abs = get(kwargs, :xtol_abs, 0.0)
  stopval = get(kwargs, :stopval, sqrt(eps()))
  timelimit = get(kwargs, :timelimit, 600)
  initialnumvertices = get(kwargs, :initialnumvertices, default_initialnumvertices)
  minnumvertices = get(kwargs, :minnumvertices, min(initialnumvertices,
                                                    default_minnumvertices))
  helper = get(kwargs, :helper, Helper(initialnumvertices, minnumvertices))
  scale(fraction) = A .+ (B .- A) .* fraction
  f_local(fraction) = f(scale(fraction))
  lower = Point(zeros(Rational, 2))
  upper = Point(ones(Rational, 2))
  vertices = distributevertices(f_local, lower, upper, helper)
  boundingboxes!(solutions, f_local, lower, upper, vertices, helper,
                        xtol_rel, stopval, timelimit)
  return [(scale(solution[1]), scale(solution[2])) for solution ∈ solutions]
end

max_stars_repo_name: jwscook/WindingNumbers

max_issues_repo_name: jwscook/WindingNumbers

max_forks_repo_name: jwscook/WindingNumbers


text:
# quadratic function

export Quadratic

### ABSTRACT TYPE

abstract type Quadratic <: ProximableFunction end

is_convex(f::Quadratic) = true
is_smooth(f::Quadratic) = true
is_quadratic(f::Quadratic) = true

fun_name(f::Quadratic) = "Quadratic function"

### CONSTRUCTORS

"""
**Quadratic function**

    Quadratic(Q, q; iterative=false)

For a matrix `Q` (dense or sparse, symmetric and positive semidefinite) and a vector `q`, returns the function
```math
f(x) = \\tfrac{1}{2}\\langle Qx, x\\rangle + \\langle q, x \\rangle.
```
By default, a direct method (based on Cholesky factorization) is used to evaluate `prox!`.
If `iterative=true`, then `prox!` is evaluated approximately using an iterative method instead.
"""
function Quadratic(Q::M, q::V; iterative=false) where {M, V}
    if iterative == false
        QuadraticDirect(Q, q)
    else
        QuadraticIterative(Q, q)
    end
end

### INCLUDE CONCRETE TYPES

include("quadraticDirect.jl")
include("quadraticIterative.jl")


text:
import Random.rand!
using TensorCast

export EuclideanRandomFeatures

"""
    EuclideanRandomFeatures

A set of Euclidean random features, parameterized by frequency and phase,
together with a set of associated basis weights.
"""
mutable struct EuclideanRandomFeatures{A<:AbstractArray{<:Any,3},M<:AbstractMatrix} <: PriorBasis
  frequency  :: A
  phase      :: M
  weights    :: M
end

Flux.trainable(f::EuclideanRandomFeatures) = ()
Flux.@functor EuclideanRandomFeatures

"""
    rand!(self::EuclideanRandomFeatures, k::EuclideanKernel,
          num_features::Int = size(self.frequency,ndims(self.frequency)))

Draw a new set of random features, by randomly sampling a new frequencies
from the spectral measure, and new phases uniformly from ``(0, 2\\pi)``.
Does NOT automatically resample the GP containing the features.
"""
function rand!(self::EuclideanRandomFeatures, k::EuclideanKernel; num_features::Int = size(self.frequency,ndims(self.frequency)))
  (id,od) = k.dims
  (_,s) = size(self.weights)
  Fl = eltype(self.frequency)
  self.frequency = spectral_distribution(k; num_samples = num_features)
  self.phase = Fl(2*pi) .* rand!(similar(self.phase, (od,num_features)))
  self.weights = randn!(similar(self.weights, (num_features, s)))
  nothing
end

"""
    EuclideanRandomFeatures(k::EuclideanKernel, num_features::Int)

Create a set of Euclidean random features with eigenvalues given by the
spectral distribution given by ``k``.
"""
function EuclideanRandomFeatures(k::EuclideanKernel, num_features::Int)
  (id,od) = k.dims
  frequency = zeros(id,od,num_features)
  phase = zeros(od,num_features)
  weights = zeros(num_features,1)
  features = EuclideanRandomFeatures(frequency, phase, weights)
  rand!(features, k)
  features
end

"""
    (self::EuclideanRandomFeatures)(x::AbstractMatrix, w::AbstractMatrix,
                                    k::EuclideanKernel)

Evaluate the ``f(x)`` where ``f`` is a Gaussian process with kernel ``k``,
and ``x`` is the data, using the random features.
"""
function (self::EuclideanRandomFeatures)(x::AbstractMatrix, k::EuclideanKernel)
  Fl = eltype(self.frequency)
  l = size(self.frequency, ndims(self.frequency))
  (outer_weights, inner_weights) = spectral_weights(k, self.frequency)
  @cast rescaled_x[ID,N] := x[ID,N] / inner_weights[ID]
  @matmul basis_fn_inner_prod[OD,L,N] := sum(ID) self.frequency[ID,OD,L] * rescaled_x[ID,N]
  @cast basis_fn[OD,L,N] := cos(basis_fn_inner_prod[OD,L,N] + self.phase[OD,L])
  basis_weight = outer_weights .* sqrt(Fl(2)) ./ sqrt(Fl(l)) .* self.weights
  @matmul output[OD,N,S] := sum(L) basis_fn[OD,L,N] * basis_weight[L,S]
  output
end

"""
    (self::EuclideanRandomFeatures)(x::AbstractMatrix, w::AbstractMatrix,
                                    k::GradientKernel{<:EuclideanKernel})

Evaluate ``(\\nabla g)(x)`` where ``g`` is a Gaussian process with kernel ``k``,
``\\nabla`` is the gradient inter-domain operator, and ``x`` is the data, using
the random features.
"""
function (self::EuclideanRandomFeatures)(x::AbstractMatrix, k::GradientKernel{<:EuclideanKernel})
  Fl = eltype(self.frequency)
  l = size(self.frequency, ndims(self.frequency))
  (outer_weights, inner_weights) = spectral_weights(k, self.frequency)
  @cast rescaled_x[ID,N] := x[ID,N] / inner_weights[ID]
  @matmul basis_fn_inner_prod[OD,L,N] := sum(ID) self.frequency[ID,OD,L] * rescaled_x[ID,N]
  @cast basis_fn_grad_outer[OD,L,N] := -sin(basis_fn_inner_prod[OD,L,N] + self.phase[OD,L])
  @cast basis_fn_grad[ID,OD,L,N] := basis_fn_grad_outer[OD,L,N] * self.frequency[ID,OD,L] / inner_weights[ID]
  basis_weight = outer_weights .* sqrt(Fl(2)) ./ sqrt(Fl(l)) .* self.weights
  @matmul output[ID,OD,N,S] := sum(L) basis_fn_grad[ID,OD,L,N] * basis_weight[L,S]
  dropdims(output; dims=2)
end

"""
    (self::EuclideanRandomFeatures)(x::AbstractArray{<:Any,3},w::AbstractMatrix,
                                    k::GradientKernel{<:EuclideanKernel})

Evaluate ``(\\nabla g)(x)`` where ``g`` is a Gaussian process with kernel ``k``,
``\\nabla`` is the gradient inter-domain operator, and ``x`` is the batched
data, using the random features.
"""
function (self::EuclideanRandomFeatures)(a::AbstractArray{<:Any,3}, k::GradientKernel{<:EuclideanKernel})
  (d,n,s) = size(a)
  x = reshape(a, (:,n*s))
  Fl = eltype(self.frequency)
  l = size(self.frequency, ndims(self.frequency))
  (outer_weights, inner_weights) = spectral_weights(k, self.frequency)
  @cast rescaled_x[ID,N] := x[ID,N] / inner_weights[ID]
  @matmul basis_fn_inner_prod[OD,L,N] := sum(ID) self.frequency[ID,OD,L] * rescaled_x[ID,N]
  @cast basis_fn_grad_outer[OD,L,N] := -sin(basis_fn_inner_prod[OD,L,N] + self.phase[OD,L])
  @cast basis_fn_grad[ID,OD,L,N] := basis_fn_grad_outer[OD,L,N] * self.frequency[ID,OD,L] / inner_weights[ID]
  basis_weight = outer_weights .* sqrt(Fl(2)) ./ sqrt(Fl(l)) .* self.weights
  basis_fn_grad_batched = reshape(basis_fn_grad, (d,l,n,s))
  @reduce output[ID,N,S] := sum(L) basis_fn_grad_batched[ID,L,N,S] * basis_weight[L,S]
  output
end


text:
@testset "Matrix: $TRef" for (TMat, T, TRef) in [
    (ArbMatrix, Arb, Arb),
    (AcbMatrix, Acb, Acb),
    (ArbRefMatrix, Arb, ArbRef),
    (AcbRefMatrix, Acb, AcbRef),
]
    @testset "Basic" begin
        M = TMat(4, 4, prec = 128)
        @test size(M) == (4, 4)
        @test precision(M) == 128

        x = T(1.5)
        M[2, 2] = x
        @test M[2, 2] == x
        @test M[2, 2] isa TRef
        @test precision(M[2, 2]) == 128

        Arblib.add!(M, M, M)
        M2 = TMat(4, 4, prec = 128)
        M2[2, 2] = T(3)
        @test M == M2

        A = TMat([T(i + j) for i = 1:4, j = 1:4])
        @test A[4, 4] == T(8)
        A = TMat([i + j for i = 1:4, j = 1:4]; prec = 96)
        @test A[4, 4] == T(8)
        @test precision(A) == precision(T(8; prec = 96)) == 96

        @test ref(A, 3, 3) isa Union{ArbRef,AcbRef}
    end

    @testset "arithmetic" begin
        AInt = [1 2; 3 4]
        BInt = [5 6; 7 8]
        A = TMat(AInt; prec = 96)
        B = TMat(BInt; prec = 96)
        @test A - B == AInt - BInt
        @test precision(A - B) == 96
        @test (A - B) isa TMat
        @test -B == -BInt
        @test -B isa TMat
        @test -B + A == A - B
        @test precision(-B + A) == 96
        @test A * B == AInt * BInt
        @test A * B isa TMat
    end

    @testset "LinearAlgebra" begin
        A = TMat(rand(3, 3))
        b = TMat(rand(3))
        c = TMat(3, 1)

        # lu factorization
        ldiv!(c, A, b)
        c′ = A \ b
        @test Arblib.overlaps(c, c′) == 1
        ldiv!(c, lu(A), b)
        @test Arblib.overlaps(c, c′) == 1
        ldiv!(c, lu!(copy(A)), b)
        @test Arblib.overlaps(c, c′) == 1
        d = copy(b)
        ldiv!(lu(A), d)
        @test Arblib.overlaps(d, c′) == 1

        # inv
        id = copy(A)
        Arblib.one!(id)
        @test Bool(Arblib.contains(inv(A) * A, id))

        # mul
        A = TMat(rand(3, 2))
        B = TMat(rand(2, 4))
        B_wrong = TMat(rand(3, 4))
        C = TMat(3, 4)
        @test_throws DimensionMismatch("Matrix sizes are not compatible.") A * B_wrong
        LinearAlgebra.mul!(C, A, B)
        @test C == A * B
    end

    @testset "indexing" begin
        A = TMat(3, 1)
        A[3, 1][] = 4
        @test A[3] == A[3, 1]
        B = TMat(reshape(1:15, 3, 5))
        @test B[1:15] == 1:15
        C = TMat(reshape(1:15, 5, 3))
        @test C[1:15] == 1:15
    end

    @testset "similar:" begin
        A = TMat(3, 5; prec = 96)

        a = similar(A)
        @test a isa TMat
        @test precision(a) == precision(A)

        a = similar(A, TRef, (2, 3))
        @test a isa TMat
        @test precision(a) == precision(A)

        for (ElT, VT, MT) in (
            (Arb, ArbVector, ArbMatrix),
            (Acb, AcbVector, AcbMatrix),
            (ArbRef, ArbRefVector, ArbRefMatrix),
            (AcbRef, AcbRefVector, AcbRefMatrix),
        )
            a = similar(A, ElT, 3)
            @test a isa VT
            @test precision(a) == precision(A)

            a = similar(A, ElT, (3, 2))
            @test a isa MT
            @test precision(a) == precision(A)
        end
    end

    @testset "set Rational" begin
        A = TMat(2, 2)
        A[1, 1] = 5 // 8
        @test !iszero(A[1, 1])
    end
end

@testset "MatrixRef: $T" for (T, TRef) in
                             [(ArbMatrix, ArbRefMatrix), (AcbMatrix, AcbRefMatrix)]
    A = T(2, 3; prec = 96)
    A[1, 2] = 3

    B = TRef(A)
    @test B isa TRef
    @test precision(B) == 96
    B[1, 2] = 4
    @test A[1, 2] == 4
    @test B[1, 2] == 4
    C = T(B)
    @test C == A
    @test C[1, 2] == 4
end


text:
"""Basic functions for computing models."""

"""
All foreground models inherit from this type.
"""
abstract type AbstractForegroundModel end

"""
    get_cosmology(::Type{T}; h=0.69, Neff=3.04, OmegaK=0.0,
        OmegaM=0.29, OmegaR=nothing, Tcmb=2.7255, w0=-1, wa=0)

Construct a background cosmology. This function duplicates the cosmology() function
in Cosmology.jl, but with the numeric type specified.

# Arguments:
- `::Type{T}`: numerical type to use for calculations

# Keywords
- `h` - Dimensionless Hubble constant
- `OmegaK` - Curvature density (Ω_k)
- `OmegaM` - Matter density (Ω_m)
- `OmegaR` - Radiation density (Ω_r)
- `Tcmb` - CMB temperature in Kelvin; used to compute Ω_γ
- `Neff` - Effective number of massless neutrino species; used to compute Ω_ν
- `w0` - CPL dark energy equation of state; `w = w0 + wa(1-a)`
- `wa` - CPL dark energy equation of state; `w = w0 + wa(1-a)`

# Example
```julia-repl
julia> get_cosmology(Float32; h=0.7)
Cosmology.FlatLCDM{Float32}(0.7f0, 0.7099147f0, 0.29f0, 8.5307016f-5)
```
"""
function get_cosmology(::Type{T}; h=0.69,
                   Neff=3.04,
                   OmegaK=0.0,
                   OmegaM=0.29,
                   OmegaR=nothing,
                   Tcmb=2.7255,
                   w0=-1,
                   wa=0) where T

    if OmegaR === nothing
        OmegaG = 4.48131e-7*Tcmb^4/h^2
        OmegaN = Neff*OmegaG*(7/8)*(4/11)^(4/3)
        OmegaR = OmegaG + OmegaN
    end

    OmegaL = 1 - OmegaK - OmegaM - OmegaR

    if !(w0 == -1 && wa == 0)
        return Cosmology.WCDM{T}(h, OmegaK, OmegaL, OmegaM, OmegaR, w0, wa)
    end

    if OmegaK < 0
        return Cosmology.ClosedLCDM{T}(h, OmegaK, OmegaL, OmegaM, OmegaR)
    elseif OmegaK > 0
        return Cosmology.OpenLCDM{T}(h, OmegaK, OmegaL, OmegaM, OmegaR)
    else
        return Cosmology.FlatLCDM{T}(h, OmegaL, OmegaM, OmegaR)
    end
end
get_cosmology(; h=0.69, Neff=3.04, OmegaK=0.0, OmegaM=0.29, OmegaR=nothing, Tcmb=2.7255,
    w0=-1, wa=0) where T = get_cosmology(Float32; h=h, Neff=Neff, OmegaK=OmegaK,
        OmegaM=OmegaM, OmegaR=OmegaR, Tcmb=Tcmb, w0=w0, wa=wa)

"""
    build_r2z_interpolator(min_z::T, max_z::T,
        cosmo::Cosmology.AbstractCosmology; n_bins=2000) where T

Construct a fast r2z linear interpolator.
"""
function build_r2z_interpolator(min_z::T, max_z::T,
    cosmo::Cosmology.AbstractCosmology; n_bins=2000) where T

    zrange = LinRange(min_z, max_z, n_bins)
    rrange = zero(zrange)
    for i in 1:n_bins
        rrange[i] = ustrip(T, u"Mpc",
            Cosmology.comoving_radial_dist(u"Mpc", cosmo, zrange[i]))
    end
    r2z = LinearInterpolation(rrange, zrange; extrapolation_bc=Line());
    return r2z
end

"""
    mz2c(m::T, z::T, cosmo::Cosmology.FlatLCDM{T}) where T

Compute concentration factor from Duffy et al. 2008.
"""
function mz2c(m::T, z::T, cosmo::Cosmology.FlatLCDM{T}) where T
    return T(7.85) * (m / (T(2.0e12)/cosmo.h))^(T(-0.081)) / (T(1.0)+z)^T(0.71)
end

"""
    m2r(m::T, cosmo::Cosmology.FlatLCDM{T}) where T

Convert virial mass to virial radius.
"""
function m2r(m::T, cosmo::Cosmology.FlatLCDM{T}) where T
    rho = T(2.78e11) * cosmo.Ω_m * cosmo.h^2
    return (T(3.0/(4π)) * m / rho)^T(1.0/3.0)
end

function f_nfw(x::T) where T
    return log1p(x) - x/(one(T)+x)
end
g(x, c, m) = m - f_nfw(x)/f_nfw(c)

function solve_r(c::T, lnm::T) where T
    m = exp(lnm)
    return find_zero( x -> g(x,c,m) , (T(0.0),T(100.0)), Bisection()) / c
end

"""
    build_c_lnm2r_interpolator(;cmin::T=1f-3, cmax::T=25.0f0,
        mmin::T=-7.1f0, mmax::T=0.0f0, nbin=100) where T

Generate a LinearInterpolation object that turns concentration
and ln(M_halo) into satellite radius.
"""
function build_c_lnm2r_interpolator(;cmin::T=1f-3, cmax::T=25.0f0,
        mmin::T=-7.1f0, mmax::T=0.0f0, nbin=100) where T
    # these defaults imply a minimum fractional mass of exp(-7)

    cs  = LinRange(cmin,cmax,nbin)
    ms  = LinRange(mmin,mmax,nbin)
    r = [solve_r(c_,m_) for c_ in cs, m_ in ms]

    return LinearInterpolation((T.(cs), T.(ms)), T.(r))
end

random_phi(::Type{Float32}) = rand(Float32) * Float32(2π)
random_theta(::Type{Float32}) = acos(2.0f0*rand(Float32)-1.0f0)
random_phi(t::DataType) = rand(t) * t(2π)
random_theta(t::DataType) = acos( t(2.0)*rand(t)-one(t))

"""
Inverse square law with redshift dependence.
"""
function l2f(luminosity::T, r_comoving::T, redshift::T) where T
    return luminosity / (T(4π) * r_comoving^2 * (one(T) + redshift) )
end

function flux2map!(result_map::Map{T_map,RingOrder}, fluxes, theta, phi) where {T_map, T}

    pixel_array = result_map.pixels
    fill!(pixel_array, zero(T))  # prepare the frequency map

    flux_to_Jy = ustrip(u"Jy", 1u"W/Hz*Mpc^-2")

    source_offset_I = generate_subhalo_offsets(sources.nsources_I)
    source_offset_II = generate_subhalo_offsets(sources.nsources_II)
    N_halo = size(sources.halo_mass, 1)
    total_n_sat_I = size(sources.L_I_151, 1)
    total_n_sat_II = size(sources.L_II_151, 1)

    flux_I = Array{T, 1}(undef, total_n_sat_I)
    flux_II = Array{T, 1}(undef, total_n_sat_II)
    redshift_I = Array{T, 1}(undef, total_n_sat_I)
    redshift_II = Array{T, 1}(undef, total_n_sat_II)
    θ_I = Array{T, 1}(undef, total_n_sat_I)
    θ_II = Array{T, 1}(undef, total_n_sat_II)
    ϕ_I = Array{T, 1}(undef, total_n_sat_I)
    ϕ_II = Array{T, 1}(undef, total_n_sat_II)

    Threads.@threads for i_halo in 1:N_halo
        nu = (1 + sources.redshift[i_halo]) * nu_obs
        hp_ind = sources.halo_hp_ind[i_halo]

        # Generate FR I fluxes
        for source_index in 1:sources.nsources_I[i_halo]
            # index of satellite in satellite arrays
            i_sat = source_offset_I[i_halo] + source_index
            L_c, L_l = get_core_lobe_lum(
                sources.L_I_151[i_sat], nu, model.I_R_int, model.I_γ,
                model.a_0,
                sources.a_coeff_I[1,i_sat], sources.a_coeff_I[2,i_sat],
                sources.a_coeff_I[3,i_sat], sources.cosθ_I[i_sat])

            flux_I[i_sat] = l2f(L_c+L_l, sources.dist[i_halo],
                sources.redshift[i_halo]) * flux_to_Jy
            redshift_I[i_sat] = sources.redshift[i_halo]
            θ_I[i_sat] = sources.θ[i_halo]
            ϕ_I[i_sat] = sources.ϕ[i_halo]
            pixel_array[hp_ind] += flux_I[i_sat]
        end

        # Generate FR II fluxes
        for source_index in 1:sources.nsources_II[i_halo]
            # index of satellite in satellite arrays
            i_sat = source_offset_II[i_halo] + source_index
            L_c, L_l = get_core_lobe_lum(
                sources.L_II_151[i_sat], nu, model.II_R_int, model.II_γ,
                model.a_0,
                sources.a_coeff_II[1,i_sat], sources.a_coeff_II[2,i_sat],
                sources.a_coeff_II[3,i_sat], sources.cosθ_II[i_sat])

            flux_II[i_sat] = l2f(L_c+L_l, sources.dist[i_halo],
                sources.redshift[i_halo]) * flux_to_Jy
            redshift_II[i_sat] = sources.redshift[i_halo]
            θ_II[i_sat] = sources.θ[i_halo]
            ϕ_II[i_sat] = sources.ϕ[i_halo]
            pixel_array[hp_ind] += flux_II[i_sat]
        end
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(result_map.resolution.nside)
    pixel_array .*= per_pixel_steradian

    # maps are in Jansky per steradian, fluxes are in Jansky
    if return_fluxes
        return flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II
    end

end


text:
using DEBTrait
using CSV, DataFrames, Statistics
using DifferentialEquations, Plots

n_polymers              = 0
n_monomers              = 6
n_microbes              = 1
n_enzymes               = 1
n_minerals              = 0

df_isolates             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/microbes2traits.csv", DataFrame, missingstring = "N/A")
df_isolates             = filter(row -> row[:Rhizosphere_response] == "positive" || row[:Rhizosphere_response] == "negative", df_isolates)
df_monomers             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_class.csv",DataFrame, missingstring = "N/A")

isolate_id              = collect(1:27)
monomer_id              = collect(1:6)
monomer_concentration   = [1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1.0]
#
BGE                     = zeros(27, 6, 13)
r                       = zeros(27, 6, 13)
m_E                     = zeros(27, 6, 13)
x                       = zeros(27, 6, 13)
#
function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1] < 1e-8
end
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)

i = 1
id_isolate              = isolate_id[i]
#
L_DNA                   = df_isolates.Genome_size[id_isolate]*ones(n_microbes)
rRNA                    = df_isolates.rRNA_genes[id_isolate]*ones(n_microbes)
min_gentime             = df_isolates.Min_gen_time[id_isolate]*ones(n_microbes)
#
z_sugars_p              = df_isolates.z_sugars[id_isolate]*ones(1,1)
z_organics_p            = df_isolates.z_organic_acids[id_isolate]*ones(1,1)
z_aminos_p              = df_isolates.z_amino_acids[id_isolate]*ones(1,1)
z_fattys_p              = df_isolates.z_fatty_acids[id_isolate]*ones(1,1)
z_nucleos_p             = df_isolates.z_nucleotides[id_isolate]*ones(1,1)
z_auxins_p              = df_isolates.z_auxins[id_isolate]*ones(1,1)
genome_distr_p          = vcat(z_sugars_p, z_organics_p, z_aminos_p, z_fattys_p, z_nucleos_p, z_auxins_p)
    #
j = 6
n_monomers              = 1
id_monomer              = monomer_id[j]
Isolate                 = BaseGenome(L_DNA, rRNA, min_gentime, genome_distr_p[id_monomer]*ones(n_monomers,n_microbes))
Monomer                 = BaseMonomer(df_monomers.D_S[id_monomer]*ones(n_monomers), df_monomers.N_C[id_monomer]*ones(n_monomers), df_monomers.N_N[id_monomer]*ones(n_monomers), df_monomers.y_DE[id_monomer]*ones(n_monomers))
Soil                    = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
#
MetabolismC             = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
MetabolismC.y_EM        = ones(n_microbes)
MetabolismC.y_EX        = ones(n_microbes).*1.4
MetabolismC.α           = zeros(n_microbes)

if df_p.Rhizosphere_response[1] .== "positive"
    MetabolismC.α .= 0.01
else
    MetabolismC.α .= 0.3
end
#

AssimilationC           = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
Setup                   = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
p                       = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
#########################
p.n_polymers              = 0
p.n_monomers              = 1
p.n_microbes              = 1
p.n_enzymes               = 1
p.n_minerals              = 0
#
u0 = zeros(p.dim)
        for k in 1:13
            init_D0 = monomer_concentration[k]
            u0[1+p.n_polymers:p.n_polymers+p.n_monomers] = init_D0*ones(p.n_monomers)
            u0[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes] = p.Setup.u0_reserves*1e-6
            u0[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes] = p.Setup.u0_structures*1e-6
            #
            tspan = (0.0,  150.0)
            prob  = ODEProblem(rhs_BatchC!, u0, tspan, p)
            sol   = solve(prob, alg_hints=[:stiff], callback=cb)
            #
            BGE_tmp = [DEBTrait.calc_BGE(zeros(p.dim), sol.u[i], p, 0.0)[1] for i in 1:size(sol.t,1)]
            BGE[i,j,k] = median(BGE_tmp[BGE_tmp.>1e-8])
            r_tmp   = [DEBTrait.calc_growth_rate(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            r[i,j,k] = median(r_tmp)
            E = sol.u[end][1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes]
            V = sol.u[end][1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes]
            m_E_tmp = E./V
            m_E_tmp[m_E_tmp.>=1.0] .= 0.0
            m_E[i,j,k] = m_E_tmp[1]
            x_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[1][1] for i in 1:size(sol.t,1)]
            x[i,j,k] = median(x_tmp)
        end

scatter(monomer_concentration, BGE[1,1,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,1,:], xscale=:log10)

scatter(monomer_concentration, BGE[1,2,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,2,:], xscale=:log10)

scatter(monomer_concentration, BGE[1,3,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,3,:], xscale=:log10)

scatter(monomer_concentration, BGE[1,4,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,4,:], xscale=:log10)

scatter(monomer_concentration, BGE[1,5,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,5,:], xscale=:log10)

scatter(monomer_concentration, BGE[1,6,:], xscale=:log10)
scatter!(monomer_concentration, BGE[2,6,:], xscale=:log10)

scatter(monomer_concentration, r[1,1,:], xscale=:log10)
scatter!(monomer_concentration, r[2,1,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,1,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,1,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,2,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,2,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,3,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,3,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,4,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,4,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,5,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,5,:], xscale=:log10)

scatter(monomer_concentration, m_E[1,6,:], xscale=:log10)
scatter!(monomer_concentration, m_E[2,6,:], xscale=:log10)

scatter(monomer_concentration, x[1,2,:], xscale=:log10)
scatter!(monomer_concentration, x[2,2,:], xscale=:log10)


text:
using JSON
using PipeLayout
using Plots

pyplot()

open(ARGS[1]) do f
    for line in eachline(f)
        key = strip(line)
        println(key)
        inst = PipeLayout.read_instance(".", key)
        topo = PipeLayout.read_topology(".", key)

        plot(inst)
        plot!(topo)
        savefig("$key.png")
    end
end


text:
using Pkg

Pkg.add([
    PackageSpec(
        name = "BenchmarkConfigSweeps",
        url = "https://github.com/tkf/BenchmarkConfigSweeps.jl.git",
    ),
    PackageSpec(
        name = "BenchPerf",
        # path = dirname(@__DIR__),
        url = "https://github.com/tkf/BenchPerf.jl.git",
    ),
    PackageSpec(
        name = "BenchPerfConfigSweeps",
        # path = dirname(@__DIR__),
        url = "https://github.com/tkf/BenchPerf.jl.git",
        subdir = "lib/BenchPerfConfigSweeps",
    ),
])

Pkg.develop([
    PackageSpec(
        name = "BenchPerf",
        path = dirname(@__DIR__),
        # url = "https://github.com/tkf/BenchPerf.jl.git",
    ),
    PackageSpec(
        name = "BenchPerfConfigSweeps",
        path = dirname(@__DIR__),
        # url = "https://github.com/tkf/BenchPerf.jl.git",
        subdir = "lib/BenchPerfConfigSweeps",
    ),
])

Pkg.instantiate()


text:
# ===== VSB Verifications: Slip Velocity =====

println("=== BEGINNING Verifications2.jl ===")

# === Imports ===
println("   IMPORTING PACKAGES...")
using VSB
using SimpleVPM
using LinearAlgebra
using Printf

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

# === Geometry and Grids ===
println("   GENERATING GEOMETRY...")
NPTS = 500                              # Number of Points on Boundary
del_theta = (2*pi)/(NPTS - 1)           # Step in theta around body
R = 1.0                                 # Cylinder Radius
X_coor = zeros(NPTS)                    # X-coordinates
Y_coor = zeros(NPTS)                    # Y-coordinates
X_tHat = zeros(NPTS)                    # X Tangent Vector
Y_tHat = zeros(NPTS)                    # X Tangent Vector
X_nHat = zeros(NPTS)                    # Y Normal Vector
Y_nHat = zeros(NPTS)                    # Y Normal Vector
theta = zeros(NPTS)                     # Angle from positive X-axis
for i=1:NPTS
    X_coor[i] = R*cos((i-1)*del_theta)
    Y_coor[i] = R*sin((i-1)*del_theta)

    X_tHat[i] = sin((i-1)*del_theta)
    Y_tHat[i] = -cos((i-1)*del_theta)

    X_nHat[i] = cos((i-1)*del_theta)
    Y_nHat[i] = sin((i-1)*del_theta)

    theta[i] = (i-1)*del_theta
end
body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]

body = VSB.Boundary(body_pts,t_hats,n_hats)

# === RBF Coefficients and Velocity Functions ===
U_inf = [1.0, 0.0, 0.0]
U_inf_func(X) = U_inf
print("   CALCULATING ALPHAS...")
alphas = VSB.CalcVSCoefs(body, U_slip=U_inf_func)
print("   DONE!")
println()

U_gamma(X) = VSB.CalcVSVelocityCirlce(body, alphas, X)

println("   CALCULATING VELOCITIES...")
println(" -      X     -  |  - U_SLIP TANGENT -  |  - U_SLIP NORMAL -  |  -      U_GAMMA      -  |  - U_GAMMA TANGENT -  |  - U_GAMMA NORMAL -  |")
println("=========================================================================================================================================")
# === Calculate U_slip along surface ===
for i in 1:NPTS
    X_eval = body.bodyPTS[i]
    n_hat = body.nHats[i]
    t_hat = body.tHats[i]

    u_gamma = U_gamma(X_eval)

    u_slip_normal = dot(U_inf_func(X_eval), n_hat)
    u_gamma_normal = dot(U_gamma(X_eval), n_hat)

    u_slip_tangent = dot(U_inf_func(X_eval), t_hat)
    u_gamma_tangent = dot(U_gamma(X_eval), t_hat)

    if i % 10 == 0
        println("  [",@sprintf("%5.2f",X_eval[1]),",",@sprintf("%5.2f",X_eval[2]),"]",
                "  |  ",@sprintf("%18.3f",u_slip_tangent),
                "  |  ",@sprintf("%17.3f",u_slip_normal),
                "  |  [",@sprintf("%5.2f",u_gamma[1]),", ",@sprintf("%5.2f",u_gamma[2]),", ",@sprintf("%5.2f",u_gamma[3]),"]",
                "  |  ",@sprintf("%19.3f",u_gamma_tangent),
                "  |  ",@sprintf("%18.3f",u_gamma_normal),
                "  |")
    end
end
println("=========================================================================================================================================")

println("=== END OF Verifications2.jl ===")


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()
using PackageCompiler
PackageCompiler.create_sysimage([:Plots, :Pluto, :PlutoUI, :MLJ, :Distributions,
                                 :Flux, :Zygote, :MLJLinearModels, :Random,
                                 :NearestNeighborModels, :Statistics, :LinearAlgebra,
                                 :DataFrames, :CSV, :StatsPlots];
                                sysimage_path=joinpath(@__DIR__, "mlcourse.so"),
                                precompile_execution_file=joinpath(@__DIR__, "warmup.jl"))

max_stars_repo_name: ablaom/MLCourse

max_issues_repo_name: ablaom/MLCourse

max_forks_repo_name: ablaom/MLCourse


text:
#ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника. Робот - в произвольной клетке поля между внешней и внутренней перегородками.
#РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.
function mark_whole_border(r)
    x,y = corner_border(r)
    i = find_whole_border(r)
    frame_border(r,i)
    secondcorner(r)
    back_home(r,x,y)
end
function find_whole_border(r)
    c = 1
    cc = 1
    while isborder(r,Ost)==0
        move!(r,Ost)
        c += 1
    end
    move!(r,Nord)
    i=0
    while true
        while isborder(r,West)==0
            move!(r,West)
            cc += 1
        end
        if cc == c
            cc = 1
            move!(r,Nord)
        else
            i=1
            break
        end
        while isborder(r,Ost)==0
            move!(r,Ost)
            cc += 1
        end
        if cc == c
            cc = 1
            move!(r,Nord)
        else
            i = 3
            break
        end
    end
return i
end
function corner_border!(r)
    x=0
    y=0
    while isborder(r,Sud)==0
        move!(r,Sud)
        y+= 1
    end
    while isborder(r,West)==0
        move!(r,West)
        x += 1
    end
    while isborder(r,Sud)==0
        move!(r,Sud)
        y+= 1
    end
    while isborder(r,West)==0
        move!(r,West)
        x += 1
    end
    return x, y
end
function secondcorner(r)
    while isborder(r,Sud)==0
        move!(r,Sud)
    end
    while isborder(r,West)==0
        move!(r,West)
    end
end
function back_home(r,x,y)
    for i in 1:x
        move!(r,Ost)
    end
    for i in 1:y
        move!(r,Nord)
    end
end
function frame_border(r,i)
    putmarker!(r)
    while isborder(r,HorizonSide(i))
        move!(r,HorizonSide(0))
        putmarker!(r)
    end
    move!(r,HorizonSide(i))
    putmarker!(r)
    while isborder(r,HorizonSide(2))
        move!(r,HorizonSide(i))
        putmarker!(r)
    end
    move!(r,HorizonSide(2))
    putmarker!(r)
    i += 2
    while isborder(r,HorizonSide(i % 4))
        move!(r,HorizonSide(2))
        putmarker!(r)
    end
    move!(r,HorizonSide(i % 4))
    putmarker!(r)
    while isborder(r,HorizonSide(0))
        move!(r,HorizonSide(i % 4))
        putmarker!(r)
    end
end
mark_whole_border(r)


text:
using Test, DynamicEnergyBudgets, Unitful, OrdinaryDiffEq, DimensionalData
using Microclimate, Photosynthesis

u0 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]u"mol" # Initial value
du0 = DimArray(fill(0.0u"mol/hr", 6), (X(Val((:VS,:CS,:NS,:VR,:CR,:NR))),))

env = MicroclimControl()
maintenance_cost = -0.001u"hr^-1"
shared=SharedParams(
    resorption_pars=nothing,
    tempcorr_pars=nothing,
    su_pars = ParallelComplementarySU(),
    core_pars = DEBCore(
        j_E_mai = maintenance_cost,
        y_V_E = 1.0,
        y_E_C = 1.0,
        y_E_N = 100.0,
        n_N_V = 0.01,
        n_N_E = 0.01,
    ),
    catabolism_pars = CatabolismCNshared(
        k = 0.5u"d^-1",
    ),
)

@testset "balance with growth and maintenance is correct" begin
    plant1 = Plant(time = 0u"hr":1u"hr":8001u"hr",
        params=(
            Params(
                assimilation_pars = ConstantCAssim(
                    c_uptake = 0.0u"mol*mol^-1*hr^-1",
                ),
            ),
            Params(
                assimilation_pars = ConstantNAssim(
                    n_uptake = 0.0u"mol*mol^-1*hr^-1",
                ),
            ),
        ),
        shared=shared,
        environment=env
    )

    du1 = deepcopy(du0)
    u1 = deepcopy(u0)
    plant1(du1, u1, nothing, 1u"hr")

    @test du1[:VS] > 0.0u"mol/hr"
    @test du1[:CS] < 0.0u"mol/hr"
    @test du1[:NS] < 0.0u"mol/hr"
    @test du1[:VR] > 0.0u"mol/hr"
    @test du1[:CR] < 0.0u"mol/hr"
    @test du1[:NR] < 0.0u"mol/hr"

    # Carbon is balanced, with losses to maintenance
    @test du1[:VS] == -du1[:CS] - maintenance_cost * 1u"mol"
    @test du1[:VR] == -du1[:CR] - maintenance_cost * 1u"mol"
    @test du1[:VS] ≈ -100du1[:NS] - maintenance_cost * 1u"mol"
    @test du1[:VR] ≈ -100du1[:NR] - maintenance_cost * 1u"mol"
end

@testset "balance with assimilation is correct" begin
    plant = Plant(time = 0u"hr":1u"hr":8001u"hr",
        params=(
            Params(
                assimilation_pars = ConstantCAssim(
                    c_uptake = 1.0u"mol*mol^-1*hr^-1",
                ),
            ),
            Params(
                assimilation_pars = ConstantNAssim(
                    n_uptake = 1.0u"mol*mol^-1*hr^-1",
                ),
            ),
        ),
        shared=shared,
        environment=env,
    )

    du2 = deepcopy(du0)
    u2 = deepcopy(u0)
    plant(du2, u2, nothing, 1u"hr")

    # Carbon is balanced, with added C and N from assimilation
    @test du2[:VS] ≈ 1.0u"mol/hr" - du2[:CS] - maintenance_cost * 1u"mol"
    @test du2[:VR] ≈ -du2[:CR] - maintenance_cost * 1u"mol"
    @test du2[:VS]/100 ≈ -du2[:NS] - maintenance_cost * 1u"mol"/100
    @test du2[:VR]/100 ≈ 1.0u"mol/hr" - du2[:NR] - maintenance_cost * 1u"mol"/100
end

@testset "balance with resorption is correct" begin
    plant3 = Plant(time = 0u"hr":1u"hr":8001u"hr",
        params=(
            Params(
                assimilation_pars = ConstantCAssim(
                    c_uptake = 0.0u"mol*mol^-1*hr^-1",
                ),
            ),
            Params(
                assimilation_pars = ConstantNAssim(
                    n_uptake = 0.0u"mol*mol^-1*hr^-1",
                ),
            ),
        ),
        shared=SharedParams(
            resorption_pars=LosslessResorption(),
            tempcorr_pars=nothing,
            su_pars = ParallelComplementarySU(),
            core_pars = DEBCore(
                j_E_mai = maintenance_cost,
                y_V_E = 1.0,
                y_E_C = 1.0,
                y_E_N = 100.0,
                n_N_V = 0.01,
                n_N_E = 0.01,
            ),
            catabolism_pars = CatabolismCNshared(
                k = 0.5u"d^-1",
            ),
        ),
        environment=env,
    )

    du3 = deepcopy(du0)
    u3 = deepcopy(u0)
    plant3(du3, u3, nothing, 1u"hr")

    # Carbon is balanced, with losses to maintenance and returned resorption
    @test du3[:VS] == -du3[:CS] - maintenance_cost * 1u"mol"
    @test du3[:VR] == -du3[:CR] - maintenance_cost * 1u"mol"
    @test du3[:VS] ≈ -100du3[:NS] - maintenance_cost * 1u"mol"
    @test du3[:VR] ≈ -100du3[:NR] - maintenance_cost * 1u"mol"
end


text:
using Wavelets
using Base.Test

print("util ...\n")

@test detailn(0) == 1
@test detailn(1) == 2
@test detailindex(0,1) == 2

j = 5
@test detailn(j)+1 == detailrange(j)[1]
@test detailindex(j,3) == detailrange(j)[3]
J = 7
@test nscales(2^J) == J
@test detailn(maxlevel(2^J)) == 2^(J-1)
@test nscales(rand(2^J)) == J
@test maxlevel(rand(2^J)) == J-1

@test tl2level(512,1) == 8
@test level2tl(512,9) == 0
@test tl2level(rand(2^9),1) == 8
@test level2tl(rand(2^9),9) == 0

# UTILITY FUNCTIONS

@test mirror([1]) == [1]
@test mirror([2,3]) == [2,-3]
@test mirror([2,3,4]) == [2,-3,4]
@test mirror([4.9,5,6,7]) == [4.9,-5,6,-7]

@test upsample([1]) == [1,0]
@test upsample([1],1) == [0,1]
@test upsample([1,2]) == [1,0,2,0]
@test upsample([1,2],1) == [0,1,0,2]
@test downsample([1,2]) == [1]
@test downsample([1,2],1) == [2]
@test downsample([1,2,3,4]) == [1,3]
@test downsample([1,2,3,4],1) == [2,4]

@test iscube(rand(4)) == true
@test iscube(rand(4,4,4,4)) == true
@test iscube(rand(4,5)) == false
@test isdyadic(rand(4)) == true
@test isdyadic(rand(32,32)) == true
@test isdyadic(rand(6)) == false

n = 64
@test wcount(randn(n)) == n
@test wcount(rand(n),0) == n
@test wcount(rand(n),1.01) == 0
@test wcount(rand(n).-2.5,1.01) == n
@test wcount([0,2,6,7,8],6.6) == 2
@test wcount([10,-11,6,7,8.0,-5,-8,0],3,level=2) == 3
@test wcount(randn(n,n)) == n*n
@test wcount(rand(n,n),1.01) == 0
@test wcount([-1 2;3 4],2.5) == 2

n = 32
x = randn(n)
y = randn(n)
for sh in (-43, -5, 0, 1, 28, 41)
    @test circshift(x,sh) == circshift!(copy(x),sh)
    @test circshift(x,-sh) == circshift!(copy(x),-sh)
    @test x == circshift!(circshift!(copy(x),sh),-sh)
    @test circshift(x,sh) == circshift!(similar(x),x,sh)
    @test circshift(x,-sh) == circshift!(similar(x),x,-sh)
end

a0 = [1,2]
exp = [1,2]
n = 2
@test split!(copy(a0)) == exp
@test split!(copy(a0), n, similar(a0)) == exp
@test split!(similar(a0), copy(a0), n) == exp
@test merge!(copy(a0)) == exp
@test merge!(copy(a0), n, similar(a0)) == exp
@test merge!(similar(a0), copy(a0), n) == exp

a0 = [-1,2,3.3,4]
exp = [-1,3.3,2,4]
n = 4
@test split!(copy(a0)) == exp
@test split!(copy(a0), n, similar(a0)) == exp
@test split!(similar(a0), copy(a0), n) == exp
a0 = [-1,3.3,2,4]
exp = [-1,2,3.3,4]
@test merge!(copy(a0)) == exp
@test merge!(copy(a0), n, similar(a0)) == exp
@test merge!(similar(a0), copy(a0), n) == exp

a0 = randn(128)
@test a0 == merge!(split!(a0))

# with strides
a0 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
ia = 3
inca = 2
@test split!(similar(a0), copy(a0), ia, inca, 2)[1:2] == [3,5]
@test split!(similar(a0), copy(a0), ia, inca, 4)[1:4] == [3,7,5,9]
@test split!(similar(a0), copy(a0), 1, 3, 4)[1:4] == [1,7,4,10]
@test split!(similar(a0), copy(a0), 1, 1, 8)[1:8] ==  split!(similar(a0), copy(a0), 8)[1:8]

@test merge!(similar(a0), ia, inca, copy(a0), 2)[ia:inca:ia+(2-1)*inca] == [1,2]
@test merge!(similar(a0), ia, inca, copy(a0), 4)[ia:inca:ia+(4-1)*inca] == [1,3,2,4]
ia = 1
inca = 3
@test merge!(similar(a0), ia, inca, copy(a0), 4)[ia:inca:ia+(4-1)*inca] == [1,3,2,4]
@test merge!(similar(a0), 1, 1, copy(a0), 8)[1:8] ==  merge!(similar(a0), copy(a0), 8)[1:8]

n = 128
x = randn(n)
for L = 0:nscales(n)
    for st in (:full, :dwt)
        @test isvalidtree(x, maketree(n, L, st))
    end
end
tree = maketree(n, 4)
tree[5] = false
@test !isvalidtree(x, tree)
tree = maketree(n, 4)
tree[18] = true
@test isvalidtree(x, tree)
tree[7] = false
@test !isvalidtree(x, tree)
tree[7] = true
tree[9] = false
@test !isvalidtree(x, tree)

EE = Exception
@test_throws EE maketree(n, 4, :foo)


text:
using Test
using AutoPages: split_by_camel_case, gather_pages, replace_reverse

@testset "Split by CamelCase" begin
    @test split_by_camel_case("AlTeRnAtInG") == ["Al", "Te", "Rn", "At", "In", "G"]

    @test split_by_camel_case("CamelCase") == ["Camel", "Case"]

    @test split_by_camel_case("UPPERThenLowerCase") == ["UPPER", "Then", "Lower", "Case"]

    @test split_by_camel_case("ALLCAPS") == ["ALLCAPS"]

    @test split_by_camel_case("MixOfCASES") == ["Mix", "Of", "CASES"]

    @test split_by_camel_case("LowercasethenUPPERCASE") == ["Lowercasethen", "UPPERCASE"]

    @test split_by_camel_case("DGMethods") == ["DG", "Methods"]
end

@testset "Replace reverse" begin
    @test replace_reverse("abc_abc", "a"=>"x"; count=1) == "abc_xbc"
    @test replace_reverse("abc_abc", "a"=>"x"; count=2) == "xbc_xbc"
    @test replace_reverse("abc_abc", "a"=>"x") == "xbc_xbc"
end

function tutorials_list(leading_path_sep)
  return String[
  leading_path_sep * joinpath("path_to", "file1.jl"),
  leading_path_sep * joinpath("path_to", "file2.jl"),
  leading_path_sep * joinpath("path_to", "SubPath1", "SUBPath2", "file3.jl"),
  leading_path_sep * joinpath("path_to", "SubPath1", "SUBPath2", "file4.jl"),
  leading_path_sep * joinpath("path_to", "SubPath1", "SUBPath3", "file5.jl"),
  leading_path_sep * joinpath("path_to", "SubPath1", "SUBPath3", "file6.jl"),
  ]
end

@testset "Pages Array" begin
    for r in ["", Base.Filesystem.path_separator]
        tutorials = tutorials_list(r)

        arr, _ = gather_pages(;filenames=tutorials)

        arr_base = arr[1].second
        @test arr[1].first == "path_to"

        @test arr_base[1].first == "File1"
        @test arr_base[1].second == joinpath("path_to", "file1.jl")
        @test arr_base[2].first == "File2"
        @test arr_base[2].second == joinpath("path_to", "file2.jl")

        @test arr_base[3].first == "Sub Path1"
        @test arr_base[3].second[1].first == "SUB Path2"
        @test arr_base[3].second[2].first == "SUB Path3"

        @test arr_base[3].second[1].second[1].first == "File3"
        @test arr_base[3].second[1].second[2].first == "File4"
        @test arr_base[3].second[2].second[1].first == "File5"
        @test arr_base[3].second[2].second[2].first == "File6"

        @test arr_base[3].second[1].second[1].second == joinpath("path_to", "SubPath1", "SUBPath2", "file3.jl")
        @test arr_base[3].second[1].second[2].second == joinpath("path_to", "SubPath1", "SUBPath2", "file4.jl")
        @test arr_base[3].second[2].second[1].second == joinpath("path_to", "SubPath1", "SUBPath3", "file5.jl")
        @test arr_base[3].second[2].second[2].second == joinpath("path_to", "SubPath1", "SUBPath3", "file6.jl")
    end
end

@testset "Pages Array Remove First Level" begin
    for r in ["", Base.Filesystem.path_separator]
        tutorials = tutorials_list(r)

        arr, _ = gather_pages(;filenames=tutorials, remove_first_level = true)

        arr_base = arr

        @test arr_base[1].first == "File1"
        @test arr_base[1].second == joinpath("path_to", "file1.jl")
        @test arr_base[2].first == "File2"
        @test arr_base[2].second == joinpath("path_to", "file2.jl")

        @test arr_base[3].first == "Sub Path1"
        @test arr_base[3].second[1].first == "SUB Path2"
        @test arr_base[3].second[2].first == "SUB Path3"

        @test arr_base[3].second[1].second[1].first == "File3"
        @test arr_base[3].second[1].second[2].first == "File4"
        @test arr_base[3].second[2].second[1].first == "File5"
        @test arr_base[3].second[2].second[2].first == "File6"

        @test arr_base[3].second[1].second[1].second == joinpath("path_to", "SubPath1", "SUBPath2", "file3.jl")
        @test arr_base[3].second[1].second[2].second == joinpath("path_to", "SubPath1", "SUBPath2", "file4.jl")
        @test arr_base[3].second[2].second[1].second == joinpath("path_to", "SubPath1", "SUBPath3", "file5.jl")
        @test arr_base[3].second[2].second[2].second == joinpath("path_to", "SubPath1", "SUBPath3", "file6.jl")
    end
end


text:
#standard theme (to make plots consistent and allow for more control for publication-quality plots)
function BAtheme(;args...)
    #font = "'PT Sans','Helvetica Neue','Helvetica',sans-serif"
    #font = "Computer Modern Math"
    font = "'Latin Modern Math','Latin-Modern',serif"
    #font = "Latin Modern"
    #font = "Helvetica"
    #font = "Times New Roman"
    return Theme(line_width = 2pt,# Width of lines in the line geometry. (Measure)
                minor_label_font = font,#: Font used for minor labels such as guide entries and labels. (String)
                #minor_label_font_size: Font size used for minor labels. (Measure)
                major_label_font = font,# Font used for major labels such as titles and axis labels. (String)
                major_label_font_size = 12pt,# Font size used for major labels. (Measure)
                key_title_font = font,# Font used for titles of keys. (String)
                key_title_font_size = 11pt, # Font size used for key titles. (Measure)
                key_label_font = font,# Font used for key entry labels. (String)
                key_label_font_size = 10pt,# Font size used for key entry labels. (Measure)
                bar_spacing = 1pt,# Spacing between bars in Geom.bar. (Measure)
                #boxplot_spacing: Spacing between boxplots in Geom.boxplot. (Measure)
                ;args...
                )
end

#standard continous color scale used by visualization functions in this file
function BAcontinuouscolorscale(scalename::AbstractString)
    #you can pick your own N colors for the gradient with a specified color
    #(or an array of colors) as a seed (i.e. these colors will be included)
    #the other colors will be maximally distinguishable: search-range can be specified:
    #  lchoices - from 0-100
    #  cchoices - from 0-100
    #  hchoices - from 0-360
    #-----------------------
    #colors = distinguishable_colors(2,[color("darkblue"),color("white")],
    #lchoices = linspace(0, 100, 15),cchoices = linspace(0, 100, 15),hchoices = linspace(0, 360, 20))

    #alternatively, use a built-in colormap; they have been designed (scientifically) for
    #most accurately displaying data using a color-coding.
    colors = colormap(scalename)

    #you can transform the colors to simulate certain visual deficiencies
    #with an additional (float) parameter, you could set the strength of the deficiency (from 0.0 to 1.0)
    #colors = deuteranopic(colors) #green-deficient, most common
    #colors = protanopic(colors)   #red-deficient
    #colors = tritanopic(colors)   #blue-yellow deficiency, least common
    return Scale.lab_gradient(colors...)
end

#standard color scale for visualizing matrices)
function BAmatrixvisscale()
    return Scale.ContinuousColorScale(BAcontinuouscolorscale("Blues"))
end

#standard color scale for visualizing probabilities (i.e. values ∈ (0,1))
function BAprobabilityvisscale()
    return Scale.ContinuousColorScale(BAcontinuouscolorscale("Reds"),minvalue=0.0,maxvalue=1)
end

#standard discrete color-scale used by visualization functions in this file
function BAdiscretecolorscale(ncolors::Integer)
    return BAdiscretecolorscale(ncolors,0)
end

#standard discrete color-scale used by visualization functions in this file
#offset is the number of colors that are skipped
function BAdiscretecolorscale(ncolors::Integer, offset::Integer)
    if ncolors < 1
        error("Less than 1 colors requested - invalid operation.")
    end
    if offset < 0
        error("Offset value must be positive but negative value was provided.")
    end
    #these are the Gadfly default colors (for continuous colors)
    colors = [Scale.color_continuous().f(p) for p in linspace(0, 1, (ncolors+offset))]
    #colors = ["blue","purple"]
    navailable = size(colors,1)

    if (offset+ncolors)>navailable
        error("More colors requested than available - extend the manual color scale in BAdiscretecolorscale(...) if possible")
    else
        return Scale.color_discrete_manual(colors[(offset+1):(offset+ncolors)]...)
    end
end

#rectbin plot of p(a) (the marginal)
function visualizeBAmarginal(pa_df::DataFrame, avec::Vector; alabel="Action a", legendlabel="p*(a)", theme_args...)
    #check if strings are provided (the check below is a bit ugly,
    #but there seems to be a bug/problem with [:a_sting] in names(pa_df))
    use_strings = true
    if sum([:a_string].==names(pa_df)) == 0
        use_strings = false
    end

    #do the plotting - using a rectbin plot with discrete scales on both axes
    av = vec(zeros(size(pa_df,1),1))
    if(use_strings)
        plt = plot(pa_df, x=av, y="a_string", color="p_a", Geom.rectbin,
                   Scale.x_discrete, Scale.y_discrete,
                   Guide.colorkey(legendlabel),
                   Guide.xticks(label=false), Guide.xlabel(nothing, orientation=:horizontal),
                   Guide.ylabel(alabel, orientation=:vertical),
                   BAtheme(;theme_args...), BAprobabilityvisscale() )
    else
        plt = plot(pa_df, x=av, y="a", color="p_a", Geom.rectbin,
                   Scale.x_discrete, Scale.y_discrete,
                   Guide.colorkey(legendlabel),
                   Guide.xticks(label=false), Guide.xlabel(nothing, orientation=:horizontal),
                   Guide.ylabel(alabel, orientation=:vertical),
                   BAtheme(;theme_args...), BAprobabilityvisscale() )
    end

    return plt
end

#2D rectbin plot of p(a) (the marginal)
function visualizeBAmarginal(pa::Vector, avec::Vector; alabel="Action a", legendlabel="p*(a)", theme_args...)
    pa_df = BAmarginal2DataFrame(pa,avec)
    plt = visualizeBAmarginal(pa_df, avec, alabel=alabel, legendlabel=legendlabel; theme_args...)
    return plt
end

#2D rectbin plot of p(a) (the marginal)
function visualizeBAmarginal{T<:AbstractString}(pa::Vector, avec::Vector, a_strings::Vector{T};
                             alabel="Action a", legendlabel="p*(a)", theme_args...)

    pa_df = BAmarginal2DataFrame(pa,avec,a_strings)
    plt = visualizeBAmarginal(pa_df, avec, alabel=alabel, legendlabel=legendlabel; theme_args...)
    return plt
end

#2D rectbin plot of p(a|w) (the conditional)
function visualizeBAconditional(pagw_df::DataFrame, avec::Vector, wvec::Vector;
                                alabel="Action a", wlabel="World state o", legendlabel="p*(a|w)",
                                useprob_colorscale::Bool=true, theme_args...)
    #check if strings are provided (the check below is a bit ugly,
    #but there seems to be a bug/problem with [:a_sting] in names(pagw_df))
    #if either one of the strings is missing, don't use both -
    #only providing a string-representation for one variable is not provided
    use_strings = true
    if sum([:a_string].==names(pagw_df)) == 0
        use_strings = false
    end
    if sum([:w_string].==names(pagw_df)) == 0
        use_strings = false
    end

    if useprob_colorscale
        colorscale = BAprobabilityvisscale()
    else
        colorscale = BAmatrixvisscale()
    end

    #do the plotting - using a rectbin plot with discrete scales on both axes
    if(use_strings)
        plt = plot(pagw_df, x="w_string", y="a_string", color="p_agw", Geom.rectbin,
                   Scale.x_discrete, Scale.y_discrete,
                   Guide.colorkey(legendlabel),Guide.xticks(orientation=:vertical),
                   Guide.xlabel(wlabel, orientation=:horizontal), Guide.ylabel(alabel, orientation=:vertical),
                   BAtheme(;theme_args...), colorscale )
    else
        plt = plot(pagw_df, x="w", y="a", color="p_agw", Geom.rectbin,
                   Scale.x_discrete, Scale.y_discrete,
                   Guide.colorkey(legendlabel),
                   Guide.xlabel(wlabel, orientation=:horizontal), Guide.ylabel(alabel, orientation=:vertical),
                   BAtheme(;theme_args...), colorscale )
    end

    return plt
end

#2D rectbin plot of p(a|w) (the conditional)
function visualizeBAconditional(pagw::Matrix, avec::Vector, wvec::Vector;
                                alabel="Action a", wlabel="World state w", legendlabel="p*(a|w)",
                                useprob_colorscale::Bool=true, theme_args...)

    pagw_df = BAconditional2DataFrame(pagw,avec,wvec)
    plt = visualizeBAconditional(pagw_df, avec, wvec, alabel=alabel, wlabel=wlabel,
                                 legendlabel=legendlabel, useprob_colorscale=useprob_colorscale; theme_args...)
    return plt
end

#2D rectbin plot of p(a|w) (the conditional)
function visualizeBAconditional{T1<:AbstractString, T2<:AbstractString}(pagw::Matrix, avec::Vector, wvec::Vector,
                                a_strings::Vector{T1}, w_strings::Vector{T2};
                                alabel="Action a", wlabel="World state w", legendlabel="p*(a|w)",
                                useprob_colorscale::Bool=true, theme_args...)

    pagw_df = BAconditional2DataFrame(pagw,avec,wvec,a_strings,w_strings)
    plt = visualizeBAconditional(pagw_df, avec, wvec, alabel=alabel, wlabel=wlabel,
                                 legendlabel=legendlabel, useprob_colorscale=useprob_colorscale; theme_args...)
    return plt
end

#visualize arbitraty 2D matrix
function visualizeMatrix(M_xy::Matrix, xvec::Vector, yvec::Vector;
                         xlabel="x", ylabel="y", legendlabel="", theme_args...)

    return visualizeBAconditional(M_xy, yvec, xvec, alabel=ylabel, wlabel=xlabel, legendlabel=legendlabel,
                                  useprob_colorscale=false; theme_args...)
end

#visualize arbitraty 2D matrix
function visualizeMatrix{T1<:AbstractString, T2<:AbstractString}(M_xy::Matrix, xvec::Vector, yvec::Vector,
                                x_strings::Vector{T1}, y_strings::Vector{T2};
                                xlabel="x", ylabel="y", legendlabel="", theme_args...)

    return visualizeBAconditional(M_xy, yvec, xvec, y_strings, x_strings, alabel=ylabel, wlabel=xlabel,
                                 legendlabel=legendlabel, useprob_colorscale=false; theme_args...)
end

#plot a 2-dim slice of a 3-dim matrix like p(a|o,w=wk) with a dropdown box that allows to select wk
#
#legendlabelfunc is a function that takes the current w (Integer) and returns a string-representation
#for the legend - the default below simply prints "p(a|o,w=3)" in case w_vis=3 (i.e. the third entry in
#the dropdown was selected)
function visualizeBA_double_conditional(pagow, avec::Vector, ovec::Vector;
                                        alabel="Action a", olabel="Observation o",
                                        legendlabelfunc=(w)->("p*(a|o,w=$w)"),
                                        dropdownlabel = "World state w",
                                        useprob_colorscale::Bool=true, theme_args...)
    if useprob_colorscale
        colorscale = BAprobabilityvisscale()
    else
        colorscale = BAmatrixvisscale()
    end

    #create a dropdown box for selecting wk in p(a|o,w=wk)
    w_vals = [1:size(pagow,3)]  #all valid w-values

    #dropdown box for selecting w
    dropdown_w = dropdown(w_vals,label=dropdownlabel)

    #use lift to connect the actual plotting-code to the dropdown box
    plt_pagow_vis = lift(w_vis->begin
                                pagow_vis = pagow[:,:,w_vis]
                                visualizeBAconditional(pagow_vis, avec, ovec,
                                alabel=alabel, wlabel=olabel, legendlabel=legendlabelfunc(w_vis),
                                useprob_colorscale=useprob_colorscale; theme_args...)
                                end, dropdown_w)

    return dropdown_w, plt_pagow_vis

end

#plot a 2-dim slice of a 3-dim matrix like p(a|o,w=wk) with a dropdown box that allows to select wk
#
#legendlabelfunc is a function that takes the current w (Integer) and returns a string-representation
#for the legend - the default below simply prints "p(a|o,w=3)" in case w_vis=3 (i.e. the third entry in
#the dropdown was selected)
function visualizeBA_double_conditional{T1<:AbstractString, T2<:AbstractString, T3<:AbstractString}(pagow, avec::Vector, ovec::Vector,
                                        a_strings::Vector{T1}, o_strings::Vector{T2}, w_strings::Vector{T3};
                                        alabel="Action a", olabel="Observation o",
                                        legendlabelfunc=(w)->("p*(a|o,w=$w)"),
                                        dropdownlabel = "World state w",
                                        useprob_colorscale::Bool=true, theme_args...)
    if useprob_colorscale
        colorscale = BAprobabilityvisscale()
    else
        colorscale = BAmatrixvisscale()
    end

    #create a dropdown box for selecting wk in p(a|o,w=wk)
    #w_vals = [1:size(pagow,3)]  #all valid w-values

    #dropdown box for selecting w
    dropdown_w = dropdown(w_strings,label=dropdownlabel)

    #use lift to connect the actual plotting-code to the dropdown box
    plt_pagow_vis = lift(w_vis->begin
                                w_index = find(w_strings.==w_vis)[1] #use this syntax to get a scalar index,
                                                                     #otherwise it will be a 1-element array
                                pagow_vis = pagow[:,:,w_index]
                                visualizeBAconditional(pagow_vis, avec, ovec, a_strings, o_strings,
                                alabel=alabel, wlabel=olabel, legendlabel=legendlabelfunc(w_vis),
                                useprob_colorscale=useprob_colorscale; theme_args...)
                                end, dropdown_w)

    return dropdown_w, plt_pagow_vis
end

#visualization of both the marginal and the conditional
function visualizeBAsolution(pa, pagw, avec::Vector, wvec::Vector;
                             alabel="Action a", wlabel="World state w",
                             legendlabel_marginal="p*(a)", legendlabel_conditional="p*(a|w)", suppress_vis::Bool=false,
                             theme_args...)

    plt_marg = visualizeBAmarginal(pa, avec, alabel=alabel, legendlabel=legendlabel_marginal; theme_args...)
    plt_cond = visualizeBAconditional(pagw, avec, wvec, alabel=alabel, wlabel=wlabel,
                                      legendlabel=legendlabel_conditional; theme_args...)

    if suppress_vis == false
        plt_stack = hstack(plt_marg, plt_cond)
        display(plt_stack)
    end

    return plt_marg, plt_cond
end

#visualization of both the marginal and the conditional
function visualizeBAsolution{T1<:AbstractString, T2<:AbstractString}(pa::Vector, pagw::Matrix, avec::Vector, wvec::Vector,
                             a_strings::Vector{T1}, w_strings::Vector{T2};
                             alabel="Action a", wlabel="World state w",
                             legendlabel_marginal="p*(a)", legendlabel_conditional="p*(a|w)", suppress_vis::Bool=false,
                             theme_args...)

    plt_marg = visualizeBAmarginal(pa, avec, a_strings, alabel=alabel, legendlabel=legendlabel_marginal; theme_args...)
    plt_cond = visualizeBAconditional(pagw, avec, wvec, a_strings, w_strings, alabel=alabel, wlabel=wlabel,
                                      legendlabel=legendlabel_conditional; theme_args...)

    if suppress_vis == false
        plt_stack = hstack(plt_marg, plt_cond)
        display(plt_stack)
    end

    return plt_marg, plt_cond
end

#plots the different prob-distributions (this one is quite restrictive in terms of strings for labels, etc.)
function visualize_three_var_BAsolution{T1<:AbstractString, T2<:AbstractString, T3<:AbstractString}(po::Vector, pa::Vector,
                                        pogw::Matrix, pago::Matrix, pagow, pagw::Matrix,
                                        ovec::Vector, avec::Vector, wvec::Vector,
                                        o_strings::Vector{T1}, a_strings::Vector{T2}, w_strings::Vector{T3};
                                        olabel_string="o", alabel_string="a", wlabel_string="w", theme_args...)

    plt_po = visualizeBAmarginal(po, ovec, o_strings, alabel="Observation $olabel_string",
                                 legendlabel="p*($olabel_string)"; theme_args...)

    plt_pa = visualizeBAmarginal(pa, avec, a_strings, alabel="Action $alabel_string",
                                 legendlabel="p*($alabel_string)"; theme_args...)

    plt_pogw = visualizeBAconditional(pogw, ovec, wvec, o_strings, w_strings,
                                      alabel="Observation $olabel_string", wlabel="World state $wlabel_string",
                                      legendlabel="p*($olabel_string|$wlabel_string)"; theme_args...)

    plt_pago = visualizeBAconditional(pago, avec, ovec, a_strings, o_strings,
                                      alabel="Action $alabel_string", wlabel="Observation $olabel_string",
                                      legendlabel="p*($alabel_string|$olabel_string)"; theme_args...)

    plt_pagw = visualizeBAconditional(pagw, avec, wvec, a_strings, w_strings,
                                      alabel="Action $alabel_string", wlabel="World state $wlabel_string",
                                      legendlabel="p*($alabel_string|$wlabel_string)"; theme_args...)

    dpdown, plt_pagow_vis = visualizeBA_double_conditional(pagow, avec, ovec, a_strings, o_strings, w_strings,
                                                           alabel="Action $alabel_string", olabel="Observation $olabel_string",
                                                           legendlabelfunc=(w)->("p*($alabel_string|$olabel_string,$wlabel_string=$w)"),
                                                           dropdownlabel = "World state $wlabel_string"; theme_args...)

    return plt_po, plt_pa, plt_pogw, plt_pago, plt_pagw, dpdown, plt_pagow_vis

end

#plots the evolution of I(A;W), H(A), H(A|W), E[U] and the rate distortion objective as a function of β
function plotperformancemeasures(I::Vector, Ha::Vector, Hagw::Vector, EU::Vector, RDobj::Vector, β_vals::Vector;
                                 suppress_vis::Bool=false, xlabel_perf="β", theme_args...)
    #turn results into data frame
    perf_res = performancemeasures2DataFrame(I, Ha, Hagw, EU, RDobj)
    return plotperformancemeasures(perf_res, β_vals, suppress_vis = suppress_vis; theme_args...)
end

#plots the evolution of I(A;W), H(A), H(A|W), E[U] and the rate distortion objective as a function of β
function plotperformancemeasures(perf_dataframe::DataFrame, β_vals::Vector; suppress_vis::Bool=false,
                                 xlabel_perf="β", theme_args...)
    #append inv. temp. column to data frame
    perf_dataframe[:β] = β_vals

    #------- plot evolution of I(A;W), H(A), H(A|W), E[U(a,w)] and E[U(a,w)]-1/β I(A;W)
    #realign columns of data frame for easy plotting
    #(each column will become a line in the plot - values will depend on the β column)
    #entropic variables (in bits) go in one plot
    perf_res_entropic = stack(perf_dataframe, [:I_aw, :H_a, :H_agw], :β)
    #expected utility and the overall objective (in utils) go in another plot
    perf_res_utils = stack(perf_dataframe, [:E_U, :RD_obj], :β)

    #since Gadfly (currently) does not support setting the legend-strings only (colorkey strings)
    #add a column to the data-frame that has a nice string-representation of the variable-name
    nrows = size(perf_res_entropic,1)
    perf_res_entropic[:variable_str] = ["" for x in 1:nrows]
    perf_res_entropic[perf_res_entropic[:variable].==:I_aw,:variable_str] = "I(A;W)"
    perf_res_entropic[perf_res_entropic[:variable].==:H_a,:variable_str] = "H(A)"
    perf_res_entropic[perf_res_entropic[:variable].==:H_agw,:variable_str] = "H(A|W)"

    nrows = size(perf_res_utils,1)
    perf_res_utils[:variable_str] = ["" for x in 1:nrows]
    perf_res_utils[perf_res_utils[:variable].==:E_U,:variable_str] = "E[U]"
    perf_res_utils[perf_res_utils[:variable].==:RD_obj,:variable_str] = "RU_obj"

    #create the two plots
    plt_entropic = plot(perf_res_entropic,x="β",y="value",color="variable_str",Geom.line,
    BAtheme(key_position = :bottom; theme_args...), Guide.ylabel("[bits]"),Guide.colorkey(""))

    plt_utils = plot(perf_res_utils,x="β",y="value",color="variable_str",Geom.line,
    BAtheme(key_position = :bottom; theme_args...), Guide.ylabel("[utils]"),Guide.colorkey(""),BAdiscretecolorscale(2))
    #----------------------------------------------------------

    #------- plot rate-utility curve (infeasible regsion is shaded)
    nvals = size(perf_dataframe,1)
    ymax_val = maximum(perf_dataframe[:E_U])
    ymax = ones(nvals)*ymax_val
    plt_rateutility = plot(perf_dataframe,x="I_aw",y="E_U", ymin="E_U", ymax=ymax, Geom.line, Geom.ribbon,
    Guide.xlabel("I(A;W) [bits]"),Guide.ylabel("E[U]"), BAtheme(;theme_args...))
    #----------------------------------------------------------

    if suppress_vis == false
        plt_performance = vstack(plt_entropic, plt_utils) #stack plots vertically
        display(plt_performance)
        display(plt_rateutility)
    end

    return plt_entropic, plt_utils, plt_rateutility
end

#plots evolution across BA iterations (i.e. convergence)
function plot_three_var_BA_convergence(performance_df::DataFrame; suppress_vis::Bool=false,
                                       xlabel_perf="Iteration", theme_args...)
    iter = collect(1:size(performance_df,1))

    #add a column with the iteration index to the performance data frame
    performance_df[:iteration] = iter

    #------- plot evolution of I(O;W), I(A;O), I(A;W|O), I(A|W)
    #realign columns of data frame for easy plotting
    #(each column will become a line in the plot - values will depend on the iter column)
    #entropic variables (in bits) go in one plot
    perf_res_MI = stack(performance_df, [:I_ow, :I_ao, :I_awgo, :I_aw], :iteration)
    #expected utility and the overall objective (in utils) go in another plot
    perf_res_utils = stack(performance_df, [:E_U, :Objective_value], :iteration)

    #since Gadfly (currently) does not support setting the legend-strings only (colorkey strings)
    #add a column to the data-frame that has a nice string-representation of the variable-name
    nrows = size(perf_res_MI,1)
    perf_res_MI[:variable_str] = ["" for x in 1:nrows]
    perf_res_MI[perf_res_MI[:variable].==:I_ow,:variable_str] = "I(O;W)"
    perf_res_MI[perf_res_MI[:variable].==:I_ao,:variable_str] = "I(A;O)"
    perf_res_MI[perf_res_MI[:variable].==:I_awgo,:variable_str] = "I(A;W|O)"
    perf_res_MI[perf_res_MI[:variable].==:I_aw,:variable_str] = "I(A;W)"

    nrows = size(perf_res_utils,1)
    perf_res_utils[:variable_str] = ["" for x in 1:nrows]
    perf_res_utils[perf_res_utils[:variable].==:E_U,:variable_str] = "E[U]"
    perf_res_utils[perf_res_utils[:variable].==:Objective_value,:variable_str] = "Objective J"

    plot_convergence_MI = plot(perf_res_MI, x=:iteration, y="value", color="variable_str", Geom.line,
                               BAtheme(key_position = :bottom; theme_args...), Guide.ylabel("[bits]"),
                               Guide.xlabel(xlabel_perf), Guide.colorkey(""))

    plot_convergence_Util = plot(perf_res_utils, x=:iteration, y="value", color="variable_str", Geom.line,
                                 BAtheme(key_position = :bottom; theme_args...), Guide.ylabel("[utils]"),
                                 Guide.xlabel(xlabel_perf), Guide.colorkey(""), BAdiscretecolorscale(2))

    if suppress_vis == false
        plt_convergence = vstack(plot_convergence_MI, plot_convergence_Util) #stack plots vertically
        display(plt_convergence)
    end

    return plot_convergence_MI, plot_convergence_Util

end

#plots mutual informations, entropies and EU, J as (stacked) bars
function plot_three_var_performancemeasures(performance_df::DataFrame, max_utility, β1, β2, β3; theme_args...)
    #extract relevant fields in correct order from data frame
    bitval_all = [performance_df[end,:I_ow], performance_df[end,:I_ao], performance_df[end,:I_awgo], performance_df[end,:I_aw],
                  performance_df[end,:H_ogw], performance_df[end,:H_ago], performance_df[end,:H_agow], performance_df[end,:H_agw]]

    #the color values will specify the legend entries
    color_label_bitvals = ["I(O;W)", "I(A;O)", "I(A;W|O)", "I(A;W)",
                           "H(O|W)", "H(A|O)", "H(A|O,W)", "H(A|W)"]

    #the x-values are used for stacking the bars and denoting the corresponding entropic term
    x_label_bitvals = ["H(O)", "H(A)", "H(A|O)", "H(A) ",
                       "H(O)","H(A)","H(A|O)", "H(A) "]  #be careful, two of the H(A) strings have an additional space, to
                                                         #make sure that it corresponds to two different bars on the x-axis
                                                         #(though the stacked height of both bars should be equal)

    #manually set y-limits, because automatic setting fails in Gadfly when using stacked bars
    maxy = maximum(bitval_all[1:4]+bitval_all[5:end])

    #plot stacked bars for the entropic terms, that are a sum of a mutual information term and
    #a conditional term
    p_composed = plot(x=x_label_bitvals, y=bitval_all, color=color_label_bitvals, Geom.bar(position=:stack),
                      Guide.ylabel("[bits]"), Guide.xlabel(""), Guide.title("Composition of entropic terms"),
                      Guide.colorkey(""), Scale.y_continuous(minvalue=0, maxvalue=maxy), BAtheme(;theme_args...))

    subs = [1,2,3,4]
    p_MI = plot(x=color_label_bitvals[subs], y=bitval_all[subs], color=color_label_bitvals[subs], Geom.bar(),
                Guide.ylabel("[bits]"), Guide.xlabel(""), Guide.title("Mutual information terms"),
                Scale.y_continuous(minvalue=0), BAtheme(key_position = :none; theme_args...))

    #extract expected utility and value of objective and mutual information terms
    subs = [1,2,3]
    if(β3==0)
        #sequential case
        β3=eps()  #use tiny value for β3 to prevent NaN - I(A;W|O) should equal zero in the sequential case anyway
    end

    inv_betas = [1/β1, 1/β2, 1/β3]
    perf = [performance_df[end,:E_U]; bitval_all[subs] .* inv_betas; performance_df[end,:Objective_value]]

    #these labels will define the legend entries
    label_perf = ["E[U]", "1/β1 I(O;W)", "1/β2 I(A;O)", "1/β3 I(A;W|O)", "J = E[U] - 1/β ∑ I"]
    x_label_perf = ["E[U]", "E[U] ", "E[U] ", "E[U] ", "E[U] "]  #be careful, four of the entries have an additonal space,
                                                                 #to make sure that there are two stacked bars (that will
                                                                 #have the same height)

    cont_colors = [Scale.color_continuous().f(p) for p in linspace(0, 1, 2)]
    disc_colors = Scale.color_discrete_hue().f(3)
    all_colors = [cont_colors[1]; disc_colors; cont_colors[2]]
    colscale = Scale.color_discrete_manual(all_colors...)

    p_perf = plot(x=x_label_perf, y=perf, color=label_perf, Geom.bar(position=:stack),
                  Guide.ylabel("[utils]"), Guide.xlabel(""), Guide.colorkey(""),
                  Guide.title("Utility vs. processing-cost trade-off"),
                  colscale, Scale.y_continuous(minvalue=0, maxvalue=max_utility), BAtheme(;theme_args...))

    return p_MI, p_composed, p_perf
end


text:
using Printf
using Plots

n =200
del, E, F ,B = init(n)

x = collect(0:n+1)*pi/(n+1)
U0 = 2.0*sin.(x)
x =x[2:n+1]
U0 = U0[2:n+1]

Ut = zeros(n)
Ux= 2.0*cos.(x)

w1 = del
w2 = del.*(E.+1.0)
w0 = hcat(w1,w2)

function runSim(w0,U0,Ut,Ux,tEnd)

tStart = 0;
t1 =tStart;
deltaT = 0.005;
#tEnd = 1.5;
w = zeros(n,2)
sols = zeros(n,2)
dt =0.001
while t1 < tEnd

w, dt, j1 ,j2 = FVMIBL(w0,U0,Ut,Ux);
ww = w
del = ww[:,1]
E = (ww[:,2]./ww[:,1]) .- 1.0
ncell = length(del)
Csep = zeros(ncell-2)

Csep = abs.(diff(del[2:end]))./abs.(diff(del[1:end-1]));
#xtrunc = x[2:end-1]
sep = abs.(diff(Csep))
for i =1:length(j2)-1
   if j1[i] > 0.0002
        println("singularity (separation) detected at t=$t1, x=$(x[i]./pi), j1=$(j1[i])")
    end
end

# the plots of del and E

#display(plot(x/pi,[del, E], xticks = 0:0.1:1, layout=(2,1), legend = false))

# convergence of the Eigen-value
p1 = plot(x[1:end-1]/pi,j1)
p2 = plot(x[1:end-1]/pi,j2)

p3 = plot(x/pi,del)
p4 = plot(x/pi,E)
display(plot(p1, p2, p3, p4, xticks = 0:0.1:1, layout=(2,2), legend = false))

#display(plot(sep, xticks = 0:10:200, legend = false))
sleep(0.05)
w0 = w;

@printf("Time :%1.10f , Time step size %1.10f \n", t1, dt);
#map(x -> @sprintf("Seperation occure at :%1.10f  \n",x), xtrunc[seperation]./pi);

t1 = t1+dt

sols = w
end

return w
end

max_stars_repo_name: desanga/IBMFVM

max_issues_repo_name: desanga/IBMFVM

max_forks_repo_name: desanga/IBMFVM


text:
mutable struct CholeskyType{T}
    c::Vector{T} # diagonal entries.
    a::Vector{T} # offdiagonal entries.
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
module Terms

using ..is_call
using ..ContrastsMatrix

import StreamModels.name

using ArgCheck

abstract type Term end

struct Intercept <: Term end

struct Interaction{Ts} <: Term
    terms::Ts
end

struct Continuous <: Term
    name::Symbol
end

struct Categorical{N} <: Term
    name::Symbol
    contrasts::ContrastsMatrix
    invindex::Dict
end

struct Eval <: Term
    name::Symbol
end

# TODO: capture the original expression (and possibly the symbol for the original call)
struct FunctionTerm{F,Forig}  <: Term where {F <: Function, Forig <: Function}
    f::F
    forig::Forig
    name::Expr
end

struct FormulaTerm{L,R}
    lhs::L
    rhs::R
end

"""
    haslhs(ft::FormulaTerm{L})

Return true is the LHS type is not Nothing
"""
haslhs(ft::FormulaTerm{L}) where L = L!==Nothing

function name(t::Term) end
name(t::Union{Eval, Continuous, Categorical}) = t.name

Base.string(t::Terms.Term) = string(name(t))

"""
    ismatrixterm(T::Type{<:Term})
    ismatrixterm(t::T) where T<:Term

Return true if this type should be included in a single model matrix (e.g., a
random effect term should return false).
"""
ismatrixterm(::T) where T<:Term = false

for T in [FunctionTerm, Categorical, Continuous, Intercept, Interaction]
    @eval ismatrixterm(::S) where S<:$T = true
end

"""
    flatten_formula(term::FormulaTerm)

Extract a vector of "flattened" terms.  If present, the first element will be
the LHS.  The next will be a vector of all the "matrix terms", and any other
terms will be passed through unchanged

"""
function flatten_formula(term::FormulaTerm)
    # create a tuple of terms where we:
    # 1. pull out LHS is present
    # 2. collect LHS matrix terms
    # 3. pass the rest on alone
    new_terms = []
    matrix_terms = []
    for t in term.rhs
        if ismatrixterm(t)
            push!(matrix_terms, t)
        else
            push!(new_terms, t)
        end
    end
    !isempty(matrix_terms) && unshift!(new_terms, tuple(matrix_terms...))
    haslhs(term) && unshift!(new_terms, term.lhs)
    tuple(new_terms...)
end

"""
    term_ex_from_formula_ex(x)

As part of the `@formula` macro, take a part of the formula and generate an
expression to construct the appropriate `Term` type.

# Extending the formula language

When `x` is an `Expr(:call, head, args...)`, this function calls the method
`term_ex_from_formula_ex(Val(head), x)`.  This allows for the formula DSL to be extended
by adding methods for your own operators or functions.  They should return an
`Expr` that constructs a `Term` based on the arguments of the `Expr`.

# Capturing arbitrary functions

For call `Expr`s that don't have specialized methods for `Val{head}`, a
`FunctionTerm` is generated, which wraps an anonymous function that takes a
single named tuple argument and calls the original function, replacing symbols
in the original `Expr` with fields of the named tuple.
"""
term_ex_from_formula_ex(::Nothing) = nothing
term_ex_from_formula_ex(i::Integer) = (@argcheck(i==1); :(Terms.Intercept()))
term_ex_from_formula_ex(s::Symbol) = Expr(:call, :(Terms.Eval), Meta.quot(s))

# calls dispatch on Val{head} for extensibility:
function term_ex_from_formula_ex(ex::Expr)
    @argcheck is_call(ex)
    term_ex_from_formula_ex(Val(ex.args[1]), ex)
end

term_ex_from_formula_ex(::Val{:~}, ex::Expr) =
    Expr(:call, :(Terms.FormulaTerm),
         term_ex_from_formula_ex(ex.args[2]),
         term_ex_from_formula_ex(ex.args[3]))
term_ex_from_formula_ex(::Val{:+}, ex::Expr) =
    Expr(:vect, [term_ex_from_formula_ex(x) for x in ex.args[2:end]]...)
term_ex_from_formula_ex(::Val{:&}, ex::Expr) =
    Expr(:call, :(Terms.Interaction),
         Expr(:tuple, [term_ex_from_formula_ex(x) for x in ex.args[2:end]]...))

# generic: capture calls as anonymous functions of named tuple.  so take
# something like log(1+a) and convert to (tup) -> log(1+tup[:a])
function term_ex_from_formula_ex(::Val{<:Any}, ex::Expr)
    tup_sym = gensym()
    anon_expr = Expr(:(->), tup_sym, replace_symbols!(copy(ex), tup_sym))
    f_orig = ex.args[1]
    Expr(:call, :(Terms.FunctionTerm), anon_expr, f_orig, Meta.quot(ex))
end

replace_symbols!(x, tup::Symbol) = x
replace_symbols!(x::Symbol, tup::Symbol) = Expr(:ref, tup, Meta.quot(x))
function replace_symbols!(ex::Expr, tup::Symbol)
    if is_call(ex)
        ex.args[2:end] .= [replace_symbols!(x, tup) for x in ex.args[2:end]]
    end
    ex
end

end # module Terms


text:
###############################################################################
# Extract designs from Sort-seq run
# -----------------------------------------------------------------------------
# Author:  Thomas E. Gorochowski <thomas.gorochowski@bristol.ac.uk>
# Licence: MIT
###############################################################################

using CSV
using DataFrames
using Query
using Glob
using Printf
using BioSequences
using BioAlignments
using FASTX

###############################################################################
# PARAMETERS – Should be optimised for particular use case
###############################################################################

# Minimum alignment score for start and end of memory register
# Previously 300, 100
global MIN_ALN_SCORE_REGION = 20
global MIN_ALN_SCORE_SPACER = 0

# Alignemnt parameters (these have been tweaked for 25 bp barcodes)
#global ALN_MATCH = 5
#global ALN_MISMATCH = -4
#global ALN_GAP_OPEN = -4
#global ALN_GAP_EXTEND = -2

# Parameters from Matt's paper (need to adjust the thresholds)
global ALN_MATCH = 2
global ALN_MISMATCH = -3
global ALN_GAP_OPEN = -5
global ALN_GAP_EXTEND = -2

global BASE_FRACTION_CALL_INSERT = 0.5

###############################################################################
# FUNCTIONS TO DO THE HEAVY LIFTING
###############################################################################

function load_references(filename)
    ref_seqs = Dict{String,LongDNASeq}()
    reader = FASTA.Reader(open(filename, "r"))
    record = FASTA.Record()
    read_count = 0
    while !eof(reader)
        read!(reader, record)
        ref_seqs[FASTA.identifier(record)] = FASTA.sequence(record)
        read_count = read_count + 1
    end
    close(reader)
    @info string("Processed ", filename, " containing ", read_count, " sequences")
    return ref_seqs
end

refs = load_references("data/references.fasta")

function ref_alignment_range(aln)
    start_pos = -1
    end_pos = -1
    aln_anchors = aln.a.aln.anchors
    first_match = true
    for a in aln_anchors
        if a.op == OP_SEQ_MATCH
            if first_match == true
                start_pos = a.refpos
                first_match = false
                end_pos = a.refpos
            else
                end_pos = a.refpos
            end
        end
    end
    return (start_pos, end_pos)
end

function check_length_fix_orientation(test_seq, start_seq, end_seq)
    # Alignment type used and scoring model
    problem = SemiGlobalAlignment()
    scoremodel = AffineGapScoreModel(
                  match=ALN_MATCH,
                  mismatch=ALN_MISMATCH,
                  gap_open=ALN_GAP_OPEN,
                  gap_extend=ALN_GAP_EXTEND
                )
    test_seq_rc = reverse_complement(test_seq)
    alignment_start = pairalign(problem, start_seq, test_seq, scoremodel)
    alignment_start_rc = pairalign(problem, start_seq, test_seq_rc, scoremodel)
    alignment_end = pairalign(problem, end_seq, test_seq, scoremodel)
    alignment_end_rc = pairalign(problem, end_seq, test_seq_rc, scoremodel)
    score_start = score(alignment_start)
    score_start_rc = score(alignment_start_rc)
    score_end = score(alignment_end)
    score_end_rc = score(alignment_end_rc)
    # Test to see what start and end are most likely and return sequence
    if score_start > score_start_rc && score_end > score_end_rc &&
        score_start > MIN_ALN_SCORE_REGION && score_end > MIN_ALN_SCORE_REGION
        (s_idx0, s_idx1) = ref_alignment_range(alignment(alignment_start))
        (e_idx0, e_idx1) = ref_alignment_range(alignment(alignment_end))
        return test_seq[s_idx0:e_idx1]
    elseif score_start_rc > score_start && score_end_rc > score_end &&
        score_start_rc > MIN_ALN_SCORE_REGION && score_end_rc > MIN_ALN_SCORE_REGION
        (s_idx0, s_idx1) = ref_alignment_range(alignment(alignment_start_rc))
        (e_idx0, e_idx1) = ref_alignment_range(alignment(alignment_end_rc))
        return test_seq_rc[s_idx0:e_idx1]
    else
        return nothing
    end
end

function find_spacer_list(seq, ref_seqs)
    spacers = []
    # Alignment type used and scoring model
    problem = SemiGlobalAlignment()
    scoremodel = AffineGapScoreModel(
                  match=ALN_MATCH,
                  mismatch=ALN_MISMATCH,
                  gap_open=ALN_GAP_OPEN,
                  gap_extend=ALN_GAP_EXTEND
                )
    seq_len = length(seq)
    seq_rc = reverse_complement(seq)
    for (ref_key, ref_seq) in ref_seqs
        # Check reference is for a spacer
        if ref_key != "start" && ref_key != "end"
            spacer_aln = pairalign(problem, ref_seq, seq, scoremodel)
            spacer_aln_rc = pairalign(problem, ref_seq, seq_rc, scoremodel)
            spacer_score = score(spacer_aln)
            spacer_score_rc = score(spacer_aln_rc)
            # Check alignments meet minimum score
            if spacer_score > spacer_score_rc && spacer_score > MIN_ALN_SCORE_SPACER
                (s_idx0, s_idx1) = ref_alignment_range(alignment(spacer_aln))
                push!(spacers, [ref_key, s_idx0])
            elseif spacer_score_rc > spacer_score && spacer_score_rc > MIN_ALN_SCORE_SPACER
                (s_idx0, s_idx1) = ref_alignment_range(alignment(spacer_aln_rc))
                push!(spacers, [string(ref_key, "_r"), seq_len-s_idx1])
            end
        end
    end
    # Sort the found spacers on their position
    sort!(spacers, lt=(x,y)->isless(x[2], y[2]))
    # Generate string for configuration

    #=
    # Only return results where all spacers are present
    if length(spacers) == length(ref_seqs) - 2
        spacer_string = ""
        for el in spacers
            if spacer_string == ""
                spacer_string = el[1]
            else
                spacer_string = string(spacer_string, " ", el[1])
            end
        end
        if spacers == []
            return "None"
        else
            return spacer_string
        end
    else
        return "None"
    end
    =#

    spacer_string = ""
    for el in spacers
        if spacer_string == ""
            spacer_string = el[1]
        else
            spacer_string = string(spacer_string, " ", el[1])
        end
    end
    if spacers == []
        return "None"
    else
        return spacer_string
    end
end

function find_spacers(fastq_filename, ref_seqs)
    @info string("Started processing ", fastq_filename, "...\n")
    results = Dict{String, Int}()
    try
        reader = FASTQ.Reader(open(fastq_filename, "r"))
        record = FASTQ.Record()
        read_count = 0
        processed_read_count = 0
        while !eof(reader)
            # Comment out if processing full data
            if read_count >= 10000
                break
            end
            if mod(read_count, 1000) == 0
                print(".")
            end
            read!(reader, record)
            fixed_seq = check_length_fix_orientation(FASTQ.sequence(record), ref_seqs["start"], ref_seqs["end"])
            if fixed_seq != nothing
                # Search for order and orientation of spacers and add to list
                cur_key = find_spacer_list(fixed_seq, ref_seqs)
                if haskey(results, cur_key)
                    results[cur_key] += 1
                else
                    results[cur_key] = 1
                end
                processed_read_count = processed_read_count + 1
            end
            read_count = read_count + 1
        end
        close(reader)
        @info string("Processed ", fastq_filename, " containing ", read_count, " reads (",
                     processed_read_count, " full length memory registers)\n")
    catch y
        @info string("Error in ", fastq_filename, "\n")
        println(y)
    end
    return results
end

function write_results_to_file(results, filename)
    open(filename, "w") do file
        write(file, "spacer_config,count\n")
        for (key, val) in results
            write(file, string(key, ",", val, "\n"))
        end
    end
end

###############################################################################
# RUN THE ANALYSIS
###############################################################################

# Load the metadata for the run
metadata = DataFrame(CSV.File("metadata.csv", header=1, delim=","))

# Create mapping from barcode to sample name using the metadata for run 1
bc_map1 = @from i in metadata begin
    @where i.run == 1
    @select i.barcode=>i.name
    @collect Dict
end

# Generate the input and output file names
read_filenames = []
output_filenames = []

for (bc, name) in bc_map1
    read_filename = "fastq/barcode$(@sprintf("%02d", bc)).fastq"
    output_filename = "output/$(name).csv"
    push!(read_filenames, read_filename)
    push!(output_filenames, output_filename)
end

# Load the spacer sequences to search for
spacer_refs_filename = "spacers.fasta"
ref_seqs = load_spacers(spacer_refs_filename)

# Extract the spacer configurations from each sample
Threads.@threads for idx in 1:length(read_filenames)
    read_filename = read_filenames[idx]
    output_filename = output_filenames[idx]
    results = find_spacers(read_filename, ref_seqs)
    write_results_to_file(results, output_filename)
end
=#

max_stars_repo_name: BiocomputeLab/tx-valve-tools

max_issues_repo_name: BiocomputeLab/tx-valve-tools

max_forks_repo_name: BiocomputeLab/tx-valve-tools