text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################
"""
    assign_start_to_variables!(data)
This function gives the measurement value in the data dictionary
as a starting value to its associated variable.
"""
function assign_start_to_variables!(data::Dict{String, Any})
    for (_,meas) in data["meas"]
        if all(i -> i ∈ [_DST.Normal{Float64}, Float64], typeof.(meas["dst"]))
            data[string(meas["cmp"])][string(meas["cmp_id"])]["$(string(meas["var"]))_start"] = _DST.mean.(meas["dst"])
        end
    end
end
"""
    assign_start_to_variables!(data, start_values_source)
This function assigns start values to the problem variables based on a dictionary where they are collected: start_values_source.
This dictionary must have the form of a powerflow solution dictionary.
"""
function assign_start_to_variables!(data::Dict{String, Any}, start_values_source::Dict{String, Any})
    for (_,meas) in data["meas"]
        msr_cmp = string(meas["cmp"])
        cmp_id = string(meas["cmp_id"])
        msr_var = string(meas["var"])
        if haskey(start_values_source["solution"][msr_cmp][cmp_id], msr_var)
            if msr_var != "w"
                data[msr_cmp][cmp_id]["$(msr_var)_start"] = start_values_source["solution"][msr_cmp][cmp_id][msr_var]
            else
                data[msr_cmp][cmp_id]["$(msr_var)_start"] = start_values_source["solution"][msr_cmp][cmp_id][msr_var][1]
            end
        else
            @warn "$(msr_var) is not in $(start_values_source["solution"][msr_cmp][cmp_id]), possible formulation mismatch"
        end
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_price_sensitive_load!(
    model::JuMP.Model,
    ps::PriceSensitiveLoad,
)::Nothing
    loads = _init(model, :loads)
    net_injection = _init(model, :expr_net_injection)
    for t in 1:model[:instance].time
        # Decision variable
        loads[ps.name, t] =
            @variable(model, lower_bound = 0, upper_bound = ps.demand[t])

        # Objective function terms
        add_to_expression!(model[:obj], loads[ps.name, t], -ps.revenue[t])

        # Net injection
        add_to_expression!(
            net_injection[ps.bus.name, t],
            loads[ps.name, t],
            -1.0,
        )
    end
    return
end

ext: jl


text:
__precompile__(true)

module StableMeander

  using Distributions
  using StatsBase
  using LambertW
  using FastGaussQuadrature
  using SpecialFunctions
  using DataFrames
  using DataFramesMeta

  export
    # distribution types
    Stable,
    PositiveStable,
    StableMeander,
    PreciseStableMeander,
    LocStableMeander,
    MvStableMeander,
    PreciseMvStableMeander,
    LocMvStableMeander,
    ConvexMinorantStableMeander,
    CrossConvexMinorantStableMeander,
    ConvexMinorantWeaklyStable,
    CrossConvexMinorantWeaklyStable,

    # methods
    rand,
    maximum,
    minimum,
    insupport,
    pdf,
    cdf,
    mgf,
    cf,
    mellin,
    mean,
    params,
    length

  # Source Files

  include("stable.jl")
  include("positivestable.jl")
  include("fdstablemeander.jl")
  include("convexminorantstable.jl")

end


text:
type NewtonSolver <: Solver1D
  solverInfo::SolverInfo
end
NewtonSolver(maxEvals::Int = 100, lowerBoundEnforced::Bool = false, upperBoundEnforced::Bool = false, lowerBound::Float64 = 0.0, upperBound::Float64 = 0.0) =
  NewtonSolver(SolverInfo(maxEvals, lowerBoundEnforced, upperBoundEnforced, lowerBound, upperBound))

# Newton Solver (safe)
function _solve(solver::NewtonSolver, f::Function, accuracy::Float64, xMin::Float64, xMax::Float64, fxMin::Float64, fxMax::Float64, root::Float64, eval_num::Int)
  # Orient the search so that f(xl) < 0
  if fxMin < 0.0
    xl = xMin
    xh = xMax
  else
    xh = xMin
    xl = xMax
  end

  # the stepsize before last
  dxold = xMax - xMin

  # and the last step
  dx = dxold

  froot = f(root)
  dfroot = f(root, Derivative()) # derivative, hurray multiple dispatch!
  eval_num += 1

  while eval_num < solver.solverInfo.maxEvals
    # bisect if out of range or not decreasing fast enough
    if ((root - xh) * dfroot - froot) * ((root - xl) * dfroot - froot) > 0.0 || abs(2.0 * froot) > abs(dxold * dfroot)
      dxold = dx
      dx = (xh - xl) / 2.0
      root = xl + dx
    else
      dxold = dx
      dx = froot / dfroot
      root -= dx
    end

    # Convergence criterion
    if abs(dx) < accuracy
      # root = f(root) # check this
      eval_num += 1
      return root
    end

    froot = f(root)
    dfroot = f(root, Derivative())
    eval_num += 1
    if froot < 0.0
      xl = root
    else
      xh = root
    end
  end

  error("Maximum number of function evals exceeded!")
end


text:
# This example highlights the use of a linear minimization oracle
# using an LP solver defined in MathOptInterface
# we compare the performance of the two LMOs, in- and out of place
#
# to get accurate timings it is important to run twice so that the compile time of Julia for the first run
# is not tainting the results

include("activate.jl")

using JSON
using LaTeXStrings
results = JSON.Parser.parsefile("movielens_result.json")

ref_optimum = results["trajectory_arr_lazy_ref"][end][2]

iteration_list = [
    [x[1] + 1 for x in results["trajectory_arr_fw"]],
    [x[1] + 1 for x in results["trajectory_arr_lazy"]],
    collect(1:1:length(results["function_values_gd"])),
]
time_list = [
    [x[5] for x in results["trajectory_arr_fw"]],
    [x[5] for x in results["trajectory_arr_lazy"]],
    results["timing_values_gd"],
]
primal_gap_list = [
    [x[2] - ref_optimum for x in results["trajectory_arr_fw"]],
    [x[2] - ref_optimum for x in results["trajectory_arr_lazy"]],
    [x - ref_optimum for x in results["function_values_gd"]],
]
test_list =
    [results["fw_test_values"], results["lazy_test_values"], results["function_values_test_gd"]]

label = [L"\textrm{FW}", L"\textrm{L-CG}", L"\textrm{GD}"]

FrankWolfe.plot_results(
    [primal_gap_list, primal_gap_list, test_list, test_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    [L"\textrm{Iteration}", L"\textrm{Time}", L"\textrm{Iteration}", L"\textrm{Time}"],
    [
        L"\textrm{Primal Gap}",
        L"\textrm{Primal Gap}",
        L"\textrm{Test Error}",
        L"\textrm{Test Error}",
    ],
    xscalelog=[:log, :identity, :log, :identity],
    legend_position=[:bottomleft, nothing, nothing, nothing],
    filename="movielens_result.pdf",
)


text:
fun mkcounter() {
    var i = 0;
    fun count() {
        i = i + 1;
        print i;
    }
    return count;
}

var counter = mkcounter();
counter();
counter();


text:
"""
    fixed_knockoffs(X::Matrix{T}; method=:sdp)

Creates fixed knockoffs based on equation (2.2)-(2.4) of
"Controlling the false discovery rate via Knockoffs" by Barber and Candes (2015)

# Inputs
+ `X`: A `n × p` numeric matrix, each row is a sample, and each column is normalized
to mean 0 variance 1 with unit norm.
+ `method`: :equi for equi-distant knockoffs (eq 2.3) or :sdp for SDP
knockoffs (eq 2.4)

# Output
+ `Knockoff`: A struct containing the original `X` and its knockoff `X̃`,
in addition to other variables (e.g. `s`)
"""
function fixed_knockoffs(X::Matrix{T}, method::Symbol) where T <: AbstractFloat
    n, p = size(X)
    n ≥ 2p || error("fixed_knockoffs: only works for n ≥ 2p case! sorry!")
    # compute gram matrix using full svd
    U, σ, V = svd(X, full=true)
    Σ = V * Diagonal(σ)^2 * V'
    Σinv = V * inv(Diagonal(σ)^2) * V'
    # compute s vector using the specified method
    if method == :equi
        λmin = typemax(T)
        for σi in σ
            σi^2 < λmin && (λmin = σi^2)
        end
        s = min(1, 2λmin) .* ones(size(Σ, 1))
    elseif method == :sdp
        svar = Variable(p)
        problem = maximize(sum(svar), svar ≥ 0, 1 ≥ svar, 2Σ - Diagonal(svar) in :SDP)
        solve!(problem, () -> SCS.Optimizer(verbose=false))
        s = clamp.(evaluate(svar), 0, 1) # for numeric stability
    else
        error("fixed_knockoffs: method can only be :equi or :sdp")
    end
    # compute Ũ such that Ũ'X = 0
    Ũ = @view(U[:, p+1:2p])
    # compute C such that C'C = 2D - D*inv(Σ)*D via eigendecomposition (cholesky not stable)
    D = Diagonal(s)
    γ, P = eigen(2D - D*Σinv*D)
    clamp!(γ, 0, typemax(T)) # numerical stability
    C = Diagonal(sqrt.(γ)) * P
    # compute knockoffs
    X̃ = X * (I - Σinv*D) + Ũ * C
    return Knockoff(X, X̃, s, Σ, Σinv)
end


text:
function clear(obj::StorelessUnivariateStatistic)
    return jcall(obj, "clear", void, ())
end

function copy(obj::StorelessUnivariateStatistic)
    return jcall(obj, "copy", StorelessUnivariateStatistic, ())
end

function get_result(obj::StorelessUnivariateStatistic)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::StorelessUnivariateStatistic, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
"""

Helper infrastructure package to compile and sample models using Stan's `cmdstan`.
Not really intended to be called directly by a user.

"""
module StanBase

using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using Unicode, DelimitedFiles
using DataFrames, CSV, Parameters

using StanDump

include("common/cmdstan_home.jl")
include("common/common_definitions.jl")
include("common/handle_keywords.jl")
include("common/update_model_file.jl")
include("common/update_R_files.jl")
include("common/par.jl")

include("stansamples/stan_summary.jl")
include("stansamples/read_summary.jl")

"""
The directory which contains the cmdstan executables such as `bin/stanc` and
`bin/stansummary`.

# Extended help

Inferred from the environment variable `JULIA_CMDSTAN_HOME` or `ENV["JULIA_CMDSTAN_HOME"]`
when available.

If these are not available, use `set_cmdstan_home!` to set the value of CMDSTAN_HOME.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`

Executing `versioninfo()` will display the value of `JULIA_CMDSTAN_HOME` if defined.
"""
CMDSTAN_HOME=""

function __init__()
  global CMDSTAN_HOME = if haskey(ENV, "CMDSTAN")
    ENV["CMDSTAN"]
  elseif isdefined(Main, :JULIA_CMDSTAN_HOME)
    Main.JULIA_CMDSTAN_HOME
  elseif haskey(ENV, "JULIA_CMDSTAN_HOME")
    ENV["JULIA_CMDSTAN_HOME"]
  elseif haskey(ENV, "CMDSTAN_HOME")
    ENV["CMDSTAN_HOME"]
  else
    @warn("Environment variable CMDSTAN_HOME not set. Use set_cmdstan_home!.")
    ""
  end
end

export
  CMDSTAN_HOME,
  get_cmdstan_home,
  set_cmdstan_home,
  make_command,
  set_make_string,
  CmdStanModels,
  StanModelError,
  read_summary,
  stan_summary

end # module


text:
export Market, compute_mu, compute_delta, I, J, K2

struct Market{T<:AbstractFloat}
    # TODO: ultimately a market is going to be a view of the larger dataset so these
    # owning types should be changed to parameterised types.
    x2::Matrix{T}      # J x K2
    shares::Vector{T}  # J
    weights::Vector{T} # I
    nodes::Matrix{T}   # I x K2

    # TODO: check correct dimensions
end

Base.eltype(::Type{Market{T}}) where {T} = T

""" Individuals in market. """
num_individuals(market::Market) = length(market.weights)

""" Products in market. """
num_products(market::Market) = size(market.x2, 1)

""" Demand-side nonlinear product characteristics. """
num_nonlinear_demand_characteristics(market::Market) = size(market.x2, 2)

""" "Compute the mean utility that solves the simple logit model. """
function compute_logit_delta(market)
    # log1p(-x) ≡ log(1-x) but faster and more accurately
    return log.(market.shares) .- log1p(-sum(market.shares))
end

""" Compute the mean utility for this market that equates observed and predicted market shares. """
function compute_delta(market::Market, iteration, sigma)
    result, _, _ = compute_delta_impl(market, iteration, sigma)
    return result
end

function compute_mu(market::Market, sigma)
    # Currently just support unobserved heterogeneity.
    return market.x2 * sigma * market.nodes'
end

# This code is used for both solving the contraction mapping as-is and for
# rules for automatic differentiation.
function compute_delta_impl(market::Market, iteration, sigma)
    mu = compute_mu(market, sigma)
    log_market_shares = log.(market.shares)

    I = num_individuals(market)
    J = num_products(market)

    # Speed is imperative for the contraction mapping so avoid unnecessary memory allocations by
    # pre-allocating and then using in-place operations. Unfortunately, this makes the code a bit
    # harder to understand.
    utilities = Matrix{eltype(market)}(undef, J, I)
    probabilities = Matrix{eltype(market)}(undef, J, I)
    shares = Vector{eltype(market)}(undef, J)

    function contraction!(delta_out::AbstractVector, delta_in::AbstractVector)
        @. utilities = delta_in + mu
        choice_probabilities!(probabilities, utilities)
        mul!(shares, probabilities, market.weights)  # Compute market share by integrating over individuals
        @. delta_out = delta_in + log_market_shares - log(shares)
    end

    return fixed_point_iteration(iteration, compute_logit_delta(market), contraction!), shares, probabilities
end


text:
module GaussianWeightedMeanPrecisionTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, isProper, unsafeMean, unsafeMode, unsafeVar, unsafeCov, unsafeMeanCov, unsafePrecision, unsafeMeanPrecision, unsafeWeightedMean, unsafeWeightedMeanPrecision
using ForneyLab: SPGaussianWeightedMeanPrecisionOutNPP, VBGaussianWeightedMeanPrecisionOut

@testset "dims" begin
    @test dims(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0)) == ()
    @test dims(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=ones(2), w=diageye(2))) == (2,)
end

@testset "vague" begin
    @test vague(GaussianWeightedMeanPrecision) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=tiny)
    @test vague(GaussianWeightedMeanPrecision, (2,)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(2), w=tiny*eye(2))
end

@testset "isProper" begin
    # Univariate
    @test isProper(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0))
    @test !isProper(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=-1.0))

    # Multivariate
    @test isProper(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[0.0], w=mat(1.0)))
    @test isProper(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=ones(2), w=diageye(2)))
    @test !isProper(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[0.0], w=mat(-1.0)))
end

@testset "==" begin
    # Univariate
    @test ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0)
    @test ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0) == ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)

    # Multivariate
    @test ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[0.0], w=mat(1.0)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[0.0], w=mat(1.0))
    @test ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[0.0], w=mat(1.0)) == ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0))
end

@testset "unsafe statistics" begin
    # Univariate
    @test unsafeMean(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 0.5
    @test unsafeMode(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 0.5
    @test unsafeVar(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 0.25
    @test unsafeCov(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 0.25
    @test unsafeMeanCov(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == (0.5, 0.25)
    @test unsafePrecision(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 4.0
    @test unsafeMeanPrecision(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == (0.5, 4.0)
    @test unsafeWeightedMean(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == 2.0
    @test unsafeWeightedMeanPrecision(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)) == (2.0, 4.0)

    # Multivariate
    @test unsafeMean(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == [0.5]
    @test unsafeMode(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == [0.5]
    @test unsafeVar(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == [0.25]
    @test unsafeCov(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == mat(0.25)
    @test unsafeMeanCov(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == ([0.5], mat(0.25))
    @test unsafePrecision(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == mat(4.0)
    @test unsafeMeanPrecision(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == ([0.5], mat(4.0))
    @test unsafeWeightedMean(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == [2.0]
    @test unsafeWeightedMeanPrecision(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))) == ([2.0], mat(4.0))
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=2.0), 1.0), -0.5723649429247)
    @test isapprox(logPdf(ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0, 2.0], w=[2.0 0.0; 0.0 2.0]), [1.0, 0.0]), -2.1447298858494)
end

@testset "convert" begin
    @test convert(ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}, ProbabilityDistribution(Univariate, GaussianMeanPrecision, m=0.5, w=4.0)) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)
    @test convert(ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.5, v=0.25)) == ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=2.0, w=4.0)
    @test convert(ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}, ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=[0.5], w=mat(4.0))) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))
    @test convert(ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.5], v=mat(0.25))) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[2.0], w=mat(4.0))
    @test convert(ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}, ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=2.0)) == ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(2.0))
end

#-------------
# Update rules
#-------------

@testset "SPGaussianWeightedMeanPrecisionOutNPP" begin
    @test SPGaussianWeightedMeanPrecisionOutNPP <: SumProductRule{GaussianWeightedMeanPrecision}
    @test outboundType(SPGaussianWeightedMeanPrecisionOutNPP) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(SPGaussianWeightedMeanPrecisionOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])
    @test !isApplicable(SPGaussianWeightedMeanPrecisionOutNPP, [Message{PointMass}, Nothing, Message{PointMass}])

    @test ruleSPGaussianWeightedMeanPrecisionOutNPP(nothing, Message(Univariate, PointMass, m=1.0), Message(Univariate, PointMass, m=2.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=2.0)
    @test ruleSPGaussianWeightedMeanPrecisionOutNPP(nothing, Message(Multivariate, PointMass, m=[1.0]), Message(MatrixVariate, PointMass, m=mat(2.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(2.0))
end

@testset "VBGaussianWeightedMeanPrecisionOut" begin
    @test VBGaussianWeightedMeanPrecisionOut <: NaiveVariationalRule{GaussianWeightedMeanPrecision}
    @test outboundType(VBGaussianWeightedMeanPrecisionOut) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(VBGaussianWeightedMeanPrecisionOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianWeightedMeanPrecisionOut, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBGaussianWeightedMeanPrecisionOut(nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=1.0, v=2.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=1.0, w=3.0)
    @test ruleVBGaussianWeightedMeanPrecisionOut(nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, PointMass, m=mat(3.0))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[1.0], w=mat(3.0))
end

end # module


text:
# This file calculates the roots of a sincfun in barycentric form
# based on the work of P.W. Lawrence, SIAM J. MATRIX ANAL. APPL., 34:1277--1300, 2013.

## TODO: This is simply the implementation via Julia's eig() routine. It could be faster
## using the specific algorithm in the paper. Also, once a differentiation command is
## implemented, it will be possible to find roots in BigFloat precision using Newton iteration.

export roots

function complexroots{D<:Domain,T<:Float64}(sf::sincfun{D,T})
    wv,xv,fv = (-one(T)).^(collect(0:sf.n-1)).*sf.ωv, ψ(sf.domain,convert(T,π)/2*sinh.(sf.jh)), sf.fϕv.*sf.ϕpv
    cutoff = abs.(wv) .≥ 10eps(T)
    wv,xv,fv = wv[cutoff],xv[cutoff],fv[cutoff]
    n = length(wv)
    wv,fv = wv/norm(wv),fv/norm(fv)
    A = sparse([zero(T) -fv'; wv diagm(xv)])
    s = [one(T);T[abs(fv[i]) == zero(T) ? one(T) : sqrt(abs(wv[i])/abs(fv[i])) for i=1:n]]
    S = sparse(collect(1:n+1),collect(1:n+1),s)
    B = diagm([zero(T);ones(T,n)])
    rts,V = eig(full(S\A*S),B)
    rts = rts[abs.(rts).<Inf]
end

function roots{D<:Domain,T<:Float64}(sf::sincfun{D,T})
    rts = complexroots(sf)
    rts = sort(real(rts[imag(rts) .== zero(T)]))
end

function roots{D<:Finite,T<:Float64}(sf::sincfun{D,T})
    rts = complexroots(sf)
    rts = real(rts[imag(rts) .== zero(T)])
    a,b = sf.domain.ab
    rts = sort(rts[a + 10sf.ωscale*eps(T) .< abs.(rts) .< b - 10sf.ωscale*eps(T)])
end

# This is a first attempt to use Lawrence's O(n^2) algorithm to transform the arrowhead matrix
# to upper Hessenberg.

#=
function complexroots{D<:Domain,T<:Float64}(sf::sincfun{D,T})
    wv,xv,fv = (-one(T)).^([0:sf.n-1]).*sf.ωv, sf.jh, sf.fϕpv
    cutoff = abs.(wv) .≥ 10eps(T)
    wv,xv,fv = wv[cutoff],xv[cutoff],fv[cutoff]
    #wvint,xvint,fvint = interlace2(wv),interlace2(xv),interlace2(fv)
    t,d,c = toupperhessenberg(wv,xv,fv)
    A = diagm([zero(T),d])+diagm(t,1)+diagm(t,-1)
    [A[1,i] += c[i] for i=1:length(c)]
    B = diagm([zero(T),ones(T,length(xv))])
    rts,V = eig(A,B)
    rts = rts[abs.(rts).<Inf]
    rts = ψ.(sf.domain,convert(T,π)/2*sinh(rts))
end

function toupperhessenberg{T<:Number}(w::Vector{T},x::Vector{T},f::Vector{T})
    @assert length(w) == length(x) == length(f)
    n = length(w)
    w,f = w/norm(w),f/norm(f)
    t,d = zeros(T,n),zeros(T,n)
    q = zeros(T,n,n)
    q[:,1] = w
    t[1] = norm(q[:,1])
    q[:,1] /= t[1]
    d[1] = dot(q[:,1],x.*q[:,1])
    q[:,2] = (x-d[1]).*q[:,1]
    t[2] = norm(q[:,2])
    q[:,2] /= t[2]
    for i=2:n-1
        d[i] = dot(q[:,i],x.*q[:,i])
        q[:,i+1] = (x-d[i]).*q[:,i]-t[i]*q[:,i-1]
        t[i+1] = norm(q[:,i+1])
        q[:,i+1] /= t[i+1]
    end
    d[n] = dot(q[:,n],x.*q[:,n])
    c = [zero(T), vec(-f'*q)]
    c[2] -= t[1]
    return t,d,c
end
=#
# Base.LinAlg.givensAlgorithm(f::T,g::T) returns cs, sn, r


text:
struct Point{T}
    x::T
    y::T
end
function convert(::Type{Point},x::Vector,y::Vector) where T
    @assert length(x) == length(y)
    [Point(x[ix],y[ix]) for ix in 1:length(x)]
end
eltype(p::Point) = eltype(p.x)

"get x coordinates from a vector of `Point`"
getx(v::Vector{Point{T}})  where T = T[v[i].x for i in 1:length(v)]

"get y coordinates from a vector of `Point`"
gety(v::Vector{Point{T}})  where T = T[v[i].y for i in 1:length(v)]

"get x-y coordinates from a vector of `Point`"
coords(v::Vector{Point{T}}) where T = [getx(v) gety(v)]

function show(io::IO, ::MIME"text/plain", p::Point{T}) where T
    print(io,"Point of type $T:\n")
    print(io,"      x = $(p.x)\n")
    print(io,"      y = $(p.y)\n")
end

show(io::IO,p::Point{T}) where T = print(io,"($(p.x),$(p.y))")

# arithmetics
(+)(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y)
(-)(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y)
(*)(n::Number, p::Point) = Point(n*p.x, n*p.y)
(*)(p::Point, n::Number) = n*p
(/)(p::Point, n::Number) = Point(p.x/n, p.y/n)

# comparison
(==)(p1::Point, p2::Point) = (p1.x == p2.x) && (p1.y == p2.y)
(==)(x::T, p2::Point{T}) where T = (x == p2.x)  # is x-coord part of this point?
isapprox(p1::Point, p2::Point) = isapprox(p1.x, p2.x) && isapprox(p1.y, p2.y)
isless(p1::Point,p2::Point) = p1.x < p2.x   # caution we sort on x only here!
in(x::T, p2::Vector{Point{T}}) where T = any(x .== p2)

# Promotion
Base.zero(::Type{Point{T}}) where {T} = Point(zero(T),zero(T))
Base.promote_rule(::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{Point{T1}}, ::Type{T2}) where {T1,T2<:Number} = Point{promote_type(T1,T2)}
Base.promote_op(::typeof(*), ::Type{T1}, ::Type{Point{T2}}) where {T1<:Number,T2} = Point{promote_type(T1,T2)}

"""
# MLine

A `MLine` is a vector of `Point`[@ref]. The x-coordinates of the points can be irregularly spaced.

## Fields

* `v`: Vector of `Point`s
* `iextrap`: indices of points that have been extrapolated
"""
mutable struct MLine{T} <: AbstractArray{T,1}
    v::Vector{Point{T}}
    iextrap::Vector{Int}
    function MLine(v::Vector{Point{T}}; extrapolated=nothing) where {T<:Number}
        this = new{T}()
        this.v = v
        if isnothing(extrapolated)
            this.iextrap = Int[]
        else
            this.iextrap = extrapolated
        end
        return this
    end
    function MLine(x::Vector{T},y::Vector{T}; extrapolated=nothing) where {T<:Number}
        this = new{T}()
        n = length(x)
        @assert n == length(y)
        this.v = [Point(x[i],y[i]) for i in 1:n]
        if isnothing(extrapolated)
            this.iextrap = Int[]
        else
            this.iextrap = extrapolated
        end
        return this
    end
end

# typemin
Base.eltype(l::MLine) = eltype(l.v)
function typemin(L::MLine{T}) where {T<:Number}
    typemin(eltype(L))
end

# indexing
Base.size(l::MLine) = (length(l.v),)
Base.length(l::MLine) = length(l.v)

function getindex(l::MLine,i::Int)
    # @boundscheck checkbounds(l,i)
    l.v[i]
end
Base.IndexStyle(::Type{<:MLine}) = IndexLinear()
getindex(l::MLine, I...) = l.v[I]
endof(l::MLine) = length(l.v)

# iteration
Base.iterate(L::MLine, state = 1) = state > length(L.v) ? nothing : (L[state], state + 1)

# printing
function show(io::IO, ::MIME"text/plain", L::MLine{T}) where {T<:Number}
    xvec = [i.x for i in L.v]
    print(io,"$T MLine\n")
    print(io,"number of points: $(length(L.v))\n")
    print(io,"xrange: $(round.(extrema(getx(L)),digits = 2))\n")
    ys= gety(L)
    print(io,"yrange: $(round.(extrema(ys[isfinite.(ys)]),digits = 2))\n")
end
show(io::IO,L::MLine{T}) where {T<:Number} = print(io,"$(length(L.v)) point $T MLine")

getv(l::MLine{T}) where T = l.v
getx(l::MLine{T}) where T = getx(l.v)
gety(l::MLine{T}) where T = gety(l.v)
getex(l::MLine{T}) where T = l.iextrap

coords(l::MLine{T}) where T = coords(l.v)
function unique!(l::MLine{T}) where T
    l.v = unique(l.v)
end

function floory!(l::MLine{T},yy::T) where {T<:Number}
    ix = findall(gety(l) .< yy)
    if length(ix) > 0
        newp = l.v[ix]
        for i in 1:length(newp)
            splice!(l.v,ix[i],[Point(newp[i].x,yy)])
        end
    end
end

# Array of MLine
function floory!(L::Array{MLine{T}},x::T) where {T<:Number}
    for i in eachindex(L)
        floory!(L[i],x)
    end
end
function gety(L::Vector{MLine{T}}) where {T<:Number}
    out = fill(zero(T),(length(L),length(L[1].v)))
    for i in eachindex(L)
        out[i,:] = gety(L[i])
    end
    out
end

function max_x(l::MLine{T}) where {T<:Number}
    lo = typemin(T)
    for i in l
        if i.x > lo
            lo = i.x
        end
    end
    return lo
end
function min_x(l::MLine{T}) where {T<:Number}
    hi = typemax(T)
    for i in l
        if i.x < hi
            hi = i.x
        end
    end
    return hi
end
function max_y(l::MLine{T}) where {T<:Number}
    lo = typemin(T)
    for i in l
        if i.y > lo
            lo = i.y
        end
    end
    return lo
end
function findmax_y(l::MLine{T}) where {T<:Number}
    lo = typemin(T)
    ix = 0
    imax = 0
    for i in l
        ix += 1
        if i.y > lo
            lo = i.y
            imax = ix
        end
    end
    return (lo,imax)
end
function min_y(l::MLine{T}) where {T<:Number}
    hi = typemax(T)
    for i in l
        if i.y < hi
            hi = i.y
        end
    end
    return hi
end

# interpolating a MLine

"""
   interp(l::MLine{T},ix::Vector{T};Bool::extrap = false) where {T<:Number}

Interpolate a `MLine` on a vector of values `x`.
Importantly, this returns a new vector of `Point` (i.e. tuples of (x,y), and not just a function value y).
"""
function interp(l::MLine{T},ix::Vector{T};extrap::Bool = true) where {T<:Number}
    # # whenever
    # xex = extrema(ix)
    # @debug(logger,"interpolating $ix ")
    xvec = getx(l)
    if !issorted(xvec)
        # println(xvec)
    end
    xrange = extrema(xvec)

    fi = findall((ix .< xrange[1]) .| (ix .> xrange[2]))

    # by default extrapolate all lines
    # but record which points have been extrapolated
    if length(fi) > 0
        if extrap
            itp = extrapolate(interpolate((xvec,),l.v,Gridded(Linear())),Line())
            out = MLine(itp(ix),extrapolated = fi)
        else
            out = MLine(zeros(T,length(ix)),zeros(T,length(ix)))
            for xi in fi
                out.v[xi] = Point(ix[xi],typemin(T))
            end
            for xi in setdiff(1:length(ix),fi)
            # for xi in collect((1:length(xvec)))[Not(fi)]
                out.v[xi] = interpolate((xvec,),l.v,Gridded(Linear()))(ix[xi])
            end
        end
    else
        itp = interpolate((xvec,),l.v,Gridded(Linear()))
        out = MLine(itp(ix))
    end

    return out
end

function interp(e::Array{MLine{T}},ix::Vector{T};extrap::Bool = true) where {T<:Number}
    [interp(e[i],ix,extrap=extrap) for i in eachindex(e)]
end

"""
    linemax(e::Array{MLine{T}}) where {T<:Number}

For an array of `MLine`s on identical support `xx`,
computes the index in `e` of the `MLine` where the `y`-value is highest for each `xx`.
One can imagine `e` as a matrix where each row represents the `y`-values from a
different `MLine`; this function returns the index of the column-wise maximum.

"""
function linemax(e::Array{MLine{T}}) where {T<:Number}
    out = zeros(Int,length(e[1]))
    rows = length(e)  # how many lines
    for j in eachindex(e[1])   # across columns
        ic = 0  # current index
        iv = 0  # index of best value
        v = typemin(T)
        for row in 1:rows
            ic += 1
            # if (e[row][j].y > v) && ( !(j ∈ e[row].iextrap) )   # if best value and not extrapolated
            if (e[row][j].y > v)    # if best value
                v = e[row][j].y
                iv = ic
            end
        end
        out[j] = iv
    end
    return out
end
# appending, prepending , deleting and splitting at

"prepend `Point`s to a `MLine`[@ref]"
function prepend!(m::MLine{T},p::Vector{Point{T}}) where T
    prepend!(m.v,p)
    # reconfigure!(m)
end

"delete an index"
function delete!(m::MLine,idx)
    deleteat!(m.v,idx)
    # reconfigure!(m)
end

"append points to a MLine"
function append!(m::MLine{T},p::Vector{Point{T}}) where T
    append!(m.v,p)
    # reconfigure!(m)
end

"insert a single value at an interior index"
function insert!(m::MLine{T},v::Point{T},idx::Int) where T
    insert!(m.v,idx,v)
    # reconfigure!(m)
end

"""
    splitat(m::MLine,j::Int,repeat_boundary::Bool=true)

Splits a `MLine` object after given index and returns 2 new `MLine`s as a tuple.
If `repeat_boundary` is true, then the separating index is the first point of the second new `MLine`.
Propagates the `iextrap` vector of indices to the newly formed `MLine`s.
"""
function splitat(m::MLine,j::Int,repeat_boundary::Bool=true)
    iex = getex(m)
    if length(iex) > 0
        # m had any extrapolated points
        # keep track of those in the new object
        m1 = MLine(m.v[1:j],extrapolated = any(iex .< j+1) ? iex[iex .< j+1] : nothing)

        if repeat_boundary
            m2 = MLine(m.v[j:end],extrapolated = any(iex .>= j) ? iex[iex .>= j] .- (j-1) : nothing)

        else
            m2 = MLine(m.v[j+1:end],extrapolated = any(iex .> j) ? iex[iex .> j] .- j : nothing)
        end
    else
        m1 = MLine(m.v[1:j])
        if repeat_boundary
            m2 = MLine(m.v[j:end])
        else
            m2 = MLine(m.v[j+1:end])
        end
    end
    return (m1,m2)
end

"sort a `MLine` along x-grid"
function sortx!(m::MLine)
    sort!(m.v)  # sorts v
    # reconfigure!(m)
end


text:
# FASTA Index
# ===========
#
# Index for random access to FASTA files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# http://www.htslib.org/doc/faidx.html
struct Index
    # offset for the record's sequence by header: See above specification
    names::Dict{String, Int}
    lengths::Vector{Int}
    offsets::Vector{Int}
    linebases::Vector{Int}
    linewidths::Vector{Int}
end

function Index(filepath::AbstractString)
    return open(read_faidx, filepath)
end

function read_faidx(input::IO)
    names = Dict{String, Int}()
    lengths = Int[]
    offsets = Int[]
    linebases = Int[]
    linewidths = Int[]
    for (i, line) in enumerate(eachline(input))
        values = split(chomp(line), '\t')
        name = values[1]
        length = parse(Int, values[2])
        offset = parse(Int, values[3])
        linebase = parse(Int, values[4])
        linewidth = parse(Int, values[5])
        names[name] = i
        push!(lengths, length)
        push!(offsets, offset)
        push!(linebases, linebase)
        push!(linewidths, linewidth)
    end
    return Index(names, lengths, offsets, linebases, linewidths)
end

function Base.getindex(index::Index, name::AbstractString)
    i = get(index.names, convert(String, name), nothing)
    if i === nothing
        throw(ArgumentError("sequence \"$(name)\" is not in the index"))
    elseif i == 1
        offset = 0
    # Else, we go to the previous one and calculate the length of the previous
    # sequence in bytes, then seek to right after that one.
    else
        prev_offset = index.offsets[i - 1]
        prev_len = index.lengths[i - 1]
        prev_linebase = index.linebases[i - 1]
        prev_linewidth = index.linewidths[i - 1]

        newline_len = prev_linewidth - prev_linebase
        len = cld(prev_len, prev_linebase) * newline_len + prev_len
        offset = prev_offset + len
    end
    return i
end

# Set the reading position of `input` to the starting position of the record `name`.
# TODO: Return the seek position to be consistent with Base?
function seekrecord(input::IO, index::Index, name::AbstractString)
    i = index[name]
    if i == 1
        offset = 0
    # Else, we go to the previous one and calculate the length of the previous
    # sequence in bytes, then seek to right after that one.
    else
        prev_offset = index.offsets[i - 1]
        prev_len = index.lengths[i - 1]
        prev_linebase = index.linebases[i - 1]
        prev_linewidth = index.linewidths[i - 1]

        # Note: newline_len may differ between sequences in the same file, as per
        # the specification.
        newline_len = prev_linewidth - prev_linebase
        len = cld(prev_len, prev_linebase) * newline_len + prev_len
        offset = prev_offset + len
    end
    seek(input, offset)
    return nothing
end

max_forks_repo_name: BioJulia/FASTX


text:
export reversed

using ReversedStrings
import ReversedStrings: reversed

reversed(x::Tuple) = reverse(x)

function set_capture(sequence::ReversedString, index::Int, x)
    @warn "check"
    set_capture(sequence.x,index,x)
end

## caveat!
regex_string(x::ReversedString) = regex_escape(x.representation)

export PositiveLookbehind
"""
    PositiveLookbehind(parser)

Parser that succeeds if and only if `parser` succeeds **before cursor**. Consumes no input.
The match is returned.
Useful for checks like "must be preceded by `parser`, don't consume its match".
"""
@auto_hash_equals struct PositiveLookbehind{S,T,P} <: WrappedAssertion{S,T}
    parser::P
    function PositiveLookbehind(p_,reversed_parser=true)
        p = reversed_parser ? deepmap_parser(reversed,IdDict(),parser(p_)) : parser(p_)
        new{Tuple{Int,state_type(p)},result_type(p),typeof(p)}(p)
    end
end
# result_type(p::Type{PositiveLookbehind{T}}) where T = T
regex_prefix(x::PositiveLookbehind) = "(?<="

export NegativeLookbehind
"""
    NegativeLookbehind(parser)

Parser that succeeds if and only if `parser` does not succeed **before cursor**.  Consumes no input.
`nothing` is returned as match.
Useful for checks like "must not be preceded by `parser`, don't consume its match".

```jldoctest
julia> la=NegativeLookbehind("keep")
re"(?<!keep)"

julia> parse("peek"*la,"peek")
("peek", re"(?<!keep)")
```
"""
@auto_hash_equals struct NegativeLookbehind{P} <: WrappedAssertion{MatchState,NegativeLookbehind{P}}
    parser::P
    function NegativeLookbehind(p_,reversed_parser=true)
        p = reversed_parser ? deepmap_parser(reversed,IdDict(),parser(p_)) : parser(p_)
        new{typeof(p)}(p)
    end
end
regex_prefix(x::NegativeLookbehind) = "(?<!"

export Lookbehind
"""
    Lookbehind(does_match::Bool, p)

[`PositiveLookbehind`](@ref) if `does_match==true`,
[`NegativeLookbehind`](@ref) otherwise.
"""
function Lookbehind(does_match::Bool, p)
    if does_match
        PositiveLookbehind(p)
    else
        NegativeLookbehind(p)
    end
end
@deprecate look_behind(does_match,p) Lookbehind(does_match, p)

function _iterate(t::NegativeLookbehind, str, till, posi, next_i, state::Nothing)
    rseq=reversed(str)
    next_i < 1 && return next_i, MatchState()
    r = _iterate(t.parser, rseq, till,
                 reverse_index(rseq,_prevind(str,next_i)), nothing)
    if r === nothing
        next_i,MatchState()
    else
        nothing
    end
end

_iterate(t::PositiveLookbehind, str, till, posi, next_i, state::MatchState) =
    nothing

function _iterate(t::PositiveLookbehind, str, till, posi, next_i, state)
    rseq=reversed(str)
    ri = reverse_index(rseq,_prevind(str,next_i))
    next_i < 1 && return nothing
    r = _iterate(t.parser, rseq, till, ri, tuple_pos(state,ri), tuple_state(state))
    if r === nothing
        nothing
    else
        next_i, r
    end
end

"""
    Base.get(parser::PositiveLookbehind, sequence, till, after, i, state)

get result of PositiveLookbehind

!!! note
    The result is currently for a `reversed` sequence, and
    you might find it difficult to [`map`](@ref) a lookbehind parser match.
    If you require this functionality please open an issue for discussion.

    Assertions do not consume input, so typically these input chars are parsed/mapped outside of the assertion.

```jldoctest
julia> p = Sequence(!re"a+b", PositiveLookbehind(!re"a+b"))
🗄 Sequence
├─ 🗄 Sequence |> !
│  ├─ a+  |> Repeat
│  └─ b
└─ (?<=🗄) PositiveLookbehind
   ├─ b
   └─ a+  |> Repeat
::Tuple{SubString,SubString}

julia> p("aaab")
("aaab", "baaa")
```
"""
function Base.get(parser::PositiveLookbehind, sequence, till, after, i, state)
    rseq = reversed(sequence)
    after_ = tuple_pos(state)
    get(parser.parser, rseq, till, after_, reverse_index(rseq,prevind(sequence, i)), tuple_state(state))
end

regex_inner(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    regex_inner(reversed(x.parser))

children(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    reverse(children(x.parser))

for T in [PositiveLookahead,NegativeLookahead,PositiveLookbehind,NegativeLookbehind]
    eval(quote
         deepmap_parser(f::Function,mem::AbstractDict,x_::$T,a...; kw...) =
         let x = deepmap_parser(f,mem,x_.parser,a...; kw...)
         $T(x,false)
         end
         end)
end
reversed(x::CombinedParser) = deepmap_parser(reversed, x)
deepmap_parser(f::Function,mem::AbstractDict,x::Union{AtStart,AtEnd,Always,Never},a...; kw...) = f(x)
reversed(x::Union{Always,Never}) = x
reversed(x::AtStart) = AtEnd()
reversed(x::AtEnd) =  AtStart()
reversed(x::NIndexParser{0}) = x
reversed(x::NIndexParser{1}) = x
reversed(x::ConstantParser) = x
reversed(x::ConstantParser{<:AbstractString}) =
    ConstantParser(reversed(x.parser))

deepmap_parser(::typeof(reversed),mem::AbstractDict,x::Sequence) =
    get!(mem,x) do
        Sequence(( deepmap_parser(reversed,mem,p) for p in reverse(x.parts) )...)
    end
deepmap_parser(::typeof(reversed),mem::AbstractDict,x::Atomic) =
    get!(mem,x) do
        Atomic(deepmap_parser(reversed,mem,x.parser))
    end
deepmap_parser(::typeof(reversed),mem::AbstractDict,x::NegativeLookbehind) =
    get!(mem,x) do
        NegativeLookahead(x.parser) ##deepmap_parser(reversed,@show x.parser))
    end
deepmap_parser(::typeof(reversed),mem::AbstractDict,x::NegativeLookahead) =
    get!(mem,x) do
        NegativeLookbehind(x.parser) ##deepmap_parser(reversed,x.parser))
    end
deepmap_parser(::typeof(reversed),mem::AbstractDict,x::PositiveLookbehind) =
    get!(mem,x) do
        PositiveLookahead(x.parser) ##deepmap_parser(reversed,x.parser))
    end
deepmap_parser(::typeof(reversed),mem::AbstractDict,x::PositiveLookahead) =
    get!(mem,x) do
        PositiveLookbehind(x.parser) ##deepmap_parser(reversed,x.parser))
    end


text:
## ==================================================================================================
## Hack to be removed or updated
## FIXME ME ON FIRE use the MB.jl prescribed interface method instead:
##   https://juliamanifolds.github.io/Manifolds.jl/stable/examples/manifold.html#manifold-tutorial
## ==================================================================================================

import ApproxManifoldProducts: coords, uncoords, getPointsManifold, _makeVectorManifold

export SE2E2_Manifold

# this is a hack and not fully implemented as per:
struct _SE2E2 <: MB.AbstractManifold{MB.ℝ} end

const SE2E2_Manifold = _SE2E2()

MB.manifold_dimension(::_SE2E2) = 5

AMP.coords(::Type{<:typeof(SE2E2_Manifold)}, p::Manifolds.ProductRepr) = [p.parts[1][1], p.parts[1][2], atan(p.parts[2][2,1],p.parts[2][1,1]), p.parts[3][1], p.parts[3][2]]

function AMP.uncoords(::typeof(SE2E2_Manifold), p::AbstractVector{<:Real})
  α = p[3]
  Manifolds.ProductRepr(([p[1], p[2]]), [cos(α) -sin(α); sin(α) cos(α)], ([p[4], p[5]]))
end

function AMP.getPointsManifold(mkd::ManifoldKernelDensity{M}) where {M <: typeof(SE2E2_Manifold)}
  data_ = getPoints(mkd.belief)
  [uncoords(mkd.manifold, view(data_, :, i)) for i in 1:size(data_,2)]
end

function Statistics.mean(::typeof(SE2E2_Manifold), pts::AbstractVector)
  se2_ = (d->Manifolds.ProductRepr(d.parts[1], d.parts[2])).(pts)
  mse2 = mean(Manifolds.SpecialEuclidean(2), se2_)
  e2_ = (d->Manifolds.ProductRepr(d.parts[3])).(pts)
  me2 = mean(Euclidean(2), e2_)
  Manifolds.ProductRepr(mse2.parts[1], mse2.parts[2], me2.parts[1])
end

AMP._makeVectorManifold(::M, prr::ProductRepr) where {M <: typeof(SE2E2_Manifold)} = coords(M, prr)

## =============================================================================
## Needs a home
## =============================================================================

export BearingRange_Manifold

struct _CircleEuclid <: MB.AbstractManifold{MB.ℝ} end

MB.manifold_dimension(::_CircleEuclid) = 2

const BearingRange_Manifold = _CircleEuclid()

AMP.coords(::Type{<:typeof(BearingRange_Manifold)}, p::ProductRepr) = [p.parts[1][1]; p.parts[2][1]]

function AMP.uncoords(::typeof(BearingRange_Manifold), p::AbstractVector{<:Real})
  ProductRepr(([p[1];]), ([p[2];]))
end

function AMP.getPointsManifold(mkd::ManifoldKernelDensity{M}) where {M <: typeof(BearingRange_Manifold)}
  data_ = getPoints(mkd.belief)
  [uncoords(mkd.manifold, view(data_, :, i)) for i in 1:size(data_,2)]
end

function Statistics.mean(::typeof(BearingRange_Manifold), pts::AbstractVector)
  TensorCast.@cast bearing[i] := pts[i][1]
  TensorCast.@cast range_[i] := pts[i][2]
  mc = mean(Circle(), bearing)
  mr = mean(range_)

  return [mc; mr]
end

AMP._makeVectorManifold(::M, prr::ProductRepr) where {M <: typeof(BearingRange_Manifold)} = coords(M, prr)

## ==================================================================================================


text:
@rule typeof(dot)(:in2, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass, meta::AbstractCorrection) = begin
    A  = mean(m_in1)
    out_wmean, out_prec = weightedmean_precision(m_out)

    ξ = A * out_wmean
    W = correction!(meta, A * out_prec * A')

    return convert(promote_variate_type(variate_form(m_in1), NormalWeightedMeanPrecision), ξ, W)
end


text:
# fallback methods

(A::LinearOperator)(b::Sequence) = *(A, b)

function Base.:*(A::LinearOperator, b::Sequence)
    _iscompatible(domain(A), space(b)) || return throw(DimensionMismatch)
    codomain_A = codomain(A)
    CoefType = promote_type(eltype(A), eltype(b))
    c = Sequence(codomain_A, Vector{CoefType}(undef, dimension(codomain_A)))
    _mul!(c, A, b, true, false)
    return c
end

@inline function LinearAlgebra.mul!(c::Sequence, A::LinearOperator, b::Sequence, α::Number, β::Number)
    _iscompatible(space(c), codomain(A)) & _iscompatible(domain(A), space(b)) || return throw(DimensionMismatch)
    _mul!(c, A, b, α, β)
    return c
end

@inline function _mul!(c::Sequence, A::LinearOperator, b::Sequence, α::Number, β::Number)
    domain_A, codomain_A = domain(A), codomain(A)
    space_b = space(b)
    space_c = space(c)
    if domain_A == space_b
        if codomain_A == space_c
            __mul!(coefficients(c), coefficients(A), coefficients(b), α, β)
        else
            if iszero(β)
                coefficients(c) .= zero(eltype(c))
            elseif !isone(β)
                coefficients(c) .*= β
            end
            inds_space = indices(codomain_A ∩ space_c)
            @inbounds __mul!(view(c, inds_space), view(A, inds_space, :), coefficients(b), α, true)
        end
    else
        inds_mult = indices(domain_A ∩ space_b)
        if codomain_A == space_c
            @inbounds __mul!(coefficients(c), view(A, :, inds_mult), view(b, inds_mult), α, β)
        else
            if iszero(β)
                coefficients(c) .= zero(eltype(c))
            elseif !isone(β)
                coefficients(c) .*= β
            end
            inds_space = indices(codomain_A ∩ space_c)
            @inbounds __mul!(view(c, inds_space), view(A, inds_space, inds_mult), view(b, inds_mult), α, true)
        end
    end
    return c
end

function Base.:\(A::LinearOperator, b::Sequence)
    _iscompatible(codomain(A), space(b)) || return throw(DimensionMismatch)
    return Sequence(domain(A), \(coefficients(A), coefficients(b)))
end

function LinearAlgebra.ldiv!(A::LinearOperator, b::Sequence)
    _iscompatible(codomain(A), space(b)) || return throw(DimensionMismatch)
    return Sequence(domain(A), ldiv!(coefficients(A), coefficients(b)))
end

function LinearAlgebra.ldiv!(c::Sequence, A::LinearOperator, b::Sequence)
    _iscompatible(space(c), domain(A)) & _iscompatible(codomain(A), space(b)) || return throw(DimensionMismatch)
    ldiv!(coefficients(c), coefficients(A), coefficients(b))
    return c
end

# Cartesian spaces

@inline function _mul!(c::Sequence{<:CartesianSpace}, A::LinearOperator{<:CartesianSpace,<:CartesianSpace}, b::Sequence{<:CartesianSpace}, α::Number, β::Number)
    domain_A, codomain_A = domain(A), codomain(A)
    space_b = space(b)
    space_c = space(c)
    if domain_A == space_b
        if codomain_A == space_c
            __mul!(coefficients(c), coefficients(A), coefficients(b), α, β)
        else
            if iszero(β)
                coefficients(c) .= zero(eltype(c))
            elseif !isone(β)
                coefficients(c) .*= β
            end
            inds_space = indices(codomain_A ∩ space_c)
            @inbounds __mul!(view(c, inds_space), view(A, inds_space, :), coefficients(b), α, true)
        end
    else
        m = nb_cartesian_product(domain_A)
        n = nb_cartesian_product(codomain_A)
        if iszero(β)
            coefficients(c) .= zero(eltype(c))
        elseif !isone(β)
            coefficients(c) .*= β
        end
        @inbounds for j ∈ 1:m
            bⱼ = component(b, j)
            @inbounds for i ∈ 1:n
                _mul!(component(c, i), component(A, i, j), bⱼ, α, true)
            end
        end
    end
    return c
end

@inline function _mul!(c::Sequence{<:VectorSpace}, A::LinearOperator{<:CartesianSpace,<:VectorSpace}, b::Sequence{<:CartesianSpace}, α::Number, β::Number)
    domain_A, codomain_A = domain(A), codomain(A)
    space_b = space(b)
    space_c = space(c)
    if domain_A == space_b
        if codomain_A == space_c
            __mul!(coefficients(c), coefficients(A), coefficients(b), α, β)
        else
            if iszero(β)
                coefficients(c) .= zero(eltype(c))
            elseif !isone(β)
                coefficients(c) .*= β
            end
            inds_space = indices(codomain_A ∩ space_c)
            @inbounds __mul!(view(c, inds_space), view(A, inds_space, :), coefficients(b), α, true)
        end
    else
        if iszero(β)
            coefficients(c) .= zero(eltype(c))
        elseif !isone(β)
            coefficients(c) .*= β
        end
        @inbounds for j ∈ 1:nb_cartesian_product(domain_A)
            _mul!(c, component(A, j), component(b, j), α, true)
        end
    end
    return c
end

@inline function _mul!(c::Sequence{<:CartesianSpace}, A::LinearOperator{<:VectorSpace,<:CartesianSpace}, b::Sequence{<:VectorSpace}, α::Number, β::Number)
    domain_A, codomain_A = domain(A), codomain(A)
    space_b = space(b)
    space_c = space(c)
    if domain_A == space_b
        if codomain_A == space_c
            __mul!(coefficients(c), coefficients(A), coefficients(b), α, β)
        else
            if iszero(β)
                coefficients(c) .= zero(eltype(c))
            elseif !isone(β)
                coefficients(c) .*= β
            end
            inds_space = indices(codomain_A ∩ space_c)
            @inbounds __mul!(view(c, inds_space), view(A, inds_space, :), coefficients(b), α, true)
        end
    else
        @inbounds for i ∈ 1:nb_cartesian_product(codomain_A)
            _mul!(component(c, i), component(A, i), b, α, β)
        end
    end
    return c
end


text:
using Plots

"""
    runge_kutta(f::Function, x₀::Real, y₀::Real, x::Real, h::Real=10^(-3))::Array{Float64,1}
Return an array with values of function which is solution to the differential equation
given by function f with intial condition
"""
function runge_kutta(f::Function, x₀::Real, y₀::Real, x::Real, h::Real=10^(-3))::Array{Float64,1}
    n = round(Int64, x / h)
    array = zeros(n)
    for i in 1:n
        k₁ = f(x₀, y₀)
        k₂ = f(x₀ + 1/2*h, y₀ + 1/2*h*k₁)
        k₃ = f(x₀ + 1/2*h, y₀ + 1/2*h*k₂)
        k₄ = f(x₀ + h, y₀ + h*k₃)
        y₀ = y₀ + h / 6 * (k₁ + 2*k₂ + 2*k₃ + k₄)
        x₀ = x₀ + h
        array[i] = y₀
    end
    return array
end

# logistic growth for chosen constants
f(t, y) = 0.6y*(1 - y / 10)
array = runge_kutta(f, 0, 1, 20)
𝐫 = plot(LinRange(0, 20, length(array)), array)
display(𝐫)

# https://en.smath.com/wiki/GetFile.aspx?File=Examples/RK4SystemEquations.pdf
# examples below are taken from this url

array  = [(t, y₁, y₂) -> sin(t) + cos(y₁) + sin(y₂),
        (t, y₁, y₂) -> cos(t) + sin(y₂)]
X = 0
Y = [-1., 1.]

"""
    evaluate(array::Array{Function, 1}, X::Real, Y::Array{Real, 1})::Array{Real, 1}
Return an array with values of function f(t, y₁, y₂) defined in array for given point X
"""
function evaluate(array::Array{Function, 1}, X::Real, Y::Array{T, 1})::Array{T, 1} where T <: Real
    return [i(X, Y[1], Y[2]) for i in array]
end

"""
    runge_kutta(F::Array{Function, 1}, x::Real, X::Real, Y::Array{T, 1}, h::Real=1e-3)::Array{T, 2} where T <: Real
Return an array with values of functions which are solutions to the system of differential equations defined by the array.
X and Y are initial conditions, x is the maximum point.
"""
function runge_kutta(F::Array{Function, 1}, x::Real, X::Real, Y::Array{T, 1}, h::Real=1e-3)::Array{T, 2} where T <: Real
    n = round(Int64, x / h)
    array = zeros((n, 2))
    for i in 1:n
        k₁ = evaluate(F, X, Y)
        k₂ = evaluate(F, X + 1/2 * h, Y .+ 1/2 * h * k₁)
        k₃ = evaluate(F, X + 1/2 * h, Y .+ 1/2 * h * k₂)
        k₄ = evaluate(F, X + h, Y .+ h * k₃)
        Y = @. Y + h / 6 * (k₁ + 2*k₂ + 2*k₃ + k₄)
        X = X + h
        array[i, :] = Y
    end
    return array
end

a = runge_kutta(array, 20, X, Y)

𝑝 = plot(LinRange(0, 5, length(a[:, 1])), a[:, 1])
plot!(LinRange(0, 5, length(a[:, 2])), a[:, 2])
display(𝑝)
# everything works fine

# harmonic oscillator ẍ + x = 0
# can be rewritten as
# ẋ = y, ẏ = -x
# note that y is velocity ⃗v

array2 = [(t, y₁, y₂) -> y₂, (t, y₁, y₂) -> -y₁]
Y2 = [0., 1.]
b = runge_kutta(array2, 30, 0, Y2)

𝑞 = plot(LinRange(0, 5, length(b[:, 1])), b[:, 1])
display(𝑞)
# plot!(LinRange(0, 5, length(b[:, 2])), b[:, 2])
𝑤 = plot(b[:, 1], b[:, 2])  # a kind of phase plane plot
display(𝑤)

# vector field for harmonic oscillator
xs = LinRange(0, 5, length(b[1, :]))

xs = xs[1:150:end]
xs |>println
y₁s = b[:, 1][1:150:end]
y₂s = b[:, 2][1:150:end]

plot(quiver(y₁s, y₂s, quiver=(0.2 * array2[1].(xs, y₁s, y₂s), 0.2 * array2[2].(xs, y₁s, y₂s)), color=:red))


text: @testset "zygoterules" begin end


text:
immutable Spec1D{T<:Real,D<:LineDomains,
        G<:LineGroups,W<:RedundantWaveletClass} <: AbstractSpec{T,D,G,W}
    ɛ::Float64
    ϕmeta::ΦMeta
    ψmetas::Array{ΨMeta1D,3}
    class::W
    domain::D
    log2_size::Int
    max_qualityfactor::Float64
    max_scale::Float64
    motherfrequency::Float64
    n_filters_per_octave::Int
    n_octaves::Int
    pointgroup::G
    signaltype::Type{T}
    function (::Type{Spec1D}){T,D,G,W}( ;
            class::W = Morlet(),
            pointgroup::G = TrivialGroup(),
            signaltype::Type{T} = Float32,
            domain::D = FourierDomain(1),
            ɛ = default_ɛ(signaltype),
            log2_size = 15,
            max_qualityfactor = nothing,
            max_scale = Inf,
            n_filters_per_octave = nothing,
            n_octaves = nothing,
            tuningfrequency = nothing)
        """Maximum quality factor and number of filters per octave are
        set simultaneously, so that they default to each other if either
        is present."""
        max_qualityfactor, n_filters_per_octave =
             default_max_qualityfactor(max_qualityfactor, n_filters_per_octave),
             default_n_filters_per_octave(n_filters_per_octave, max_qualityfactor)
        """The center frequency of the mother wavelet (`motherfrequency`)
        may be tuned to user-defined frequency (e.g. 440 Hz in music) or left
        as a default."""
        motherfrequency =
            tune_motherfrequency(tuningfrequency, class, n_filters_per_octave)
        """By default, the number of octaves in the filter bank is such that
        all wavelet scales do not exceed the signal length, nor the
        user-defined maximum scale `max_scale` if present."""
        n_octaves = default_n_octaves(n_octaves, class, log2_size,
            max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave)
        """The number of orientations of a one-dimensional filter bank is
        equal to 1 for real input (`pointgroup::TrivialGroup`), equal to 2 for
        complex input (`pointgroup::ReflectionGroup`)."""
        nΘs = get_n_orientations(pointgroup)
        ħ = uncertainty(class)
        """The meta-information of the wavelets is computed as a 3-d array,
        whose dimensions respectively correspond to
        * the orientation variable `θ`,
        * the chroma variable `χ`, and
        * the octave variable `j`."""
        ψmetas = Array{ΨMeta1D}(nΘs, n_filters_per_octave, n_octaves)
        for j in 0:(n_octaves-1), χ in 0:(n_filters_per_octave-1)
            γ = j * n_filters_per_octave + χ
            resolution = exp2(-γ / n_filters_per_octave)
            centerfrequency = motherfrequency * resolution
            unbounded_scale = ħ * max_qualityfactor / centerfrequency
            scale = min(unbounded_scale, max_scale)
            unbounded_q = scale * centerfrequency / ħ
            qualityfactor = clamp(unbounded_q, 1.0, max_qualityfactor)
            bandwidth = centerfrequency / qualityfactor
            for θ in 0:(nΘs-1)
                ψmetas[1+θ, 1+χ, 1+j] = ΨMeta1D(γ, θ, χ, bandwidth,
                    centerfrequency, j, qualityfactor, scale)
            end
        end
        """The bandwidth of the lowpass filter `ϕ` is such that the Fourier
        transform of `ϕ` will cover the low-frequency part of the spectrum
        that is not covered by the wavelets `ψs`."""
        ϕbandwidth = motherfrequency * exp2(-n_octaves)
        ϕscale = ħ / ϕbandwidth
        ϕmeta = ΦMeta(ϕbandwidth, ϕscale)
        spec = new{T,D,G,W}(ɛ, ϕmeta, ψmetas, class, domain, log2_size,
            max_qualityfactor, max_scale, motherfrequency, n_filters_per_octave,
            n_octaves, pointgroup, signaltype)
        """Before returning the `spec`, we call the function `checkspec`
        which enforces properties of the wavelet filter bank to satisfy null
        mean, limited spatial support, and Littlewood-Paley inequality."""
        checkspec(spec) && return spec
    end
end

checkspec(spec::Spec1D) = checkspec_super(spec) && return true

"""Given a dimensionless tuning frequency, returns the maximal admissible
mother frequency such that the subsequent wavelets will be in tune with the
tuning frequency.

For example, to tune a 12-chroma Morlet filter bank to a concert pitch of
440 Hz at a sample rate of 44,1 kHz:
    `Morlet1DSpec(n_filters_per_octave = 12, tuning_frequency = 440.0/44100.0)`
"""
function tune_motherfrequency(tuningfrequency, spectype, n_filters_per_octave)
    max_centerfrequency =
        default_motherfrequency(spectype, n_filters_per_octave)
    tuning_ratio = max_centerfrequency / tuningfrequency
    tuning_γ = floor(n_filters_per_octave*log2(tuning_ratio))
    return tuningfrequency * exp2(tuning_γ/n_filters_per_octave)
end
tune_motherfrequency(tuningfrequency::Void, spectype, n_filters_per_octave) =
    default_motherfrequency(spectype, n_filters_per_octave)


text:
module Model

using Distributions

using RCall  # Mclust

import LinearAlgebra  # Identity matrix
import Random  # shuffle, seed!
import StatsBase  # wsample, counts

include("../MCMC/MCMC.jl")
import .MCMC.Util.@namedargs

include("util.jl")
include("DICparam.jl")
include("State.jl")
include("Data.jl")
include("Constants.jl")
include("Tuners.jl")
include("repFAM/repFAM.jl")
include("update.jl")
include("genInitialState.jl")
include("fit.jl")

#=
precompile(cytof_fit, (State, Constants, Data, Int, Int,
                       Vector{Vector{Symbol}}, Vector{Int}, Bool, Int, Bool,
                       Bool, Bool))
=#

end # Model

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
# This file contains extensions for base `println`
# which make it easy to print real values and arrays

"""
    println(io::IO, value::RT)

Print an [`RT`](@ref) value.

# Usage

```jldoctest
using Scats.Internal.Prec: RT
using Scats.Internal.Extras: println

x = RT(1.0)
println(x)

# output

 1.000000000000000E+00
```

Note: this doctest requires x64 architecture.
"""
function Base.println(io::IO, value::RT)
    println(io, sprintf1(RF, value))
end

"""
    println(io::IO, array::Array{RT, 1})

Print an array containing [`RT`](@ref) values.

# Usage

```jldoctest
using Scats.Internal.Prec: RT
using Scats.Internal.Extras: println

x = Array{RT}([1.0, 2.0, 3.0])
println(x)

# output

 1.000000000000000E+00    2.000000000000000E+00    3.000000000000000E+00
```

Note: this doctest requires x64 architecture.
"""
function Base.println(io::IO, array::Vector{RT})
    println(io, join([sprintf1(RF, s) for s in array], " "^3))
end


text:
import DataStructures
using DataStructures

"""
    compute_contr

Gives the contrast of an image
Example:
A = rand(Int16, 70, 80)
"""
function compute_contr(A)
    n,m = size(A)
    sz = 0.05 * (n * m)
    # println(sz)
    pq_max = PriorityQueue()
    pq_min = PriorityQueue()
    convert(Int,sz)
    cntr=0
    for i in eachindex(A)

            cntr+=1
            if cntr<=sz
                enqueue!(pq_max, i, A[i])
            else
                enqueue!(pq_max, i, A[i])
                dequeue!(pq_max)
            end
            if cntr<=sz
                enqueue!(pq_min, i, -A[i])
            else
                enqueue!(pq_min, i, -A[i])
                dequeue!(pq_min)
            end
        end
        contrast = 0
        bright = 0
        dark = 0
        cntr=sz
        while cntr>0

            cntr-=1
            contrast +=peek(pq_max).second
            bright +=peek(pq_max).second
            contrast +=peek(pq_min).second
            dark -=peek(pq_min).second
            dequeue!(pq_max)
            dequeue!(pq_min)
        end

        contrast /= (2.0*sz)
        bright /= sz
        dark /= sz

        return bright, dark
end


text:
struct TerminalTestMDP <: MDP{Union{Int,TerminalState}, Int} end
@test terminalstate == terminalstate
@test isterminal(TerminalTestMDP(), terminalstate)
@test promote_type(Int, TerminalState) == Union{Int, TerminalState}
@test typeof([1,terminalstate]) == Vector{statetype(TerminalTestMDP())}


text:
using ParetoFront
using Test

# Suppose we are in two dimensions and want min, min to dominate.
# We know then, that the point that has a minimum value of y will dominate
# all points to the right of it, ie  at x > x(at ymin) will be dominated
# and so for all x on the Pareto front we must have x < x(at ymin)
# and also the point with the minimum value of x will dominate all
# points that are above it, ie that y > y(at xmin) will be dominated
# and so for all y on the Pareto front we must have y < y(at xmin)
function test_min_blob(n)
    v = randn(n, 2)
    ((xmin, ymin), (ix, iy)) = findmin(v, dims=1)
    ymax = v[ix[1], 2]
    xmax = v[iy[1], 1]
    pareto = get_pareto(v, min_idxs=[1,2], max_idxs=[])
    front = Matrix(reduce(hcat, collect(pareto))')
    pxmax, pymax = maximum(front, dims=1)
    return pxmax <= xmax && pymax <= ymax
end

# the test for max, max as optimal
function test_max_blob(n)
    v = randn(n, 2)
    ((xmax, ymax), (ix, iy)) = findmax(v, dims=1)
    ymin = v[ix[1], 2]
    xmin = v[iy[1], 1]
    pareto = get_pareto(v, min_idxs=[], max_idxs=[1,2])
    front = Matrix(reduce(hcat, collect(pareto))')
    pxmin, pymin = minimum(front, dims=1)
    return pxmin >= xmin && pymin >= ymin
end

@testset "ParetoFront.jl" begin
    for i in 10 .^ (1:6)
        @test test_min_blob(i)
        @test test_max_blob(i)
    end
end


text:
function FWI{T<:AbstractFloat}(vp0::Array{T,2},d::Array{T,3},wav::Array{T,1},isz::Array{Int,1},isx::Array{Int,1},igz::Array{Int,1},igx::Array{Int,1},ot::Array{T,1},fmin::T,fmax::T,nf::Int,dz::T,dx::T,dt::T,ext::Int=50,atten_max::T=2.,GNiter=5,CGiter=20,CGtol=1.0e-100)

# This function performs frequency-domain acoustic full waveform inversion
# using the Gauss-Newton method. A multi-scale approach is adopted by
# fitting the lowest frequencies first. The data is assumed to consist of
# time-domain shot gathers that have been organized into a data cube such
# that the three dimensions are number of time steps, number of receivers,
# and number of shots, respectively. It is assumed that each shot uses the
# same array of receivers (receivers are not moved) and has the same duration
# for the recording period. Additionally, it is assumed that each shot has
# the same source wavelet. Full waveform inversion is performed for only a
# few selected frequencies falling between 'fmin' and 'fmax'; the frequencies
# are selected using the scheme proposed by Sirgue and Pratt, 2004.
#
# REFERENCES: 1) Laurent Sirgue and R. Gerhard Pratt (2004), Efficient waveform
#                inversion and imaging: A strategy for selecting temporal
#                frequencies.
#
# INPUTS:     vp0       - Initial guess for the velocity model
#             d         - Data (time domain shot gathers) orgaized into a cube with dimensions nt,ng,ns
#             wav       - Source wavelet (time domain)
#             isz       - z-indices for source locations
#             isx       - x-indices for source locations
#             igz       - z-indices for receiver locations
#             igx       - x-indices for receiver locations
#             ot        - Activation time for each of the sources
#             fmin      - Minimum frequency to model
#             fmax      - Maximum frequency to model
#             nf        - Number of frequency bins; 'nf' must be greater than or equal to 'length(wav)'
#             dz        - Grid spacing in z-direction
#             dx        - Grid spacing in x-direction
#             dt        - Sampling interval (time step)
#             ext       - Thickness (number of grid points) of absorbing boundary region
#             atten_max - Maximum complex amplitude in the attenuating boundary layer
#             GNiter    - Maximum number of iterations for the Gauss-Newton method
#             CGiter    - Maximum number of iterations for conjugate gradients
#             CGtol     - Convergence tolerance for conjugate gradients
#
# OUTPUTS:    vp        - Estimated velocity model
#             cost      - Normalized cost function at each iteration (for each frequency)
#             f_inv     - A list of the frequencies that FWI was performed for

    nz = size(vp0,1) + 2*ext
    nx = size(vp0,2) + 2*ext
    nt = size(d,1)
    ng = size(d,2)
    ns = size(d,3)

    L = Laplacian(nz,nx,dz,dx)
    M = MassMatrix(vp0,nz,nx,ext)
    A = Attenuation(nz,nx,ext,atten_max)
    R = Restriction(nz,nx,ext,igz,igx)
    param = Dict(:R=>R,:A=>A)

    D = fft([d ; zeros(nf-nt,ng,ns)],1)
    WAV = fft([wav ; zeros(nf-length(wav))])
    fs = 1/dt
    df = fs/nf
    faxis = fftshift(-fs/2:df:fs/2-df)
    waxis = 2*pi*faxis

    zmax = dz*(nz-2*ext)
    hmax = 0.5*maximum(dx*(kron(isx,ones(Int,ng)) - kron(ones(Int,ns),igx)))
    alpha = 1/sqrt(1 + (hmax/zmax)^2)

    _,iwmin = findmin(abs(fmin-faxis))
    f_inv = [faxis[iwmin]]
    w_inv = [waxis[iwmin]]
    iw_inv = [iwmin]
    while f_inv[end]/alpha <= fmax
        _,iw = findmin(abs(f_inv[end]/alpha-faxis))
        push!(f_inv,f_inv[end]/alpha)
        push!(w_inv,2*pi*f_inv[end])
        push!(iw_inv,iw)
    end

    cost = Array{Float64}[]
    for iw in iw_inv
        COST = Float64[]
        w = waxis[iw]
        s = zeros(eltype(D),nz*nx,ns)
        for ishot = 1:ns
            s[:,ishot] = Source(isz[ishot],isx[ishot],ot[ishot],WAV,waxis,w,nz,nx,ext)
        end
        for iter = 1:GNiter
            H = lufact(L + w^2*M*A)
            U = zeros(eltype(D),nz*nx,ns)
            r = zeros(eltype(D),ng*ns)
            for ishot = 1:ns
                U[:,ishot] = H\s[:,ishot]
                r[(ishot-1)*ng+1:ishot*ng] = R*U[:,ishot] - D[iw,:,ishot]
            end
            param[:H] = H
            param[:U] = U
            param[:w] = w
            dM,_ = ConjugateGradients(r,[SensitivityMultiShot],[param],Niter=CGiter,mu=0.,tol=CGtol)
            M -= spdiagm(real(dM))
            push!(COST,norm(r))
        end
        push!(cost,COST/COST[1])
    end

    vp = reshape(sqrt(1./diag(M)),nz,nx)[ext+1:end-ext,ext+1:end-ext]

    return vp,cost,f_inv

end


text:
## A factorization type for upper diagonal matrices. One of:
##
## * RankOne{T,S}
## * Z{T,S}, pencil
## * Identity{T,S}
## * UnitaryDiagonal
## * UpperTriangular{T,S}

"""

An R factorization encapsupulate the `R` in a QR factorization. For the companion matrix case, this is a sparse factorization in terms of rotators. Other scenarios are different sub-types.

"""
abstract type AbstractRFactorization{T, S} <: LinearAlgebra.Factorization{S} end

# Basic interface for an RFactorization:
#
# getindex (need just k-2 <= j <= k, but as general as possible)
# passthrough!(RF, U) to pass a rotator through from left to right or right to left
# simple_passthrough!(RF, U) for a shortcut passthrough
# we also might need:
# length
# eltype
# Matrix
# zero(RF{T, Rt}) to give 0
# one(RF{T, Rt}) to give 1
Base.eltype(RF::AbstractRFactorization{T, S}) where {T,S} = S
Base.zero(RF::AbstractRFactorization) = zero(eltype(RF))
Base.one(RF::AbstractRFactorization) = one(eltype(RF))

## passthrough: RF <- Us
function passthrough!(RF::AbstractRFactorization, Us::Union{AscendingChain, DescendingChain})
    for i in  eachindex(Us.x)
        Us.x[i] =  passthrough!(RF, Us.x[i])
    end
end

# passthrough: Us -> RF
function passthrough!(Us::Union{AscendingChain, DescendingChain}, RF::AbstractRFactorization)
    for i in length(Us.x):-1:1
         Us.x[i] =  passthrough!(Us.x[i], RF)
    end
end

function passthrough!(RF::AbstractRFactorization, B::TwistedChain)
    length(B) == 0 && return nothing
    for i in iterate_rl(B.pv)
        B.x[i] = passthrough!(RF, B.x[i])
    end
end

function passthrough!(B::TwistedChain, RF::AbstractRFactorization)
    length(B) == 0 && return nothing
    for i in iterate_lr(B.pv)
        B.x[i] = passthrough!(B.x[i], RF)
    end
end

##################################################

## Rank one decomposition for amrvw algorithm to find eigenvalues of companion matrix
"""

A companion matrix will have a QR decomposition wherein R is essentially an identy plus a rank one matrix

"""
struct RFactorizationRankOne{T,S, V} <: AbstractRFactorization{T, S}
  Ct::AscendingChain{T,S,V}
  B::DescendingChain{T,S,V}
  D::SparseDiagonal{S}
end

Base.copy(RF::RFactorizationRankOne) = RFactorizationRankOne(copy(RF.Ct), copy(RF.B), copy(RF.D))
Base.length(RF::RFactorizationRankOne) = length(RF.Ct)
Base.size(RF::RFactorizationRankOne) = (length(RF)+1, length(RF)+1)

## getindex
## helper function
function _w(B, j, k)

    ck, sk = vals(B[k])
    if j == k
        return -sk
    else
        cj, sj = vals(B[j])
        s = one(sj)
        for i in (j+1):(k-1)
            ci, si = vals(B[i])
            s *= si
        end
        return conj(cj) * s *  ck
    end
end

## We have W = C * R, where W invovles B and D, and C involves Ct; This
## solves for R terms, as in the paper
## To see the W pattern
## function rotm(a::T,b, i, N) where {T}
##      r = diagm(0 => ones(T, N))
##      r[i:i+1, i:i+1] = [a b; -b conj(a)]
##      r
## end
## @vars bgc bhc bic bjc bkc blc
## @vars bgs bhs bis bjs bks bls
## @vars dg dh di dj dk dl dm
## Bg = rotm(bgc, bgs, 1, 7)
## Bh = rotm(bhc, bhs, 2, 7)
## Bi = rotm(bic, bis, 3, 7)
## Bj = rotm(bjc, bjs, 4, 7)
## Bk =  rotm(bkc, bks, 5, 7)
## Bl = rotm(blc, bls, 6, 7)
## D = diagm(0  => [dg, dh, di, dj, dk,dl, dm])
## en = [0,0,0,0, 0, 1,0]
## Ws = Bg * Bh  * Bi * Bj *  Bk *  *Bl * D   *  en
## To see the C pattern:
## @vars bc_i bs_i bc_j bs_j bc_k bs_k
## @vars cc_g cs_g cc_h cs_h cc_i cs_i cc_j cs_j cc_k cs_k
## @vars d_i d_j d_k d_l
## @vars what wg wh wi wj wk wl
## u = rotm(cc_k, cs_k, 1,2) * [what, wl]
## kk = u[1](what => solve(u[2], what)[1]) |> simplify

## u = rotm(cc_k, cs_k, 2,3) *  rotm(cc_j, cs_j, 1,3) * [what, wk, wl]
## kk_1 = u[1](what => solve(u[3], what)[1]) |> simplify

## u = rotm(cc_k, cs_k,3, 4)*rotm(cc_j, cs_j, 2,4) *  rotm(cc_i, cs_i, 1,4) * [what, wj, wk, wl]
## kk_2 = u[1](what => solve(u[4], what)[1]) |> simplify

## u = rotm(cc_k, cs_k, 4,5) * rotm(cc_j, cs_j,3, 5)*rotm(cc_i, cs_i, 2,5) *  rotm(cc_h, cs_h, 1,5) * [what, wi, wj, wk, wl]
## kk_3 = u[1](what => solve(u[5], what)[1]) |> simplify

## u = rotm(cc_k, cs_k, 5,6) * rotm(cc_j, cs_j, 4,6) * rotm(cc_i, cs_i,3, 6)*rotm(cc_h, cs_h, 2,6) *  rotm(cc_g, cs_g, 1,6) * [what, wh, wi, wj, wk, wl]
## kk_4 = u[1](what => solve(u[6], what)[1]) |> simplify
function Base.getindex(RF::RFactorizationRankOne, j, k)
    # neeed to compute Cts and Ws(B,D)
    Ct = RF.Ct
    B = RF.B
    D = RF.D
    N = length(Ct.x)

    ## return 0 if request is out of bounds
    ## or in lower triangular part (k < j)
    (k < j || j == 0 || k > N+1) && return zero(RF)

    cj, sj  = vals(Ct[N+1-j])
    s = sj
    par = -1

    dk = D[k]
    w  =  _w(B, j,  k) * dk

    tot = w/s

    for i in (j+1):k
        w =  _w(B,  i,  k)
        ci, si = vals(Ct[N+1-i])
        s  = s*si
        tot += par *  cj*conj(ci)/s * w * dk
        par *= -1
    end

    tot
end

function Base.Matrix(RF::RFactorizationRankOne{T,S}) where {T,S}

    n = length(RF) + 1

    M = diagm(0 => ones(S, n))
    e1 = vcat(1, zeros(S, n-1))
    en1 = vcat(zeros(S,n-1), 1)
    en = vcat(zeros(S,n-2), 1, 0)

    ## compute yt
    Ct, B = RF.Ct, RF.B
    D = Matrix(RF.D)
    D = D * M
    rho = (en1' * (Ct * M) * e1)
    yt = -(1/rho * en1' * (Ct*(B*D)))

    # Compute R
    R =   (Ct * (B * D)) +  (Ct * (e1 * yt))

end

## Pass a rotator through Rfactorization from left or right; return modified U
function passthrough!(RF::RFactorizationRankOne, U::AbstractRotator)

    U = passthrough!(RF.D, U)
    U = passthrough!(RF.B, U)
    U = passthrough!(RF.Ct, U)

    U

end

function passthrough!(U::AbstractRotator, RF::RFactorizationRankOne)

    U = passthrough!(U, RF.Ct)
    U = passthrough!(U, RF.B)
    U = passthrough!(U, RF.D)

    U

end

##################################################
##
## This is useful if a unitary Hessenberg matrix is factored into rotators:
## Un Un_1 ... U2 U1 * Q = R so that Q = U1' U2' ... Un' * R
## R factor will have this structure
"""

For the case where the QR decomposion has R as a diagonal matrix that is unitary

"""
struct RFactorizationUnitaryDiagonal{T, S} <: AbstractRFactorization{T, S}
D::SparseDiagonal{S}
RFactorizationUnitaryDiagonal(D::SparseDiagonal{S}) where {S} = new{real(S),S}(D)
function RFactorizationUnitaryDiagonal(xs::Vector{S}) where {S}
    D = SparseDiagonal(xs)
    RFactorizationUnitaryDiagonal(D)
end
end
#XXX
Base.copy(RF::RFactorizationUnitaryDiagonal) = RFactorizationUnitaryDiagonal(RF.D.x)
Base.size(RF::RFactorizationUnitaryDiagonal) = size(RF.D)
Base.getindex(RF::RFactorizationUnitaryDiagonal{T, S}, i, j) where {T, S} = RF.D[i,j]
Base.Matrix(RF::RFactorizationUnitaryDiagonal) = Matrix(RF.D)
Base.length(RF::RFactorizationUnitaryDiagonal) = error("No dimension known")
passthrough!(RF::RFactorizationUnitaryDiagonal, U::AbstractRotator) = passthrough!(RF.D, U)
passthrough!(U::AbstractRotator, RF::RFactorizationUnitaryDiagonal) = passthrough!(U, RF.D)

simple_passthrough!(RF::RFactorizationUnitaryDiagonal, args...) = false

##################################################

# hold upper triangular matrix as full matrix
"""

For the case where the QR decomposition has R as a full, upper-triangular matrix

"""
struct RFactorizationUpperTriangular{T, S, Rt <: AbstractArray{S,2}} <: AbstractRFactorization{T, S}
R::Rt
RFactorizationUpperTriangular{T, S, Rt}(M) where {T, S, Rt}= new(M)
RFactorizationUpperTriangular(M::AbstractArray{S}) where {S} = RFactorizationUpperTriangular{real(S), S, typeof(M)}(M)
end

Base.copy(RF::RFactorizationUpperTriangular) = RFactorizationUpperTriangular(copy(RF.R))
Base.length(RF::RFactorizationUpperTriangular) = size(RF.R)[1]
Base.size(RF::RFactorizationUpperTriangular, args...) = size(RF.R, args...)
function Base.getindex(RF::RFactorizationUpperTriangular, i, j)
    if i > 0 && j > 0
        RF.R[i,j]
    else
        zero(eltype(RF.R))
    end
end
Base.Matrix(RF::RFactorizationUpperTriangular) = RF.R

## Pass rotator through RFactorization
## return modified U
function passthrough!(U::Rt, RF::RFactorizationUpperTriangular) where {Rt <: AbstractRotator}

    R = RF.R
    n = size(R, 2)

    c,s = vals(U)
    i = idx(U); j = i+1

    rik, rjk = R[i,i], R[j,i] # k = i
    R[i,i] = c *  rik + s *  rjk
    rji = - conj(s) * rik + conj(c) * rjk
    for k in j:n
        rik, rjk =  R[i,k],  R[j,k]
        R[i,k] = c * rik + s * rjk
        R[j,k] = -conj(s) * rik + conj(c) * rjk
    end

    c, s, r = givensrot(R[j,j], rji)
    c = conj(c)

    Vt = Rt(c, s,i)

    for k in 1:i
        rki, rkj =  R[k,i], R[k,j]
        R[k, i] = c * rki - conj(s) * rkj
        R[k, j] = s * rki + conj(c) * rkj
    end
    R[j,j] = s * rji + conj(c) * R[j,j]
    R[j,i] = zero(eltype(R))

    Vt'

end

function passthrough!(RF::RFactorizationUpperTriangular, V::Rt) where {Rt <: AbstractRotator}

    R = RF.R
    n = size(R)[2]

    c, s = vals(V)
    i = idx(V); j = i+1
    for k in 1:i
        rki, rkj =  R[k,i], R[k,j]
        R[k,i] = c * rki - conj(s) * rkj
        R[k,j] = s * rki + conj(c) * rkj
    end
    rji    = c * R[j,i] - conj(s) * R[j,j] # avoid R[j,i], as R could be upper triangular
    R[j,j] = s * R[j,i] + conj(c) *  R[j,j]

    c,s,r = givensrot(R[i,i], rji)
    Ut = Rt(c, s, i)

    # k = i
    rik, rjk = R[i,i], rji
    R[i,i] = c  * rik + s * rjk
    R[j,i] = 0

    for k in j:n
        rik, rjk = R[i,k], R[j,k]
        R[i, k]  = c * rik + s * rjk
        R[j, k]  = -conj(s) * rik + conj(c) * rjk
    end

    return Ut'

end

simple_passthrough!(RF::RFactorizationUpperTriangular, Us...) = false

##################################################
# Wrapper for I
struct RFactorizationIdentity{T, S} <: AbstractRFactorization{T, S}
RFactorizationIdentity{T, S}() where {T, S}= new()
end

Base.copy(RF::RFactorizationIdentity) = RF
Base.size(RF::RFactorizationIdentity) = (-1, -1)
Base.getindex(RF::RFactorizationIdentity{T, S}, i, j) where {T, S} = i == j ? one(S) : zero(S)
Base.Matrix(::RFactorizationIdentity) = I
Base.length(RF::RFactorizationIdentity) = error("No dimension known")
passthrough!(RF::RFactorizationIdentity, U::AbstractRotator) = U
passthrough!(U::AbstractRotator, RF::RFactorizationIdentity) = U
# this bypasses some more general one for speed
passthrough!(RF::RFactorizationIdentity, C::DescendingChain) = nothing
passthrough!(RF::RFactorizationIdentity, C::AscendingChain) = nothing
passthrough!(RF::RFactorizationIdentity, C::TwistedChain) = nothing
passthrough!(C::DescendingChain,RF::RFactorizationIdentity) = nothing
passthrough!(C::AscendingChain,RF::RFactorizationIdentity) = nothing
passthrough!(C::TwistedChain,RF::RFactorizationIdentity) = nothing

simple_passthrough!(RF::RFactorizationIdentity, Us...) = true


text:
PlotlyBase.trace_map(p::SyncPlot, axis) = trace_map(p.plot, axis)
JSON.lower(sp::SyncPlot) = sp.plot

PlotlyBase._is3d(p::SyncPlot) = _is3d(p.plot)

# subplot methods on syncplot
Base.hcat(sps::SyncPlot...) = SyncPlot(hcat(Plot[sp.plot for sp in sps]...))
Base.vcat(sps::SyncPlot...) = SyncPlot(vcat(Plot[sp.plot for sp in sps]...))
Base.hvcat(rows::Tuple{Vararg{Int}}, sps::SyncPlot...) =
    SyncPlot(hvcat(rows, Plot[sp.plot for sp in sps]...))

function PlotlyBase.add_recession_bands!(p::SyncPlot; kwargs...)
    new_shapes = add_recession_bands!(p.plot; kwargs...)
    relayout!(p, shapes=new_shapes)
    new_shapes
end

function mgrid(arrays...)
    lengths = collect(length.(arrays))
    uno = ones(Int, length(arrays))
    out = []
    for i in 1:length(arrays)
       repeats = copy(lengths)
       repeats[i] = 1

       shape = copy(uno)
       shape[i] = lengths[i]
       push!(out, reshape(arrays[i], shape...) .* ones(repeats...))
    end
    out
end


text:
"""
Implementation of types to describe the behavior of a material from a dependent
variable.
"""

#= #############################################################################
Material Properties "Law" Types
############################################################################# =#

"""
MaterialPropertyOfX:
describes a whatever (on x) depended material property (p)

We use several subtypes to allow for more accurate material property fitting
"""
abstract type MaterialPropertyOfX end;
abstract type BoundedMaterialPropertyOfX <: MaterialPropertyOfX end;

"""
ConstantMaterialPropertyOfX:
property that has not dependency on X
"""
struct ConstantMaterialPropertyOfX{T <: AbstractFloat} <: MaterialPropertyOfX
    name::String
    unit::String
    val::T

    ConstantMaterialPropertyOfX(val::T;
                                name::String = "empty",
                                unit::String = "empty"
                                ) where T <: AbstractFloat =
                                new{T}(name,unit,val)
end

"""
TabulatedMaterialPropertyOfX:
tabulated material property

Using tabulated material properties does not make many sense right now as the
computation of material properties itself is very cheap (even polymial)
"""
struct TabulatedMaterialPropertyOfX{T1 <: LinRange{<:AbstractFloat},
    T2 <: AbstractFloat, T3 <: Vector{<:AbstractFloat}} <: MaterialPropertyOfX
    name::String
    unit::String
    range::T1 # .len .stop .start
    rangestep::T2
    val_arr::T3

    TabulatedMaterialPropertyOfX(range::T1,
                                val_arr::T2;
                                name::String = "empty",
                                unit::String = "empty"
                                ) where {T1 <: LinRange{<:AbstractFloat},
                                T2 <: Vector{<:AbstractFloat}} =
                                new{T1, typeof(((range.stop-range.start)/range.len)),T2}(
                                name, unit, range,
                                ((range.stop - range.start) / range.len),
                                val_arr)
end

function TabulatedMaterialPropertyOfX(MP::MaterialPropertyOfX,
        range::LinRange{<:AbstractFloat}, bounded::Bool = true)
    # we get T of MaterialPropertyOfX
    # MPtype = typeof(MP)
    # MPdatafieldtype = MPtype.parameters[1]
    # val_arr = zeros(MPdatafieldtype, range.len)

    if bounded
        val_arr= pX_bounded.(MP, range)
    else
        val_arr = pX.(MP, range)
    end

    return TabulatedMaterialPropertyOfX(range, val_arr, name = MP.name, unit = MP.unit)
end

"""
LinearInterpMaterialPropertyOfX:
Implements type for linear interpolated materialproperties. Values are lineary
interpolated between die nearest values. If bounded is true no linear
extrapolation will be done and the last value will be used as constant
extrapolation.
"""
struct LinearInterpMaterialPropertyOfX{T <: AbstractFloat} <: MaterialPropertyOfX
    name::String
    unit::String
    bounded::Bool
    x::Vector{T}
    val::Vector{T}

    LinearInterpMaterialPropertyOfX(x::Vector{T}, val::Vector{T};
                                    bounded::Bool = false,
                                    name::String = "empty",
                                    unit::String = "empty"
                                    ) where T <: AbstractFloat =
                                    new{T}(name,unit,bounded,x,val)
end

"""
LinearMaterialPropertyOfX:
Implements material property based on linear law

p(x) = m*x + y1
"""
struct LinearMaterialPropertyOfX{T <: AbstractFloat} <: BoundedMaterialPropertyOfX
    name::String
    unit::String
    m::T
    y1::T
    bounds::Union{Nothing, Tuple{T,T}}

    LinearMaterialPropertyOfX(m::T, y1::T;
                            name::String = "empty",
                            unit::String = "empty",
                            bounds::Union{Nothing, Tuple{T,T}} = nothing
                            ) where T <: AbstractFloat =
                            new{T}(name,unit,m,y1, bounds)
end

"""
PolynomialMaterialPropertyOfX:
Implements material property based on polynomial law from till degree n

p(x) = a[1] + a[2]*x + a[3] * x^2 + ... + a[n] * x^n
"""
struct PolynomialMaterialPropertyOfX{T1 <: AbstractFloat, T2 <: Integer}  <: BoundedMaterialPropertyOfX
    name::String
    unit::String
    n::T2
    a::Vector{T1}
    bounds::Union{Nothing, Tuple{T1,T1}}

    PolynomialMaterialPropertyOfX(a::Vector{T};
                                  name::String = "empty",
                                  unit::String = "empty",
                                  bounds::Union{Nothing, Tuple{T,T}} = nothing
                                ) where T <: AbstractFloat =
                                new{T, typeof(length(a))}(name,unit,length(a),a,bounds)

end

"""
ReciprocalMaterialPropertyOfX:
Reciprocal material property

p(x) = m * 1/x + y1
"""
struct ReciprocalMaterialPropertyOfX{T <: AbstractFloat} <: BoundedMaterialPropertyOfX
    name::String
    unit::String
    m::T
    y1::T
    bounds::Union{Nothing, Tuple{T,T}}

    ReciprocalMaterialPropertyOfX(m::T, y1::T;
                                    name::String = "empty",
                                    unit::String = "empty",
                                    bounds::Union{Nothing, Tuple{T,T}} = nothing
                                    ) where T <: AbstractFloat =
                                    new{T}(name,unit,m,y1,bounds)
end

"""
ExponentialMaterialPropertyOfX:
Exponential material property

p(x) = a * exp(b*x) + y1
"""
struct ExponentialMaterialPropertyOfX{T <: AbstractFloat} <: BoundedMaterialPropertyOfX
    name::String
    unit::String
    a::T
    b::T
    y1::T
    bounds::Union{Nothing, Tuple{T,T}}

    ExponentialMaterialPropertyOfX(a::T,b::T, y1::T;
                                    name::String = "empty",
                                    unit::String = "empty",
                                    bounds::Union{Nothing, Tuple{T,T}} = nothing
                                    ) where T <: AbstractFloat =
                                    new{T}(name,unit,a,b,y1,bounds)
end

"""
LogarithmicMaterialPropertyOfX:
Logarithmic material property

p(x) = a * ln(x) + y1
"""
struct LogarithmicMaterialPropertyOfX{T <: AbstractFloat} <: BoundedMaterialPropertyOfX
    name::String
    unit::String
    a::T
    y1::T
    bounds::Union{Nothing, Tuple{T,T}}

    LogarithmicMaterialPropertyOfX(a::T, y1::T;
                                    name::String = "empty",
                                    unit::String = "empty",
                                    bounds::Union{Nothing, Tuple{T,T}} = nothing
                                    ) where T <: AbstractFloat =
                                    new{T}(name,unit,a,y1,bounds)
end

#= #############################################################################
Final Material Types
############################################################################# =#

"""
IncompressibleFluid:
MaterialPropertyOfX here X should represent the temperature

"""
struct IncompressibleFluidThermalMaterial{T1 <: MaterialPropertyOfX,
        T2 <: MaterialPropertyOfX,
        T3 <: MaterialPropertyOfX}

    name::String
    rho::T1
    cp_true::T2
    lambda::T3
end

"""
SolidThermalMaterial:
MaterialPropertyOfX here X should represent the temperature

"""
struct SolidThermalMaterial{T1 <: MaterialPropertyOfX,
        T2 <: MaterialPropertyOfX, T3 <: MaterialPropertyOfX,
        T4 <: MaterialPropertyOfX}

    name::String
    rho::T1
    cp_true::T2
    lambda::T3
    epsilon::T4
end


text:
"""
This module provides utilities to build neural networks with Flux,
along with a library of standard architectures.
"""
module FluxNets

export SimpleNet, SimpleNetHP, ResNet, ResNetHP

using ..Network
using Base: @kwdef
import ..GameInterface
import ..Util

using CUDAapi

# Import CuArrays only if CUDA is installed
if has_cuda()
  try
    using CuArrays
    @show CuArrays.usage_limit[]
    @eval const CUARRAYS_IMPORTED = true
  catch ex
    @warn(
      "CUDA is installed, but CuArrays.jl fails to load.",
      exception=(ex,catch_backtrace()))
    @eval const CUARRAYS_IMPORTED = false
  end
else
  @eval const CUARRAYS_IMPORTED = false
end

import Flux

if CUARRAYS_IMPORTED
  @eval begin
    CuArrays.allowscalar(false)
    on_gpu(::Type{<:Array}) = false
    on_gpu(::Type{<:CuArray}) = true
    on_gpu(x) = on_gpu(typeof(x))
  end
else
  @eval begin
    on_gpu(x) = false
  end
end

using Flux: relu, softmax
using Flux: Chain, Dense, Conv, BatchNorm, SkipConnection

#####
##### Flux Networks
#####

abstract type FluxNetwork{Game} <: AbstractNetwork{Game} end

function Base.copy(nn::Net) where Net <: FluxNetwork
  #new = Net(Network.hyperparams(nn))
  #Flux.loadparams!(new, Flux.params(nn))
  #return new
  return Base.deepcopy(nn)
end

Network.to_cpu(nn::FluxNetwork) = Flux.cpu(nn)

function Network.to_gpu(nn::FluxNetwork)
  CUARRAYS_IMPORTED && CuArrays.allowscalar(false)
  return Flux.gpu(nn)
end

# TODO: dirty hack to cope with the fact that Flux 0.10 removed Flux.testmode!
const TEST_MODE = IdDict()

function Network.set_test_mode!(nn::FluxNetwork, mode)
  TEST_MODE[nn] = mode
end

Network.convert_input(nn::FluxNetwork, x) =
  Network.on_gpu(nn) ? Flux.gpu(x) : x

Network.convert_output(nn::FluxNetwork, x) = Flux.cpu(x)

Network.params(nn::FluxNetwork) = Flux.params(nn)

function Network.train!(nn::FluxNetwork, loss, data, lr)
  optimizer = Flux.ADAM(lr)
  Flux.train!(loss, Flux.params(nn), data, optimizer)
end

regularized_child_leaves(l) = []
regularized_child_leaves(l::Flux.Dense) = [l.W]
regularized_child_leaves(l::Flux.Conv) = [l.weight]

# Reimplementation of what used to be Flux.prefor, does not visit leafs
function foreach_flux_node(f::Function, x, seen = IdDict())
  Flux.isleaf(x) && return
  haskey(seen, x) && return
  seen[x] = true
  f(x)
  for child in Flux.trainable(x)
    foreach_flux_node(f, child, seen)
  end
end

function Network.regularized_params(net::FluxNetwork)
  ps = Flux.Params()
  foreach_flux_node(net) do p
    for r in regularized_child_leaves(p)
      any(x -> x === r, ps) || push!(ps, r)
    end
  end
  return ps
end

function Network.gc(::FluxNetwork)
  CUARRAYS_IMPORTED || return
  GC.gc(true)
  CuArrays.reclaim()
end

#####
##### Common functions between two-head neural networks
#####

abstract type TwoHeadNetwork{G} <: FluxNetwork{G} end

function Network.forward(nn::TwoHeadNetwork, board)
  # TODO: eliminate this horror when Flux catches up
  # @eval Flux.istraining() = $(!get(TEST_MODE, nn, false))
  c = nn.common(board)
  v = nn.vhead(c)
  p = nn.phead(c)
  return (p, v)
end

# Flux.@functor does not work do to Network being parametric
function Flux.functor(nn::Net) where Net <: TwoHeadNetwork
  children = (nn.common, nn.vhead, nn.phead)
  constructor = cs -> Net(nn.hyper, cs...)
  return (children, constructor)
end

Network.hyperparams(nn::TwoHeadNetwork) = nn.hyper

Network.on_gpu(nn::TwoHeadNetwork) = on_gpu(nn.vhead[end].b)

#####
##### Utilities for the networks library
#####

linearize(x) = reshape(x, :, size(x)[end])

#####
##### Include networks library
#####

include("architectures/simplenet.jl")
include("architectures/resnet.jl")

end


text:
# The following was inspired by AtariAlgos.jl
# Modifications by Dennis Wilson @d9w

using Colors
using ImageCore
using ImageTransformations
export
    Game,
    close!,
    draw,
    get_inputs,
    get_rgb,
    get_float,
    get_int,
    get_bool,
    get_string

rom_directory(x::String) = joinpath(dirname(@__FILE__), "..", "deps", "rom_files", "$x.bin")

struct Game
    ale::ALEPtr
    width::Int
    height::Int
    actions::Array{Int32}
end

function Game(romfile::String, seed::Int64)
    ale = ALE_new()
    setInt(ale, "random_seed", Cint(seed))
    loadROM(ale, string(rom_directory(romfile)))
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    actions = getMinimalActionSet(ale)
    Game(ale, w, h, actions)
end

function close!(game::Game)
    ALE_del(game.ale)
end

function draw(game::Game)
    rawscreen = Array{Cuchar}(undef, game.width * game.height * 3)
    getScreenRGB(game.ale, rawscreen)
    colorview(RGB, Float64.(reshape(rawscreen/256.,
                                    (3, game.width, game.height))))';
end

function get_inputs(game::Game)
    screen = getScreen(game.ale)/(0xff*1.0)
    screen = reshape(screen, (game.width, game.height))'
    # imresize(screen, (42, 32))/256.
    screen
end

function get_rgb(game::Game)
    rawscreen = Array{Cuchar}(undef, game.width * game.height * 3)
    getScreenRGB(game.ale, rawscreen)
    rgb = Float64.(reshape(rawscreen/256.,(3, game.width, game.height)));
    [rgb[i,:,:]' for i in 1:3]
end

# Read access to A.L.E. configuration:
get_float(game::Game, key::String) = getFloat(game.ale, key)
get_int(game::Game, key::String) = getInt(game.ale, key)
get_bool(game::Game, key::String) = getBool(game.ale, key)
get_string(game::Game, key::String) = getString(game.ale, key)


text:
using FTD2XX
using GWinstekPSP
using Base.Test

uartconfig = UARTConfiguration(2400,8,1,"n",
                               readtimeout = 1000,   # in milliseconds
                               writetimeout = 1000)  # in milliseconds
io_psp603 = open(FT_Description("PSP-603"),uartconfig)

status = psp_status(io_psp603)
show(status)

voltage = psp_readvoltage(io_psp603)
current = psp_readcurrent(io_psp603)
power = psp_readpower(io_psp603)
vlimit = psp_readvlimit(io_psp603)
ilimit = psp_readilimit(io_psp603)
plimit = psp_readplimit(io_psp603)
(isrelayon,isoverheat,isknobfine,isknobunlock,isremote,islock) = psp_readflags(io_psp603)

voltage = psp_readvoltage(io_psp603)
@test voltage == status.voltage
current = psp_readcurrent(io_psp603)
@test current == status.current
power = psp_readpower(io_psp603)
@test power == status.power
vlimit = psp_readvlimit(io_psp603)
@test vlimit == status.vlimit
ilimit = psp_readilimit(io_psp603)
@test ilimit == status.ilimit
plimit = psp_readplimit(io_psp603)
@test plimit == status.plimit
(isrelayon,isoverheat,isknobfine,isknobunlock,isremote,islock) = psp_readflags(io_psp603)
@test isrelayon == status.isrelayon
@test isoverheat == status.isoverheat
@test isknobfine == status.isknobfine
@test isknobunlock == status.isknobunlock
@test isremote == status.isremote
@test islock == status.islock

psp_setrelayoff(io_psp603)
psp_incrementvoltage(io_psp603)
psp_decrementvoltage(io_psp603)
psp_incrementvlimit(io_psp603)
psp_decrementvlimit(io_psp603)
psp_incrementilimit(io_psp603)
psp_decrementilimit(io_psp603)
psp_incrementplimit(io_psp603)
psp_decrementplimit(io_psp603)

psp_setvlimitmax(io_psp603)
psp_setilimitmax(io_psp603)
psp_setplimitmax(io_psp603)
psp_setknobfine(io_psp603)
psp_setknobnormal(io_psp603)

psp_togglerelay(io_psp603)
psp_setrelayon(io_psp603)
psp_setrelayoff(io_psp603)

psp_savepreset(io_psp603)

psp_setvoltage(io_psp603,1)
psp_setvlimit(io_psp603,10.123456)
psp_setilimit(io_psp603,2)
psp_setvlimitmax(io_psp603)
psp_setplimit(io_psp603,100)
psp_setvoltage(io_psp603,11.12345)
psp_setvlimit(io_psp603,22)
psp_setilimit(io_psp603,1)

psp_setvoltage(io_psp603,0)
psp_setvlimitmax(io_psp603)
psp_setilimitmax(io_psp603)
psp_setplimitmax(io_psp603)
psp_savepreset(io_psp603)

# restore state
psp_setvoltage(io_psp603,status.voltage)
psp_setvlimit(io_psp603,status.vlimit)
psp_setilimit(io_psp603,status.ilimit)
psp_setplimit(io_psp603,status.plimit)

close(io_psp603)