text:
"""
  OptimizationProblem(...)
A OptimizationProblem (Minimization Problem), is used to represent the mathematical problem
of finding local minima of the given objective function. The problem is defined by `objective`
which is an appropriate objective type (for example `NonDiffed`, `OnceDiffed`, ...)
for the types of algorithm to be used. The constraints of the problem are encoded
in `constraints`. See the documentation for supported types of constraints
including convex sets, and more. It is possible to explicitly state that there
are bounds constraints and manifold constraints on the inputs.

Options are stored in `options` and must be an appropriate options type. See more information
about options using `?OptimizationOptions`.
"""
struct OptimizationProblem{O, B, M<:Manifold, C, I, X} <: AbstractProblem
    objective::O
    bounds::B
    manifold::M
    constraints::C
    mstyle::I
    initial::X
end
value(prob::OptimizationProblem, x) = value(prob.objective, x)
batched_value(prob::OptimizationProblem, X) = batched_value(prob.objective, X)
batched_value(prob::OptimizationProblem, F, X) = batched_value(prob.objective, F, X)
upto_gradient(prob::OptimizationProblem, ∇f, x) = upto_gradient(prob.objective, ∇f, x)
upto_hessian(prob::OptimizationProblem, ∇f, ∇²f, x) = upto_hessian(prob.objective, ∇f, ∇²f, x)
_manifold(prob::OptimizationProblem) = prob.manifold
lowerbounds(mp::OptimizationProblem) = mp.bounds[1]
upperbounds(mp::OptimizationProblem) = mp.bounds[2]
hasbounds(mp::OptimizationProblem) = mp.bounds isa Tuple
bounds(mp::OptimizationProblem) = (lower=lowerbounds(mp), upper=upperbounds(mp))
isboundedonly(::OptimizationProblem{<:Any, <:Nothing, <:Any, <:Nothing}) = false
isboundedonly(::OptimizationProblem{<:Any, <:Nothing, <:Any, <:Any}) = false
isboundedonly(::OptimizationProblem{<:Any, <:Any, <:Any, <:Nothing}) = true
#value(p::OptimizationProblem, args...) = p.objective(args...)
constraints(p::OptimizationProblem, args...) = p.constraints(args...)
OptimizationProblem(; obj=nothing, bounds=nothing, manifold=Euclidean(0), constraints=nothing, inplace=true, initial_x=nothing) =
  OptimizationProblem(obj, bounds, manifold, constraints, inplace===true ? InPlace() : OutOfPlace(), initial_x)
OptimizationProblem(obj; inplace=true) = OptimizationProblem(obj, nothing, Euclidean(0), nothing, inplace===true ? InPlace() : OutOfPlace(), nothing)

struct ConvergenceInfo{Ts, T, O}
  solver::Ts
  info::T
  options::O
end
function Base.show(io::IO, ci::ConvergenceInfo)
  opt = ci.options
  info = ci.info

  println(io, "Results of minimization\n")
  println(io, "* Algorithm:")
  println(io, "  $(summary(ci.solver))")
  println(io)
  println(io, "* Candidate solution:")
  println(io, "  Final objective value:    $(@sprintf("%.2e", ci.info.minimum))")
  if haskey(info, :∇fz)
    println(io, "  Final gradient norm:      $(@sprintf("%.2e", opt.g_norm(info.∇fz)))")
    if haskey(info, :prob) && hasbounds(info.prob)
      ρP = opt.g_norm(info.minimizer.-clamp.(info.minimizer.-info.∇fz, info.prob.bounds...))
      println(io, "  Final projected gradient norm:  $(@sprintf("%.2e", ρP))")
    end
  end
  if haskey(info, :temperature)
    println(io, "  Final temperature:        $(@sprintf("%.2e", ci.info.temperature))")
  end
  println(io)
  println(io, "  Initial objective value:  $(@sprintf("%.2e", ci.info.f0))")
  if haskey(info, :∇f0)
    println(io, "  Initial gradient norm:    $(@sprintf("%.2e", info.∇f0))")
  end
  println(io)
  println(io, "* Convergence measures")
  if isa(ci.solver, NelderMead)
    nm_converged(r) = 0.0
    println(io, "  √(Σ(yᵢ-ȳ)²)/n         = $(@sprintf("%.2e", info.nm_obj)) <= $(@sprintf("%.2e", opt.nm_tol)) ($(info.nm_obj<=opt.nm_tol))")
  elseif isa(ci.solver, SimulatedAnnealing)
  else
    if haskey(info, :ρs)
      println(io, "  |x - x'|              = $(@sprintf("%.2e", info.ρs)) <= $(@sprintf("%.2e", opt.x_abstol)) ($(info.ρs<=opt.x_abstol))")
      println(io, "  |x - x'|/|x|          = $(@sprintf("%.2e", info.ρs/info.ρx)) <= $(@sprintf("%.2e", opt.x_reltol)) ($(info.ρs/info.ρx <= opt.x_reltol))")
    end
    if isfinite(opt.f_limit)
      println(io, "  |f(x')|               = $(@sprintf("%.2e", info.minimum)) <= $(@sprintf("%.2e", opt.f_limit)) ($(info.minimum<=opt.f_limit))")
    end
    if haskey(info, :fx)
      Δf = abs(info.fx-info.minimum)
      println(io, "  |f(x) - f(x')|        = $(@sprintf("%.2e", Δf)) <= $(@sprintf("%.2e", opt.f_abstol)) ($(Δf<=opt.f_abstol))")
      println(io, "  |f(x) - f(x')|/|f(x)| = $(@sprintf("%.2e", Δf/abs(info.fx))) <= $(@sprintf("%.2e", opt.f_reltol)) ($(Δf/abs(info.fx)<=opt.f_reltol))")
    end
    if haskey(info, :∇fz)
      ρ∇f = opt.g_norm(info.∇fz)
      if haskey(info, :prob) && hasbounds(info.prob)
        println(io, "  |x - P(x - g(x))|     = $(@sprintf("%.2e", ρP)) <= $(@sprintf("%.2e", opt.g_abstol)) ($(ρP<=opt.g_abstol))")
      end
      println(io, "  |g(x)|                = $(@sprintf("%.2e", ρ∇f)) <= $(@sprintf("%.2e", opt.g_abstol)) ($(ρ∇f<=opt.g_abstol))")
      println(io, "  |g(x)|/|g(x₀)|        = $(@sprintf("%.2e", ρ∇f/info.∇f0)) <= $(@sprintf("%.2e", opt.g_reltol)) ($(ρ∇f/info.∇f0<=opt.g_reltol))")
    end
    if haskey(info, :Δ)
      Δmin = ci.solver.Δupdate.Δmin isa Nothing ? 0 : ci.solver.Δupdate.Δmin
      if isnan(info.Δ)
        println(io, "  Δ                     = $(@sprintf("%.2e", info.Δ)) (updated radius was not finite)")
      else
        Δtest = info.Δ<=Δmin
        println(io, "  Δ                     = $(@sprintf("%.2e", info.Δ)) <= $(@sprintf("%.2e", Δmin)) ($Δtest)")
      end
    end
    if haskey(info, :prob) && hasbounds(info.prob)
      if any(iszero, info.minimizer.-info.prob.bounds[1]) || any(iszero, info.prob.bounds[2].-info.minimizer)
        println(io, "\n  !!! Solution is at the boundary !!!")
      end
    end
  end
  println(io)
  println(io, "* Work counters")
  println(io, "  Seconds run:   $(@sprintf("%.2e", info.time))")
  println(io, "  Iterations:    $(info.iter)")
end

struct OptimizationOptions{T1, T2, T3, T4, Txn, Tgn, Tlog}
  x_abstol::T1
  x_reltol::T1
  x_norm::Txn
  g_abstol::T2
  g_reltol::T2
  g_norm::Tgn
  f_limit::T3
  f_abstol::T3
  f_reltol::T3
  nm_tol::T3
  maxiter::T4
  show_trace::Bool
  logger::Tlog
end

OptimizationOptions(; x_abstol=0.0, x_reltol=0.0, x_norm=x->norm(x, Inf),
             g_abstol=1e-8, g_reltol=0.0, g_norm=x->norm(x, Inf),
             f_limit=-Inf, f_abstol=0.0, f_reltol=0.0,
             nm_tol=1e-8, maxiter=10000, show_trace=false, logger=show_trace ? ConsoleLogger() : NullLogger()) =
  OptimizationOptions(x_abstol, x_reltol, x_norm,
             g_abstol, g_reltol, g_norm,
             f_limit, f_abstol, f_reltol,
             nm_tol,
             maxiter, show_trace, logger)

struct MinResults{Tr, Tc<:ConvergenceInfo, Th, Ts, To}
  res::Tr
  conv::Tc
  history::Th
  solver::Ts
  options::To
end
convinfo(mr::MinResults) = mr.conv
function converged(MinResults)
end
function Base.show(io::IO, mr::MinResults)
  println(io, "* Status: $(any(converged(mr)))")
  println(io)
  println(io, "* Candidate solution")
  println(io, "  MinimizerMinOpt: ", minimizer(mr))
  println(io, "  Minimum:   ", minimum(mr))
  println(io)
  println(io, "* Found with")
  println(io, "  Algorithm: ", summary(mr.solver))
  println(io, "  Initial point: ", initial_point(mr))
  println(io, "  Initial value: ", initial_value(mr))

  println(" * Trace stored: ", has_history(mr))

  println(io)
  println(io, " * Convergence measures\n")
  show(io, convinfo(mr))
end
#=
* Work counters
  Seconds run:   0  (vs limit Inf)
  Iterations:    16
  f(x) calls:    53
  ∇f(x) calls:   53
=#

function prepare_variables(prob, approach, x0, ∇fz, B)
    objective = prob.objective
    z = x0
    x = copy(z)
    if isboundedonly(prob)
        !any(clamp.(x0, lowerbounds(prob), upperbounds(prob)) .!= x0) || error("Initial guess not in the feasible region")
    end

    if isa(B, Nothing)  # didn't provide a B
        if modelscheme(approach) isa GradientDescent
            # We don't need to maintain a dense matrix for Gradient Descent
            B = I
        elseif modelscheme(approach) isa LBFGS
            B = nothing
        else
            # Construct a matrix on the correct form and of the correct type
            # with the content of I_{n,n}
            B = I + abs.(0*x*x')
        end
    end
    # first evaluation
    if isa(modelscheme(approach), Newton)
        fz, ∇fz, B = upto_hessian(prob, ∇fz, B, x)
    else
        fz, ∇fz = upto_gradient(prob, ∇fz, x)
    end
    fx = copy(fz)
    ∇fx = copy(∇fz)
    Pg = ∇fz
    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=B, Pg=Pg)
end

function g_converged(∇fz, ∇f0, options)
  g_converged = options.g_norm(∇fz) ≤ options.g_abstol
  g_converged = g_converged || options.g_norm(∇fz) ≤ ∇f0*options.g_reltol
  g_converged = g_converged || any(isnan, ∇fz)
  return g_converged
end

function x_converged(x, z, options)
  x_converged = false
  if x !== nothing # if not calling from initial_converged
    y = x .- z
    ynorm = options.x_norm(y)
    x_converged = x_converged || ynorm ≤ options.x_abstol
    x_converged = x_converged || ynorm ≤ options.x_norm(x)*options.x_reltol
  end
  x_converged = x_converged || any(isnan, z)
  return x_converged
end
function f_converged(fx, fz, options)
  f_converged = false
  if fx !== nothing # if not calling from initial_converged
    y = fx - fz
    ynorm = abs(y)
    f_converged = f_converged || ynorm ≤ options.f_abstol
    f_converged = f_converged || ynorm ≤ abs(fx)*options.f_reltol
  end
  f_converged = f_converged || fz ≤ options.f_limit
  f_converged = f_converged || isnan(fz)
  return f_converged
end
function converged(approach, objvars, ∇f0, options, skip=nothing, Δkp1=nothing)
  x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
  xcon, gcon, fcon = x_converged(x, z, options), g_converged(∇fz, ∇f0, options), f_converged(fx, fz, options)
  if approach isa TrustRegion
    Δcon = isnan(Δkp1)
    Δcon = Δcon || (!(approach.Δupdate.Δmin isa Nothing) && Δkp1 < approach.Δupdate.Δmin)
    if skip == true
      # special logic for region reduction here
      return false, false, false, Δcon
    end
    return xcon, gcon, fcon, Δcon
  end
  return xcon, gcon, fcon
end

function initial_converged(approach, objvars, ∇f0, options, skip=nothing, Δkp1=nothing)
  x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
  objvars = (x=nothing, fx=nothing, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=B, Pg=Pg)
  converged(approach, objvars, ∇f0, options, skip, Δkp1)
end


text:
using AbstractPlotting: PointTrans, xyz_boundingbox, apply_transform

@testset "Basic transforms" begin
    function fpoint2(x::Point2)
        return Point2f0(x[1] + 10, x[2] - 77)
    end

    function fpoint3(x::Point3)
        return Point3f0(x[1] + 10, x[2] - 77, x[3] /  4)
    end
    trans2 = PointTrans{2}(fpoint2)
    trans3 = PointTrans{3}(fpoint3)
    points2 = [Point2f0(0, 0), Point2f0(0, 1)]
    bb = xyz_boundingbox(trans2, points2)
    @test bb == Rect(Vec3f0(10, -77, 0), Vec3f0(0, 1, 0))
    bb = xyz_boundingbox(trans3, points2)
    @test bb == Rect(Vec3f0(10, -77, 0), Vec3f0(0, 1, 0))

    points3 = [Point3f0(0, 0, 4), Point3f0(0, 1, -8)]
    bb = xyz_boundingbox(trans2, points3)
    @test bb == Rect(Vec3f0(10, -77, -8), Vec3f0(0, 1, 12))
    bb = xyz_boundingbox(trans3, points3)
    @test bb == Rect(Vec3f0(10, -77, -2.0), Vec3f0(0, 1, 3.0))

    @test apply_transform(trans2, points2) == fpoint2.(points2)
    @test apply_transform(trans3, points3) == fpoint3.(points3)

    @test_throws ErrorException PointTrans{2}(x::Int -> x)
    @test_throws ErrorException PointTrans{3}(x::Int -> x)
end

@testset "Tuple and identity transforms" begin
    t1 = sqrt
    t2 = (sqrt, log)
    t3 = (sqrt, log, log10)

    p2 = Point(2.0, 5.0)
    p3 = Point(2.0, 5.0, 4.0)

    @test apply_transform(identity, p2) == p2
    @test apply_transform(identity, p3) == p3

    @test apply_transform(t1, p2) == Point(sqrt(2.0), sqrt(5.0))
    @test apply_transform(t1, p3) == Point(sqrt(2.0), sqrt(5.0), sqrt(4.0))

    @test apply_transform(t2, p2) == Point2f0(sqrt(2.0), log(5.0))
    @test apply_transform(t2, p3) == Point3f0(sqrt(2.0), log(5.0), 4.0)

    @test apply_transform(t3, p3) == Point3f0(sqrt(2.0), log(5.0), log10(4.0))

    i2 = (identity, identity)
    i3 = (identity, identity, identity)
    @test apply_transform(i2, p2) == p2
    @test apply_transform(i3, p3) == p3

    # test that identity gives back exact same arrays without copying
    p2s = Point2f0[(1, 2), (3, 4)]
    @test apply_transform(identity, p2s) === p2s
    @test apply_transform(i2, p2s) === p2s
    @test apply_transform(i3, p2s) === p2s

    p3s = Point3f0[(1, 2, 3), (3, 4, 5)]
    @test apply_transform(identity, p3s) === p3s
    @test apply_transform(i2, p3s) === p3s
    @test apply_transform(i3, p3s) === p3s

    @test apply_transform(identity, 1) == 1
    @test apply_transform(i2, 1) == 1
    @test apply_transform(i3, 1) == 1

    @test apply_transform(identity, 1..2) == 1..2
    @test apply_transform(i2, 1..2) == 1..2
    @test apply_transform(i3, 1..2) == 1..2

    pa = Point2f0(1, 2)
    pb = Point2f0(3, 4)
    r2 = FRect2D(pa, pb .- pa)
    @test apply_transform(t1, r2) == FRect2D(apply_transform(t1, pa), apply_transform(t1, pb) .- apply_transform(t1, pa) )
end

ext: jl


text:
function default_convex_hull_algorithm(points)
    if length(points) != 0 && length(first(points)) == 2
        return "monotone_chain"
    else
        return nothing
    end
end

"""
    convex_hull(points::Vector{VN};
                [algorithm]=default_convex_hull_algorithm(points),
                [backend]=nothing
                )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}

Compute the convex hull of the given points.

### Input

- `points`    -- list of vectors
- `algorithm` -- (optional, default: depends on the dimension) the convex hull
                 algorithm, see valid options below
- `backend`   -- (optional, default: `"nothing"`) polyhedral computation backend
                 for higher-dimensional point sets

### Output

The convex hull as a list of vectors with the coordinates of the points.

### Algorithm

A pre-processing step treats the cases with `0`, `1` and `2` points in any dimension.
For more than `3` points, the algorithm used depends on the dimension.

For the one-dimensional case we return the minimum and maximum points, in that order.

The two-dimensional case is handled with a planar convex hull algorithm. The
following algorithms are available:

* `"monotone_chain"`        -- compute the convex hull of points in the plane
                               using Andrew's monotone chain method
* `"monotone_chain_sorted"` -- the same as `"monotone_chain"` but assuming that
                               the points are already sorted in counter-clockwise
                               fashion

See the reference docstring of each of those algorithms for details.

The higher dimensional case is treated using the concrete polyhedra library
`Polyhedra`, that gives access to libraries such as `CDDLib` and `ConvexHull.jl`.
These libraries can be chosen from the `backend` argument.

### Notes

For the in-place version use `convex_hull!` instead of `convex_hull`.

### Examples

Compute the convex hull of a random set of points:

```jldoctest ch_label
julia> points = [randn(2) for i in 1:30]; # 30 random points in 2D

julia> hull = convex_hull(points);

julia> typeof(hull)
Array{Array{Float64,1},1}
```

Plot both the random points and the computed convex hull polygon:

```jldoctest ch_label
julia> using Plots;

julia> plot([Tuple(pi) for pi in points], seriestype=:scatter);

julia> plot!(VPolygon(hull), alpha=0.2);
```
"""
function convex_hull(points::Vector{VN};
                     algorithm=default_convex_hull_algorithm(points),
                     backend=nothing
                     )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
    return convex_hull!(copy(points), algorithm=algorithm, backend=backend)
end

function convex_hull!(points::Vector{VN};
                      algorithm=default_convex_hull_algorithm(points),
                      backend=nothing)::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}

    m = length(points)

    # =====================
    # Treat special cases
    # =====================

    # zero or one point
    if m == 0 || m == 1
        return points
    end

    n = length(first(points))

    # two points
    if m == 2
        if n == 1
            p1, p2 = points[1], points[2]
            if p1 == p2  # check for redundancy
                pop!(points)
            elseif p1[1] > p2[1]
                points[1], points[2] = p2, p1
            end
        elseif n == 2
            # special case, see #876
            p1, p2 = points[1], points[2]
            if p1 == p2  # check for redundancy
                pop!(points)
            elseif p1 <= p2
                nothing
            else
                points[1], points[2] = p2, p1
            end
        end
        return points
    end

    # ===========================
    # Dispatch for general cases
    # ===========================

    # _convex_hull_x can assume that there are at least three points
    if n == 1
        return _convex_hull_1d!(points)
    elseif n == 2
        return _convex_hull_2d!(points, algorithm=algorithm)
    else
        return _convex_hull_nd!(points, backend=backend)
    end
end

function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
    points[1:2] = [minimum(points), maximum(points)]
    return resize!(points, 2)
end

function _convex_hull_nd!(points::Vector{VN};
                          backend=nothing
                          )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
    V = VPolytope(points)
    if backend == nothing
        backend = default_polyhedra_backend(V, N)
    end
    Vch = remove_redundant_vertices(V, backend=backend)
    m = length(Vch.vertices)
    points[1:m] = Vch.vertices
    return resize!(points, m)
end

function _convex_hull_2d!(points::Vector{VN};
                          algorithm::String="monotone_chain"
                          )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}
    if algorithm == "monotone_chain"
        return monotone_chain!(points)
    elseif algorithm == "monotone_chain_sorted"
        return monotone_chain!(points, sort=false)
    else
        error("the convex hull algorithm $algorithm is unknown")
    end
end

"""
    right_turn(O::AbstractVector{N}, A::AbstractVector{N}, B::AbstractVector{N}
              )::N where {N<:Real}

Determine if the acute angle defined by the three points `O`, `A`, `B` in the
plane is a right turn (counter-clockwise) with respect to the center `O`.

### Input

- `O` -- 2D center point
- `A` -- 2D one point
- `B` -- 2D another point

### Output

Scalar representing the rotation.

### Algorithm

The [cross product](https://en.wikipedia.org/wiki/Cross_product) is used to
determine the sense of rotation. If the result is 0, the points are collinear;
if it is positive, the three points constitute a positive angle of rotation
around `O` from `A` to `B`; otherwise they constitute a negative angle.
"""
@inline function right_turn(O::AbstractVector{N},
                            A::AbstractVector{N},
                            B::AbstractVector{N})::N where {N<:Real}
    return (A[1] - O[1]) * (B[2] - O[2]) - (A[2] - O[2]) * (B[1] - O[1])
end

"""
    monotone_chain!(points::Vector{VN}; sort::Bool=true
                   )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}

Compute the convex hull of points in the plane using Andrew's monotone chain
method.

### Input

- `points` -- list of 2D vectors; is sorted in-place inside this function
- `sort`   -- (optional, default: `true`) flag for sorting the vertices
              lexicographically; sortedness is required for correctness

### Output

List of vectors containing the 2D coordinates of the corner points of the
convex hull.

### Notes

For large sets of points, it is convenient to use static vectors to get
maximum performance. For information on how to convert usual vectors
into static vectors, see the type `SVector` provided by the
[StaticArrays](http://juliaarrays.github.io/StaticArrays.jl/stable/)
package.

### Algorithm

This function implements Andrew's monotone chain convex hull algorithm to
construct the convex hull of a set of ``n`` points in the plane in
``O(n \\log n)`` time.
For further details see
[Monotone chain](https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain)
"""
function monotone_chain!(points::Vector{VN}; sort::Bool=true
                        )::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}

    @inline function build_hull!(semihull, iterator, points, zero_N)
        @inbounds for i in iterator
            while length(semihull) >= 2 &&
                    (right_turn(semihull[end-1], semihull[end], points[i])
                         <= zero_N)
                pop!(semihull)
            end
            push!(semihull, points[i])
        end
    end

    if sort
        # sort the rows lexicographically (requires a two-dimensional array)
        # points = sortrows(hcat(points...)', alg=QuickSort) # out-of-place version
        sort!(points, by=x->(x[1], x[2]))                    # in-place version
    end

    zero_N = zero(N)

    # build lower hull
    lower = Vector{VN}()
    build_hull!(lower, axes(points)[1], points, zero_N)

    # build upper hull
    upper = Vector{VN}()
    build_hull!(upper, reverse(axes(points)[1]), points, zero_N)

    # remove the last point of each segment because they are repeated
    copyto!(points, @view(lower[1:end-1]))
    copyto!(points, length(lower), @view(upper[1:end-1]))
    return resize!(points, length(lower) + length(upper) - 2)
end


text:
using Revise
using Pkg

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using BayesMNIST

sleep(0.4)
@time BayesMNIST.modelmlp()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using VectorSphericalWaves

m, n, kr_r, kr_i, θ, ϕ = 1, 1, 1e6, 1e2, 0.2, 0.3;
kr = complex(kr_r, kr_i);

M_mn_wave(m, n, kr, θ, ϕ; kind="regular");
N_mn_wave(m, n, kr, θ, ϕ; kind="regular");

# BigFloat
m, n, kr_r, kr_i, θ, ϕ = 1, 1, BigFloat(1e6), BigFloat(1e2), BigFloat(0.2), BigFloat(0.3);
kr = complex(kr_r, kr_i);

M_mn_wave(m, n, kr, θ, ϕ; kind="regular");
N_mn_wave(m, n, kr, θ, ϕ; kind="regular");

# SVector comparison
m, n, kr_r, kr_i, θ, ϕ = 1, 1, 1e6, 1e2, 0.2, 0.3;
kr = complex(kr_r, kr_i);

θ_array, ϕ_array = LinRange(0, π, 1000), LinRange(0, π, 1000)

println("M and N returning Arrays")
@time M_mn_wave.(m, n, kr, θ_array, ϕ_array; kind="regular");
@time N_mn_wave.(m, n, kr, θ_array, ϕ_array; kind="regular");

println("M and N returning SVector, you should notice the sligth speedup")
@time VectorSphericalWaves.M_mn_wave_SVector.(m, n, kr, θ_array, ϕ_array; kind="regular");
@time VectorSphericalWaves.N_mn_wave_SVector.(m, n, kr, θ_array, ϕ_array; kind="regular");


text:
using SDatasets
using Base.Test

d = data_ball()

@test size(d) == (1, 2)

d = data_ball(10)

@test size(d) == (10, 2)

@test sum(sum(d.*d,2).>1) == 0 # test radius


text:
module test_loltools_gameclient_replay

using Test
using LOLTools.GameClient
using HTTP, Bukdu

function mock_action(verb, endpoint, path, headers=[], query=nothing, body=nothing)
    if query !== nothing
        path = string(merge(HTTP.URI(path), query=query))
    end
    if verb === HTTP.post
        bukdu_verb = post
    else
        bukdu_verb = get
    end
    Router.call(bukdu_verb, path, headers, something(body, UInt8[])).resp
end

struct ReplayController <: ApplicationController
    conn::Conn
end

function replay_game(c::ReplayController)
    render(JSON, (processID = 2544,))
end

function replay_particles(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_particles_post(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_playback(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_playback_post(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_render(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_render_post(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_recording(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_recording_post(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_sequence(c::ReplayController)
    render(JSON, (a=1,))
end

function replay_sequence_post(c::ReplayController)
    render(JSON, (a=1,))
end

function swagger_v2_swagger(c::ReplayController)
    render(JSON, (a=1,))
end

function swagger_v3_openapi(c::ReplayController)
    render(JSON, (a=1,))
end

routes() do
    get("/replay/game", ReplayController, replay_game)
    get("/replay/particles", ReplayController, replay_particles)
    get("/replay/playback", ReplayController, replay_playback)
    get("/replay/render", ReplayController, replay_render)
    get("/replay/recording", ReplayController, replay_recording)
    get("/replay/sequence", ReplayController, replay_sequence)
    get("/swagger/v2/swagger.json", ReplayController, swagger_v2_swagger)
    get("/swagger/v3/openapi.json", ReplayController, swagger_v3_openapi)
    post("/replay/particles", ReplayController, replay_particles_post)
    post("/replay/playback", ReplayController, replay_playback_post)
    post("/replay/render", ReplayController, replay_render_post)
    post("/replay/recording", ReplayController, replay_recording_post)
    post("/replay/sequence", ReplayController, replay_sequence_post)
end

@test (processID = 2544,) == GameClient.replay_game(; action=mock_action)
particles = GameClient.replay_particles(; action=mock_action)
GameClient.replay_particles(particles; action=mock_action)
playback = GameClient.replay_playback(; action=mock_action)
GameClient.replay_playback(playback; action=mock_action)
render2 = GameClient.replay_render(; action=mock_action)
GameClient.replay_render(render2; action=mock_action)
recording = GameClient.replay_recording(; action=mock_action)
GameClient.replay_recording(recording; action=mock_action)
sequence = GameClient.replay_sequence(; action=mock_action)
GameClient.replay_sequence(sequence; action=mock_action)
GameClient.swagger_v2_swagger(; action=mock_action)
GameClient.swagger_v3_openapi(; action=mock_action)

end # modu test_loltools_gameclient_replay


text:
str = collect(readline())
n = parse(Int, readline())
for _ in 1:n
    op = split(readline())
    a, b = parse.(Int, op[2:3])
    if op[1] == "print"
        for i in a+1:b+1
            print(str[i])
        end
        println()
    elseif op[1] == "reverse"
        reverse!(view(str, a+1:b+1))
    else
        substr = view(str, a+1:b+1)
        substr .= collect(op[4])
    end
end


text:
## path to TRPL folder --> set it before running the code!
path_to_trpl = "D:\\Projects\\TRPL" # some path to the TRPL folder (just an example here)

## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## loading packages
using OrdinaryDiffEq, LinearAlgebra, BlackBoxOptim, ForwardDiff,
	  Printf, Distributions, DataFrames, CSV, XLSX, NLopt

## directory
cd(joinpath(path_to_trpl,"Data\\Real"))

## output results? (true or false)
out = false

## parameters measurement data
name_measurement = "sample1"
if name_measurement == "sample1"
    N_DA = 5e16 # sample1: N_DA = 5e16
elseif name_measurement == "sample2"
    N_DA = 1e17 # sample2: N_DA = 1e17
end

## confidence level (1-α)
alpha = 0.05 # corresponds to a 95% confidence level
q = quantile(Normal(),1-alpha/2)

## load measurement data
data_measurement = Matrix(CSV.read(name_measurement*".csv",DataFrame,skipto=2))

## function for data processing (data structre is explained in README.md)
function data_processing(data)
    t = data[4:end,1]
    num_traj = size(data,2)-1
	num_pars = 10
    dims = num_pars + num_traj
    P = round.(data[1,2:end],sigdigits=3)
    ND = 10 .^ data[2,2:end]
    background = data[3,2:end]
    IPL = data[4:end,2:num_traj+1]
    tspan = (t[1], t[end])
    N = length(t)
    idx_IPL = zeros(Bool,size(IPL))
    for i = 1:num_traj
		idx_IPL[:,i] = IPL[:,i] .> 0
    end
    C_scaling = maximum(IPL)
    return t, tspan, num_traj, num_pars, dims, IPL, idx_IPL, P, ND, background, N, C_scaling
end

## process measurement data
t, tspan, num_traj, num_pars, dims, IPL, idx_IPL, P, ND, background, N, C_scaling = data_processing(data_measurement)

## definition rate equations
function rate_equations(du, u, p, t)
    du[1] = -((1 - u[2]) / p[3]) * u[1] - ((1 - u[3]) / p[4]) * u[1] + (1 / p[5]) * u[3] * 1 / p[7] -
            ((1 + u[1]) / p[1]) * u[1] - ((1 + u[1]) / (p[2] * (1 + p[9] * u[1]))) * u[1]
    du[2] = ((1 - u[2]) / (p[3])) * p[8] * u[1]
    du[3] = ((1 - u[3]) / p[4]) * p[7] * u[1] - u[3] / p[5] - u[3] / p[6]
end

## some inital parameters
# 		         τ_r0,  τ_nr0, τ_d, τ_0,  τ_1, τ_l,  N_l,  N_D,    C,  r , η_0s
init_par = vcat([   1,      1,   1,   1,    1,   1,  1e8,  1e8,  1e4, 10], repeat([1e-6],num_traj))

## ODEProblem definition with temporary u0 and parameters
u0 = [0.1, 0, 0]
prob = ODEProblem(rate_equations, u0, tspan, init_par)

## function to generate the model for plotting
function generateipl(par)
    function prob_func(prob, i, repeat)
        remake(prob, u0 = [par[num_pars+i], 0, 0])
    end
    ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
    sol = Array(solve(ensemble_simu, Rosenbrock23(), p = par[1:num_pars-1], EnsembleThreads(),
		      		  trajectories = num_traj, saveat = t, save_idxs = [1],
		      		  reltol = 1e-6, abstol = 1e-8))[1,:,:]
    ipl = zeros(size(sol))
    for i = 1:num_traj
    	ipl[:,i] = ((C_scaling * par[num_pars] / ND[i]) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
    end
    return ipl
end

## definition loss function
function loss(par)
    l = 0.0
    function prob_func(prob, i, repeat)
        remake(prob, u0 = [par[num_pars+i], 0, 0])
    end
    ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
    sol = Array(solve(ensemble_simu, Rosenbrock23(), p = par[1:num_pars-1], EnsembleThreads(),
		      		  trajectories = num_traj, saveat = t,save_idxs = [1],
		      		  reltol = 1e-6,abstol = 1e-8))[1,:,:]
    for i = 1:num_traj
		model = ((C_scaling * par[num_pars] / ND[i]) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
		idx_used = idx_IPL[:,i] .& (model .> 0) # we ignore zeros for the log
		l += sum((model[idx_used]) .- (IPL[idx_used,i] .* log.(model[idx_used])))
    end
    return l
end

## definition of loss function for NLopt
function loss_nlopt(par::Vector,grad::Vector)
    if length(grad) > 0
		grad[:] = ForwardDiff.gradient(loss,par)
    end
    return loss(par)
end

## load all optimization results
cd(joinpath(path_to_trpl,"Optimizations\\Real\\",name_measurement))
opt_outputs = Vector{String}()
for (root, dirs, files) in walkdir(pwd())
    for file in files
		if file[end-3:end] == ".csv"
			push!(opt_outputs,joinpath(root,file))
		end
    end
end

## process optimization parameters and results (needed for output)
n_runs = size(Matrix(CSV.read(opt_outputs[1],DataFrame,skipto=2)),1)-2
n_jobs = length(opt_outputs)
n_opt = Int(n_jobs * n_runs)
println("$n_opt total optimizations: $n_jobs batchjobs with $n_runs optimizations each")
param = Array{Any}(undef,(n_opt,dims))
mle_loss = Array{Any}(undef,(n_opt,1))
rand_seeds = Array{Any}(undef,(n_opt,1))
for i = 0:n_jobs-1
	param[1+(n_runs*i):n_runs+(n_runs*i),:] = Matrix(CSV.read(opt_outputs[i+1],DataFrame,skipto=2))[1:n_runs,3:end-1]
	mle_loss[1+(n_runs*i):n_runs+(n_runs*i)] = Matrix(CSV.read(opt_outputs[i+1],DataFrame,skipto=2))[1:n_runs,2]
	rand_seeds[1+(n_runs*i):n_runs+(n_runs*i)] = Matrix(CSV.read(opt_outputs[i+1],DataFrame,skipto=2))[1:n_runs,end]
end
lb = Matrix(CSV.read(opt_outputs[1],DataFrame,skipto=2))[end-1,3:end-1]
ub = Matrix(CSV.read(opt_outputs[1],DataFrame,skipto=2))[end,3:end-1]
rand_seeds = Int.([parse(BigFloat,rand_seeds[i]) for i = 1:n_opt])
mle_loss = [parse(Float64,mle_loss[i]) for i = 1:n_opt]
opt_results = hcat(mle_loss,param)
sorted_opt_results = opt_results[sortperm(opt_results[:, 1]),:]

## get best fit parameters
p_fit = sorted_opt_results[1,2:end]

## check that the loss function defined here gives the same result
println("sanity check for loss: loss(p_fit) - loss(p_csv) = "*
        "$(loss(p_fit) .- sorted_opt_results[1,1])")

## rescaling parameters for differentiation
p_r = copy(p_fit)
sf = 10 .^ -floor.(log10.(p_r)) |> Diagonal
p_rescaled = sf * p_r

## define rescaled loss function
function loss_rescaled(par)
	_par = inv(sf) * par
    l = 0.0
    function prob_func(prob, i, repeat)
        remake(prob, u0 = [_par[num_pars+i], 0, 0])
    end
    ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
    sol = Array(solve(ensemble_simu, Rosenbrock23(), p = _par[1:num_pars-1], EnsembleThreads(),
		     		  trajectories = num_traj, saveat = t,save_idxs = [1],
		      		  reltol = 1e-6,abstol = 1e-8))[1,:,:]
    for i = 1:num_traj
		model = ((C_scaling * _par[num_pars] / ND[i]) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
		idx_used = idx_IPL[:,i] .& (model .> 0) # we ignore zeros for the log
		l += sum((model[idx_used]) .- (IPL[idx_used,i] .* log.(model[idx_used])))
    end
    return l
end

## define rescaled loss for NLopt
function loss_rescaled_nlopt(par::Vector,grad::Vector)
    if length(grad) > 0
        grad[:] = ForwardDiff.gradient(loss_rescaled,par)
    end
    display(loss_rescaled(par))
    return loss_rescaled(par)
end

## sanity check for rescaled loss function
println("sanity check for loss: loss_rescaled(p_rescaled) - loss(p_csv) = "*
        "$(loss_rescaled(p_rescaled) .- sorted_opt_results[1,1])")

## local optimization without bounds around the found optimum (LBFGS: local optimization)
opt = Opt(:LD_LBFGS,dims)
opt.min_objective = loss_rescaled_nlopt
opt.ftol_abs = 1e-14
opt.maxtime = 600
(optf,p_opt,ret) = optimize(opt,p_rescaled)

## calculate gradient
g = ForwardDiff.gradient(loss_rescaled,p_opt)

## calculate hessian
H = ForwardDiff.hessian(loss_rescaled,p_opt)

## calculate eigenvalues of the hessian
e = eigvals(H)

## calculate standard errors of fit parameters
se = inv(sf) * inv(H) * inv(sf) |> diag .|> sqrt

## generate fit data for plot
IPL_opt = generateipl(inv(sf) * p_opt)

## convert parameters N_st, and N_dt back to ffective DOS
p_fit_output = copy(p_fit)
p_fit_output[7:8] = N_DA ./ p_fit_output[7:8]
p_opt_output = inv(sf) * copy(p_opt)
p_opt_output[7:8] = N_DA ./ p_opt_output[7:8]
se_output = copy(se)
se_output[7:8] = se_output[7:8] .* p_opt_output[7:8] .^2 ./ N_DA

## parameter table
# p_fit is best global optimization result
# p_opt is the final result after the reparameterized local optimization (started at p_fit)
par_table = Array{Any}(undef,7,dims+2)
par_table[1,1] = "p_fit"
par_table[1,2:end] = vcat(loss(p_fit),p_fit_output)
par_table[2,1] = "p_opt"
par_table[2,2:end] = vcat(loss(inv(sf)*p_opt),p_opt_output)
par_table[3,1] = "cihw95"
par_table[3,2:end] = vcat(["/"],q*se_output)
par_table[4,1] = "lb"
par_table[4,2:end] = vcat(["/"],lb)
par_table[5,1] = "ub"
par_table[5,2:end] = vcat(["/"],ub)
par_table[6,1] = "gradient"
par_table[6,2:end] = vcat(["/"],g)
par_table[7,1] = "hessian eigenvalues (sorted)"
par_table[7,2:end] = vcat(["/"],e)

## optimation result table
opt_table = Array{Any}(undef,n_opt,dims+3)
opt_table[:,2:end] = hcat(opt_results,rand_seeds)
opt_table[:,1] = ["Run $i" for i = 1:n_opt]

## data for output
data_table = Array{Any}(undef,size(data_measurement,1)+1,num_traj+1)
data_table[1,1] = "doping"
data_table[1,2:end] = vcat(N_DA,repeat(["/"],num_traj-1))
data_table[2:end,:] = data_measurement
data_table[2,1] = "power"
data_table[3,1] = "ND"
data_table[4,1] = "background"

## variable names
var_names = vcat(["","l_opt", "τ_r0", "τ_nr0", "τ_2t",
				  "τ_1t", "τ_1e", "τ_1d", "N_l", "N_D", "r", "C"],
	 			  [repeat("η_0$i", 1) for i = 1:num_traj])

## generate output data
df_res = DataFrame(par_table,:auto)
rename!(df_res,var_names)
df_fit = DataFrame(hcat(t,IPL_opt),:auto)
rename!(df_fit,vcat("t",["Transient $i" for i = 1:num_traj]))
df_opt = DataFrame(opt_table,:auto)
rename!(df_opt,vcat(var_names,"Seed"))
df_data = DataFrame(data_table,:auto)
rename!(df_data,vcat("/",["Transient $i" for i = 1:num_traj]))

## output
if out
	## output file path
	file_name = "$name_measurement.xlsx"
	xlsx_path = joinpath(path_to_trpl,"Results\\Real",file_name)

	## save output
	if isfile(xlsx_path)
		rm(xlsx_path)
	end
	XLSX.writetable(xlsx_path, Results = (collect(DataFrames.eachcol(df_res)), DataFrames.names(df_res)),
                    Fits = (collect(DataFrames.eachcol(df_fit)), DataFrames.names(df_fit)),
					Opts = (collect(DataFrames.eachcol(df_opt)), DataFrames.names(df_opt)),
					Data = (collect(DataFrames.eachcol(df_data)), DataFrames.names(df_data)))
end

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
using UncertaintyQuantification

σx = [1, 1.1, 0.9, 1.2, 0.8]
σω = [0.7, 1.3, 1.4, 0.6, 0.95]

X = RandomVariable.(Normal.(0, σx), [:X1, :X2, :X3, :X4, :X5])
ω = RandomVariable.(Normal.(0, σω), [:ω1, :ω2, :ω3, :ω4, :ω5])

B = Model(
    df ->
        df.X1 .* df.ω1 + df.X2 .* df.ω2 + df.X3 .* df.ω3 + df.X4 .* df.ω4 + df.X5 .* df.ω5,
    :B,
)

mc = MonteCarlo(10000)

s_mc = sobolindices([B], [X; ω], :B, mc)

# Compare with analytical solution
VB = sum(σx .^ 2 .* σω .^ 2)
analytical = (σx .^ 2 .* σω .^ 2) / VB

# Monte Carlo
@assert all(isapprox.(s_mc.FirstOrder, 0.0, atol=0.1))
@assert all(isapprox.(s_mc.TotalEffect[1:5], analytical, rtol=0.1))
@assert all(isapprox.(s_mc.TotalEffect[6:end], analytical, rtol=0.1))

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
include("splinesParallel.jl")

V=[0.1,0.2,0.3]

B1=fetch(@spawnat 2 larBernsteinBasis(sel(1),4)(V))
B2=fetch(@spawnat 3 larBernsteinBasis(sel(2),4)(V))
B3=fetch(@spawnat 4 larBernsteinBasis(sel(3),4)(V))

C=[B1,B2,B3]


text:
mutable struct TraitLikelihoodXMLElement <: MyXMLElement
    el::XMLOrNothing
    mbd_el::MBDXMLElement
    treeModel_el::TreeModelXMLElement
    extension_el::Union{ModelExtensionXMLElement, Nothing}
    μ::Vector{Float64}
    pss::Float64
    attrs::Dict{String, String}
    xml_name::String

end

function TraitLikelihoodXMLElement(mbd_el::MBDXMLElement,
                    treeModel_el::TreeModelXMLElement,
                    extension_el::Union{ModelExtensionXMLElement, Nothing})
    μ = zeros(size(mbd_el.precision, 1))
    pss = 0.001
    trait_name = treeModel_el.node_traits[1]
    attrs = deepcopy(DEFAULT_TRAITLIKELIHOOD_ATTRS)
    attrs[bn.TRAIT_NAME] = trait_name
    return TraitLikelihoodXMLElement(nothing, mbd_el, treeModel_el,
                extension_el, μ, pss, attrs,
                bn.TRAIT_DATA_LIKELIHOOD)
end

function name(tl::TraitLikelihoodXMLElement)
    return tl.xml_name
end

function get_id(tl::TraitLikelihoodXMLElement)
    return tl.attrs[bn.ID]
end

function copy(x::TraitLikelihoodXMLElement)
    return TraitLikelihoodXMLElement(
        x.el,
        x.mbd_el,
        x.treeModel_el,
        x.extension_el,
        x.μ,
        x.pss,
        copy(x.attrs),
        x.xml_name
    )
end

function get_trait_name(x::TraitLikelihoodXMLElement)
    return x.attrs[bn.TRAIT_NAME]
end

function make_xml(tl_el::TraitLikelihoodXMLElement)
    make_xml(tl_el.mbd_el)
    make_xml(tl_el.treeModel_el)
    make_xml(tl_el.extension_el)

    el = new_element(tl_el.xml_name)
    set_attributes(el, tl_el.attrs)

    add_ref_el(el, tl_el.mbd_el.el)
    add_ref_el(el, tl_el.treeModel_el.el)
    if !isnothing(tl_el.extension_el)
        add_ref_el(el, tl_el.extension_el)
    else
        tp_el = new_child(el, bn.TRAIT_PARAMETER)
        add_parameter_id(tp_el, bn.TRAIT_PARAMETER)
    end
    add_conjugate_root_prior(el, tl_el.μ, tl_el.pss)
    tl_el.el = el
    return el
end

function add_conjugate_root_prior(pel::XMLElement, μ::Vector{Float64}, pss::Float64)
    el = new_element(bn.CONJUGATE_ROOT_PRIOR)
    mean_el = new_element(bn.MEAN_PARAMETER)
    add_parameter(mean_el, value = μ)
    add_child(el, mean_el)

    pss_el = new_element(bn.PSS)
    add_parameter(pss_el, dim = 1, value = [pss])
    add_child(el, pss_el)

    add_child(pel, el)
    return el
end

function no_standardization!(el::TraitLikelihoodXMLElement)
    el.attrs[bn.SCALE_BY_TIME] = bn.FALSE
    el.attrs[bn.STANDARDIZE] = bn.FALSE
end

const DEFAULT_TRAITLIKELIHOOD_ATTRS = Dict(
        [(bn.ID, bn.TRAIT_DATA_LIKELIHOOD_ID),
        (bn.TRAIT_NAME, bn.DEFAULT_TRAIT_NAME),
        (bn.CACHE_BRANCHES, bn.TRUE),
        (bn.ALLOW_IDENTICAL, bn.TRUE),
        (bn.USE_TREE_LENGTH, bn.FALSE),
        (bn.SCALE_BY_TIME, bn.TRUE),
        (bn.REPORT_AS_MULTIVARIATE, bn.TRUE),
        (bn.INTEGRATE_INTERNAL_TRAITS, bn.TRUE),
        (bn.STANDARDIZE, bn.TRUE),
        (bn.ALLOW_SINGULAR, bn.FALSE)]
        )


text:
"""
k = Number of causal SNPs
p = Total number of SNPs
traits = Number of traits (phenotypes)
overlap = number of causal SNPs shared in each trait
"""
function simulate_random_beta(k::Int, p::Int, traits::Int; overlap::Int=0)
    true_b = zeros(p, traits)
    if overlap == 0
        causal_snps = sample(1:(traits * p), k, replace=false)
        true_b[causal_snps] = randn(k)
    else
        shared_snps = sample(1:p, overlap, replace=false)
        weight_vector = aweights(1 / (traits * (p - overlap)) * ones(traits * p))
        for i in 1:traits
            weight_vector[i*shared_snps] .= 0.0 # avoid sampling from shared snps
        end
        @assert sum(weight_vector) ≈ 1.0
        # simulate β for shared predictors
        for i in 1:traits
            true_b[shared_snps, i] = randn(overlap)
        end
        # simulate β for none shared predictors
        nonshared_snps = sample(1:(traits * p), weight_vector, k - traits * overlap, replace=false)
        true_b[nonshared_snps] = randn(k - traits * overlap)
    end

    return true_b
end

@testset "wrapper univariate" begin
    n = 1000  # number of samples
    p = 10000 # number of SNPs
    k = 10    # number of causal SNPs per trait
    for d in [Normal, Bernoulli, Poisson, NegativeBinomial]
        l = d == NegativeBinomial ? LogLink() : canonicallink(d())

        # set random seed for reproducibility
        Random.seed!(2021)

        # simulate `.bed` file with no missing data
        x = simulate_random_snparray("univariate$d.bed", n, p)
        xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)

        # intercept is the only nongenetic covariate
        z = ones(n)
        intercept = 1.0

        # simulate response y, true model b, and the correct non-0 positions of b
        y, true_b, correct_position = simulate_random_response(xla, k, d, l, Zu=z*intercept);

        # create covariate files, `.bim` and `.bam` files, and separate phenotype file
        writedlm("covariates$d.txt", z)
        make_bim_fam_files(x, y, "univariate$d")
        writedlm("univariate$d.phen", y, ',')

        result1 = iht("univariate$d", 11, d, verbose=false, summaryfile="iht.$(d)summary.txt",
            betafile="iht.$(d)beta.txt", max_iter=5)
        result2 = iht("univariate$d", 11, d, covariates="covariates$d.txt", verbose=false,
            summaryfile="iht.$(d)summary.txt", betafile="iht.$(d)beta.txt", max_iter=5)
        result3 = iht("univariate$d", 11, d, covariates="covariates$d.txt", phenotypes="univariate$d.phen",
            summaryfile="iht.$(d)summary.txt", betafile="iht.$(d)beta.txt", verbose=false, max_iter=5)

        @test all(result1.beta .≈ result2.beta .≈ result3.beta)
        @test result1.logl ≈ result2.logl ≈ result3.logl
        @test result1.iter == result2.iter == result3.iter
        @test result1.σg ≈ result2.σg ≈ result3.σg

        Random.seed!(2021)
        result1 = cross_validate("univariate$d", d, verbose=false, cv_summaryfile="cviht.$(d)summary.txt", max_iter=5)
        Random.seed!(2021)
        result2 = cross_validate("univariate$d", d, covariates="covariates$d.txt", verbose=false, max_iter=5)
        Random.seed!(2021)
        result3 = cross_validate("univariate$d", d, covariates="covariates$d.txt",
            phenotypes="univariate$d.phen", verbose=false, max_iter=5)

        @test all(result1 .≈ result2 .≈ result3)
        # [result1 result2 result3]
        # [argmin(result1) argmin(result2) argmin(result3)]

        try
            rm("univariate$d.bim", force=true)
            rm("univariate$d.bed", force=true)
            rm("univariate$d.fam", force=true)
            rm("univariate$d.phen", force=true)
            rm("covariates$d.txt", force=true)
            rm("cviht.$(d)summary.txt", force=true)
            rm("iht.$(d)summary.txt", force=true)
            rm("iht.$(d)beta.txt", force=true)
        catch
            println("Unsuccessful deletion of intermediate files generated in" *
                " unit tests! Windows users can remove these files manually " *
                "located at " * normpath(pathof(MendelIHT) * "../../../test"))
        end
    end
end

@testset "wrapper multivariate" begin
    n = 1000  # number of samples
    p = 10000 # number of SNPs
    k = 10    # number of causal SNPs
    r = 2     # number of traits
    d = MvNormal

    # set random seed for reproducibility
    Random.seed!(111)

    # simulate `.bed` file with no missing data
    x = simulate_random_snparray("multivariate_$(r)traits.bed", n, p)
    xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)

    # intercept is the only nongenetic covariate
    z = ones(n, 1)
    intercepts = [0.5 1.0] # each trait have different intercept

    # simulate β
    B = simulate_random_beta(k, p, r, overlap=2)

    # between trait covariance matrix
    Σ = random_covariance_matrix(r)

    # between sample covariance is identity + GRM (2 times because in SnpArrays grm is halved)
    Φ = 2 .* SnpArrays.grm(x)
    σg = 0.6
    σe = 0.4
    V = σg * Φ + σe * I

    # simulate Y
    μ = z * intercepts + xla * B
    Yt = rand(MatrixNormal(μ', Σ, V))

    # create covariate files (each sample occupies a row)
    writedlm("covariates.txt", z)

    # create `.bim` and `.bam` files using phenotype
    make_bim_fam_files(x, Transpose(Yt), "multivariate_$(r)traits")

    # save phenotypes in separate file (each sample occupies a row)
    writedlm("multivariate_$(r)traits.phen", Yt', ',')

    @test_throws ArgumentError iht("multivariate_$(r)traits", 17, MvNormal)

    result1 = iht("multivariate_$(r)traits", 17, MvNormal, phenotypes=[6, 7], verbose=false)
    result2 = iht("multivariate_$(r)traits", 17, MvNormal, phenotypes=[6, 7], covariates="covariates.txt", verbose=false)
    result3 = iht("multivariate_$(r)traits", 17, MvNormal, phenotypes="multivariate_$(r)traits.phen",
        covariates="covariates.txt", verbose=false);

    true_b1_idx = findall(!iszero, B[:, 1])
    true_b2_idx = findall(!iszero, B[:, 2])
    @test all(result1.beta[1, true_b1_idx] .≈ result2.beta[1, true_b1_idx] .≈ result3.beta[1, true_b1_idx])
    @test all(result1.beta[2, true_b2_idx] .≈ result2.beta[2, true_b2_idx] .≈ result3.beta[2, true_b2_idx])
    @test all(vec(result1.Σ) .≈ vec(result2.Σ) .≈ vec(result3.Σ))
    @test result1.logl ≈ result2.logl ≈ result3.logl
    @test result1.σg ≈ result2.σg ≈ result3.σg
    @test result1.iter == result2.iter == result3.iter

    Random.seed!(2020)
    result1 = cross_validate("multivariate_$(r)traits", d, phenotypes=[6, 7])
    Random.seed!(2020)
    result2 = cross_validate("multivariate_$(r)traits", d, phenotypes=[6, 7], covariates="covariates.txt")
    Random.seed!(2020)
    result3 = cross_validate("multivariate_$(r)traits", d, phenotypes="multivariate_$(r)traits.phen",
        covariates="covariates.txt")

    @test all(result1 .≈ result2 .≈ result3)

    try
        rm("multivariate_$(r)traits.bim", force=true)
        rm("multivariate_$(r)traits.bed", force=true)
        rm("multivariate_$(r)traits.fam", force=true)
        rm("multivariate_$(r)traits.phen", force=true)
        rm("covariates.txt", force=true)
        rm("cviht.summary.txt", force=true)
        rm("iht.summary.txt", force=true)
        rm("iht.beta.txt", force=true)
        rm("iht.cov.txt", force=true)
    catch
        println("Unsuccessful deletion of intermediate files generated in" *
            " unit tests! Windows users can remove these files manually " *
            "located at " * normpath(pathof(MendelIHT) * "../../../test"))
    end
end

@testset "read BGEN and VCF" begin
    cd(normpath(MendelIHT.datadir()))
    xtrue = convert(Matrix{Float64}, SnpArray("normal.bed"), center=true, scale=true, impute=true)
    @time xbgen, _ = MendelIHT.convert_bgen_gt(Float64, Bgen("normal.bgen"))
    @test all(xbgen .≈ xtrue)
    # BGEN wrapper
    @time xbgen, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = MendelIHT.parse_genotypes("normal.bgen")
    @test all(xbgen .≈ xtrue)
    # VCF wrapper
    @time xvcf, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = MendelIHT.parse_genotypes("normal.vcf.gz")
    @test all(xvcf .≈ xtrue)
    # Run IHT on wrapper
    result_plink = iht("normal", 10, Normal, phenotypes="phenotypes.txt", verbose=true)
    result_vcf = iht("normal.vcf.gz", 10, Normal, phenotypes="phenotypes.txt", verbose=true)
    result_bgen = iht("normal.bgen", 10, Normal, phenotypes="phenotypes.txt", verbose=true)
    @test all(result_plink.beta .≈ result_vcf.beta .≈ result_bgen.beta)
    @test result_plink.logl ≈ result_vcf.logl ≈ result_bgen.logl
    @test result_plink.iter == result_vcf.iter == result_bgen.iter
    @test result_plink.σg ≈ result_vcf.σg ≈ result_bgen.σg
    # clean up
    rm("iht.summary.txt", force=true)
    rm("iht.beta.txt", force=true)
end


text:
module Causal

using Spec
using DocStringExtensions
using Cassette
using ..Space
import ..Omega
using ..Omega: RandVar, constant, proj
import ..Omega: apl, ppapl, id, params, ciid
using Callbacks: donothing

include("replace.jl")
include("counterfactual.jl")
include("causes.jl")

export replace,
       cf,
       iscausebf
end


text:
function get_available_components(::Type{T}, sys::PSY.System) where {T <: PSY.Component}
    return PSY.get_components(T, sys, x -> PSY.get_available(x))
end

function get_available_components(
    ::Type{PSY.RegulationDevice{T}},
    sys::PSY.System,
) where {T <: PSY.Component}
    return PSY.get_components(
        PSY.RegulationDevice{T},
        sys,
        x -> (PSY.get_available(x) && PSY.has_service(x, PSY.AGC)),
    )
end

make_system_filename(sys::PSY.System) = "system-$(IS.get_uuid(sys)).json"
make_system_filename(sys_uuid::Union{Base.UUID, AbstractString}) = "system-$(sys_uuid).json"


text:
using DataFrames, JSON

f = open("./full_mcm.kpp")
string = replace(readstring(f),r"\h","")
close(f)

mcm3 = DataFrame(readtable("mcm3_3_1_species.csv",separator=',',header=true))

reactions =[vcat(split(i[1],":")) for i in eachmatch(r"\}(.*);",string)]
groups  =  [match(r"RO2NO\d?|RO2HO2|RO2|NO\d?|HO2|PAN|J",i[2]) for i in reactions]
allocations = [(r"\b(OH)\b","OH"),(r"\b(CH3O|HOCH2CH2O)\b","ALKOXYL"),(r"\b(HCOCO)\b","GLYOXL"),(r"\b(PPN)\b","PAN")]#[[regex,group]]

for i in 1:length(groups)
    groups[i] == nothing ? nothing : continue;
    rxn = split(reactions[i][1],'=')[1]
    groups[i] = match(r"\b(NO\d?|HO2)\b",rxn)

    if groups[i] == nothing

        for aloc in allocations
            if (ismatch( aloc[1] , rxn))
                groups[i] = aloc[2]
                @goto skip
            end
        end

        for j in split(rxn,"+")
            try
                if (ismatch(r"(\[O\-\]\[O\+\]\=C|\[C\]O\[O\]|\[O\]O\[C\])",mcm3[mcm3[:name].==j,:][:smiles][1]) )
                    groups[i] = "CRI"
                    @goto skip
                end
            catch
                groups[i] = "INORGANIC"
                @goto skip
            end
        end

        print(reactions[i][1],groups[i],"\n")
        groups[i]="DISSOCIATION"

        @label skip

    end
end

function rmrx(i)
    try
        return i.match
    end
    return i
end

groups = [rmrx(i) for i in groups] #remove RegexMatch objects

#
print("Groups are: \n", length(Set(groups)))
println(Set(groups))

#########
#convert into graph network links
dict = Dict()

for i in 1:length(groups)
    dict[reactions[i][1]] = groups[i]
end

jsondict = JSON.json(dict)

open("edgecat.json", "w") do f
        write(f, jsondict)
     end

println("\n'fini")

# arrays[1:end-3]
#=
m=match(r"(?P<hour>\d+):(?P<minute>\d+)","12:45")
RegexMatch("12:45", hour="12", minute="45")
julia> m[:minute]
=#

dict

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
int_rules_3_2_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 3.2.3*u*log(e*(f*(a+b*x)^p*(c+d*x)^q)^r)^s =#
    @apply_utils Antiderivative(~(u') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => Antiderivative(~u * log(~e * ((((~b) ^ ~p * ~f) / (~d) ^ ~p) * (~c + ~d * ~x) ^ (~p + ~q)) ^ ~r) ^ ~s, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r, ~s], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => (((~a + ~b * ~x) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s) / ~b - ~r * ~s * (~p + ~q) * Antiderivative(log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1), ~x)) + ((~q * ~r * ~s * (~b * ~c - ~a * ~d)) / ~b) * Antiderivative(log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r, ~s], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ(~p + ~q, 0) && (IGtQ(~s, 0) && LtQ(~s, 4))))
    @apply_utils Antiderivative(log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) / (~(g') + ~(h') * ~x), ~x) => ((log(~g + ~h * ~x) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / ~h - ((~b * ~p * ~r) / ~h) * Antiderivative(log(~g + ~h * ~x) / (~a + ~b * ~x), ~x)) - ((~d * ~q * ~r) / ~h) * Antiderivative(log(~g + ~h * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q, ~r], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')), ~x) => (((~g + ~h * ~x) ^ (~m + 1) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / (~h * (~m + 1)) - ((~b * ~p * ~r) / (~h * (~m + 1))) * Antiderivative((~g + ~h * ~x) ^ (~m + 1) / (~a + ~b * ~x), ~x)) - ((~d * ~q * ~r) / (~h * (~m + 1))) * Antiderivative((~g + ~h * ~x) ^ (~m + 1) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~p, ~q, ~r], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative(log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ 2 / (~(g') + ~(h') * ~x), ~x) => Antiderivative((log((~a + ~b * ~x) ^ (~p * ~r)) + log((~c + ~d * ~x) ^ (~q * ~r))) ^ 2 / (~g + ~h * ~x), ~x) + ((log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) - log((~a + ~b * ~x) ^ (~p * ~r))) - log((~c + ~d * ~x) ^ (~q * ~r))) * (2 * Antiderivative(log((~c + ~d * ~x) ^ (~q * ~r)) / (~g + ~h * ~x), ~x) + Antiderivative(((log((~a + ~b * ~x) ^ (~p * ~r)) - log((~c + ~d * ~x) ^ (~q * ~r))) + log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / (~g + ~h * ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q, ~r], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~b * ~g - ~a * ~h, 0))
    #= Antiderivative(log((~e')*((~f')*((~a')+(~b')*(~x))^(~p')*((~c')+(~d')*(~x))^(~q'))^(~r'))^2/((~g')+(~h') *(~x)),~x) :=  Antiderivative((log((a+b*x)^(p*r))+log((c+d*x)^(q*r)))^2/(g+h*x),x) + (log(e*(f*(a+b*x)^p*(c+d*x)^q)^r)-log((a+b*x)^(p*r))-log((c+d*x)^(q* r)))* Antiderivative((log((a+b*x)^(p*r))+log((c+d*x)^(q*r))+log(e*(f*(a+b*x)^p*(c+ d*x)^q)^r))/(g+h*x),x) <-- FreeQ([a,b,c,d,e,f,g,h,p,q,r],x) &&  NeQ(b*c-a*d,0) && EqQ(b*g-a*h,0) =#
    @apply_utils Antiderivative(log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ 2 / (~(g') + ~(h') * ~x), ~x) => ((log(~g + ~h * ~x) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ 2) / ~h - ((2 * ~b * ~p * ~r) / ~h) * Antiderivative((log(~g + ~h * ~x) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / (~a + ~b * ~x), ~x)) - ((2 * ~d * ~q * ~r) / ~h) * Antiderivative((log(~g + ~h * ~x) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q, ~r], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~s, ~x) => (((~g + ~h * ~x) ^ (~m + 1) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s) / (~h * (~m + 1)) - ((~b * ~p * ~r * ~s) / (~h * (~m + 1))) * Antiderivative(((~g + ~h * ~x) ^ (~m + 1) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1)) / (~a + ~b * ~x), ~x)) - ((~d * ~q * ~r * ~s) / (~h * (~m + 1))) * Antiderivative(((~g + ~h * ~x) ^ (~m + 1) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~p, ~q, ~r, ~s], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IGtQ(~s, 0) && NeQ(~m, -1)))
    @apply_utils Antiderivative(((~(s') + ~(t') * log(~(i') * (~(g') + ~(h') * ~x) ^ ~(n'))) ^ ~(m') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r'))) / (~(j') + ~(k') * ~x), ~x) => (((~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n)) ^ (~m + 1) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) / (~k * ~n * ~t * (~m + 1)) - ((~b * ~p * ~r) / (~k * ~n * ~t * (~m + 1))) * Antiderivative((~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n)) ^ (~m + 1) / (~a + ~b * ~x), ~x)) - ((~d * ~q * ~r) / (~k * ~n * ~t * (~m + 1))) * Antiderivative((~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n)) ^ (~m + 1) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~j, ~k, ~s, ~t, ~m, ~n, ~p, ~q, ~r], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~h * ~j - ~g * ~k, 0) && IGtQ(~m, 0)))
    @apply_utils Antiderivative(((~(s') + ~(t') * log(~(i') * (~(g') + ~(h') * ~x) ^ ~(n'))) * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r'))) / (~(j') + ~(k') * ~x), ~x) => ((log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) - log((~a + ~b * ~x) ^ (~p * ~r))) - log((~c + ~d * ~x) ^ (~q * ~r))) * Antiderivative((~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n)) / (~j + ~k * ~x), ~x) + Antiderivative((log((~a + ~b * ~x) ^ (~p * ~r)) * (~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n))) / (~j + ~k * ~x), ~x) + Antiderivative((log((~c + ~d * ~x) ^ (~q * ~r)) * (~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n))) / (~j + ~k * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~j, ~k, ~s, ~t, ~n, ~p, ~q, ~r], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(((~(s') + ~(t') * log(~(i') * (~(g') + ~(h') * ~x) ^ ~(n'))) ^ ~(m') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(u')) / (~(j') + ~(k') * ~x), ~x) => Unintegrable(((~s + ~t * log(~i * (~g + ~h * ~x) ^ ~n)) ^ ~m * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~u) / (~j + ~k * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~j, ~k, ~s, ~t, ~m, ~n, ~p, ~q, ~r, ~u], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(~u * log(~v) * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => With([g = Simplify(((~v - 1) * (~c + ~d * ~x)) / (~a + ~b * ~x)), h = Simplify(~u * (~a + ~b * ~x) * (~c + ~d * ~x))], (-h * PolyLog(2, 1 - ~v) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s) / (~b * ~c - ~a * ~d) + h * ~p * ~r * ~s * Antiderivative((PolyLog(2, 1 - ~v) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1)) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([g, h], ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r, ~s], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IGtQ(~s, 0) && EqQ(~p + ~q, 0)))
    @apply_utils Antiderivative(~v * log(~(i') * (~(j') * (~(g') + ~(h') * ~x) ^ ~(t')) ^ ~(u')) * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => With([k = Simplify(~v * (~a + ~b * ~x) * (~c + ~d * ~x))], (k * log(~i * (~j * (~g + ~h * ~x) ^ ~t) ^ ~u) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s + 1)) / (~p * ~r * (~s + 1) * (~b * ~c - ~a * ~d)) - ((k * ~h * ~t * ~u) / (~p * ~r * (~s + 1) * (~b * ~c - ~a * ~d))) * Antiderivative(log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s + 1) / (~g + ~h * ~x), ~x) <-- FreeQ(k, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~j, ~p, ~q, ~r, ~s, ~t, ~u], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~p + ~q, 0) && NeQ(~s, -1)))
    @apply_utils Antiderivative(~u * PolyLog(~n, ~v) * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => With([g = Simplify((~v * (~c + ~d * ~x)) / (~a + ~b * ~x)), h = Simplify(~u * (~a + ~b * ~x) * (~c + ~d * ~x))], (h * PolyLog(~n + 1, ~v) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s) / (~b * ~c - ~a * ~d) - h * ~p * ~r * ~s * Antiderivative((PolyLog(~n + 1, ~v) * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ (~s - 1)) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([g, h], ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p, ~q, ~r, ~s], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IGtQ(~s, 0) && EqQ(~p + ~q, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * log((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~(n') / (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2), ~x) => ((2 * ~e * ~g) / (~C * (~e * ~f - ~d * ~g))) * Subst(Antiderivative((~a + ~b * log(~c * ~x)) ^ ~n / ~x, ~x), ~x, sqrt(~d + ~e * ~x) / sqrt(~f + ~g * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~C, ~n], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && EqQ(~B * ~e * ~g - ~C * (~e * ~f + ~d * ~g), 0))
    @apply_utils Antiderivative((~(a') + ~(b') * log((~(c') * sqrt(~(d') + ~(e') * ~x)) / sqrt(~(f') + ~(g') * ~x))) ^ ~(n') / (~(A') + ~(C') * (~x) ^ 2), ~x) => (~g / (~C * ~f)) * Subst(Antiderivative((~a + ~b * log(~c * ~x)) ^ ~n / ~x, ~x), ~x, sqrt(~d + ~e * ~x) / sqrt(~f + ~g * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~C, ~n], ~x) && (EqQ(~C * ~d * ~f - ~A * ~e * ~g, 0) && EqQ(~e * ~f + ~d * ~g, 0))
    @apply_utils Antiderivative(~(RFx') * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')), ~x) => (~p * ~r * Antiderivative(~RFx * log(~a + ~b * ~x), ~x) + ~q * ~r * Antiderivative(~RFx * log(~c + ~d * ~x), ~x)) - ((~p * ~r * log(~a + ~b * ~x) + ~q * ~r * log(~c + ~d * ~x)) - log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r)) * Antiderivative(~RFx, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r], ~x) && (RationalFunctionQ(~RFx, ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && Not(MatchQ(~RFx, ~(u') * (~a + ~b * ~x) ^ ~(m') * (~c + ~d * ~x) ^ ~(n') <-- IntegersQ(m, n)))))
    #= Antiderivative((~RFx)*log((~e')*((~f')*((~a')+(~b')*(~x))^(~p')*((~c')+(~d')*(~x))^(~q'))^(~r')),(~x)*Symbol) := With([u=IntHide(RFx,x)], u*log(e*(f*(a+b*x)^p*(c+d*x)^q)^r) - b*p*r*Antiderivative(u/(a+b*x),x) -  d*q*r*Antiderivative(u/(c+d*x),x) <-- NonsumQ(u)) <--  FreeQ([a,b,c,d,e,f,p,q,r],x) && RationalFunctionQ(RFx,x) &&  NeQ(b*c-a*d,0) =#
    @apply_utils Antiderivative(~RFx * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => With([u = ExpandIntegrand(log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s, ~RFx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r, ~s], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~s, 0))
    @apply_utils Antiderivative(~RFx * log(~(e') * (~(f') * (~(a') + ~(b') * ~x) ^ ~(p') * (~(c') + ~(d') * ~x) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => Unintegrable(~RFx * log(~e * (~f * (~a + ~b * ~x) ^ ~p * (~c + ~d * ~x) ^ ~q) ^ ~r) ^ ~s, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q, ~r, ~s], ~x) && RationalFunctionQ(~RFx, ~x)
    @apply_utils Antiderivative(~(u') * log(~(e') * (~(f') * (~v) ^ ~(p') * (~w) ^ ~(q')) ^ ~(r')) ^ ~(s'), ~x) => Antiderivative(~u * log(~e * (~f * ExpandToSum(~v, ~x) ^ ~p * ExpandToSum(~w, ~x) ^ ~q) ^ ~r) ^ ~s, ~x) <-- FreeQ([~e, ~f, ~p, ~q, ~r, ~s], ~x) && (LinearQ([~v, ~w], ~x) && (Not(LinearMatchQ([~v, ~w], ~x)) && AlgebraicFunctionQ(~u, ~x)))
    @apply_utils Antiderivative(~(u') * log(~(e') * (~(f') * (~g + ~(v') / ~w)) ^ ~(r')) ^ ~(s'), ~x) => Antiderivative(~u * log(~e * ((~f * ExpandToSum(~v + ~g * ~w, ~x)) / ExpandToSum(~w, ~x)) ^ ~r) ^ ~s, ~x) <-- FreeQ([~e, ~f, ~g, ~r, ~s], ~x) && (LinearQ(~w, ~x) && ((FreeQ(~v, ~x) || LinearQ(~v, ~x)) && AlgebraicFunctionQ(~u, ~x)))
    #= Antiderivative(log((~g')*((~h')*((~a')+(~b')*(~x))^(~p'))^(~q'))*log((~i')*((~j')*((~c')+(~d')*(~x))^ (~r'))^(~s'))/((~e)+(~f')*(~x)),~x) :=  1/f*Subst(Antiderivative(log(g*(h*Simp(-(b*e-a*f)/f+b*x/f,x)^p)^q)*log(i*(j*Simp( -(d*e-c*f)/f+d*x/f,x)^r)^s)/x,x),x,e+f*x) <--  FreeQ([a,b,c,d,e,f,g,h,i,j,p,q,r,s],x) =#
end

max_stars_repo_path: src/rules/3 Logarithms/3.2.3 u log(e (f (a+b x)^p (c+d x)^q)^r)^s.jl


text:
using GeneralizedSampling
using Base.Test

#=
Various small functions are defined to interact with a Freq2Wave type.
=#

# ------------------------------------------------------------
# Creating matrices

J = 3
M = 2^(J+1)
N = 2^J

# Uniform samples, No Boundary
Usamples = grid(M, 0.5)
TUNB = Freq2Wave(Usamples, "haar", J)

# Uniform samples, with Boundary
wavename = "db3"
TUB = Freq2Wave(Usamples, wavename, J)

# Non-Uniform samples, No Boundary
K = N/2
NUsamples = N*rand(M) - K
TNUNB = Freq2Wave(NUsamples, "haar", J, K)

# Non-Uniform samples, with Boundary
TNUB = Freq2Wave(NUsamples, wavename, J, K)

# ------------------------------------------------------------
# Testing

@test isuniform(TUNB)
@test isuniform(TUB)
@test !isuniform(TNUNB)
@test !isuniform(TNUB)

@test !hasboundary(TUNB)
@test hasboundary(TUB)
@test !hasboundary(TNUNB)
@test hasboundary(TNUB)

@test dim(TUNB) == dim(TUB) == dim(TNUNB) == dim(TNUB) == 1

@test wscale(TUNB) == wscale(TUB) == wscale(TNUNB) == wscale(TNUB) == J

@test wsize(TUNB) == wsize(TUB) == wsize(TNUNB) == wsize(TNUB) == (N,)

@test size(TUNB) == size(TUB) == size(TNUNB) == size(TNUB) == (M,N)

@test van_moment(TUNB) == van_moment(TNUNB) == 1
@test van_moment(TUB) == van_moment(TNUB) == van_moment(wavename)


text:
using NoiseFourier
using Test

@testset "NoiseFourier.jl" begin
    N = 1024
    @test [NoiseFourier.unidimensional_index(i, N) for i in 0:N] == [j for j in 1:N+1]
    @test [NoiseFourier.unidimensional_index(i, N) for i in -N:-1] == [j for j in N+2:2*N+1]

    for i in 1:2049
        @test NoiseFourier.unidimensional_index(NoiseFourier.inverse_unidim_index(i, N), N) == i
    end
end


text:
module DataSet
    using DataStructures

    mutable struct ColumnSet
        ###parameters
        #arguments
        data::Array
        #hidden parameters
        _data_type::String
        _parent_name::String
        _col_num::Int
        _col_name::String

        ###functions
        #descriptive statistics/common
        describe::Function
        count::Function
        #descriptive statistics/numeric
        sum::Function
        avg::Function
        std::Function
        min::Function
        max::Function
        #descriptive statistics/string
        #top::Function
        unique::Function
        #descriptive statistics/invalid
        _invalid::Function

        #init
        function ColumnSet(data = [], col_num = 0, col_name = "dat", parent = nothing)
            instance = new(data) ; instance._parent_name, instance._col_num, instance._col_name = col_name = parent.name, col_num, col_name
            type_count = Dict("Numeric"=>0, "String"=>0)
            for i in eachindex(instance.data)
                item = string(instance.data[i])
                if !isnull(tryparse(Float64, item))
                    if contains(item, ".")
                        parse_type = Float64
                    else
                        parse_type = Int
                    end

                    instance.data[i] = parse(parse_type, item)
                    type_count["Numeric"] += 1
                else
                    type_count["String"] += 1
                end
            end

            #descriptive statistics/invalid
            instance._invalid = function()
                println(string("invalid for data-type is ", instance._data_type))
                return false
            end

            #set default invalid to functions
            instance.sum, instance.avg, instance.std, instance.min, instance.max, instance.unique, instance.describe = [instance._invalid for i in 1:7]
            if type_count["Numeric"] >= type_count["String"]
                instance._data_type = "Numeric"

                #descriptive statistics/numeric
                parsed_data = [parse(Float64, string(item)) for item in instance.data]
                instance.sum = function()
                    return Base.sum(instance.data)
                end

                instance.avg = function()
                    return Base.mean(parsed_data)
                end

                instance.std = function()
                    return Base.std(parsed_data)
                end

                instance.max = function()
                    Base.maximum(instance.data)
                end

                instance.min = function()
                    Base.minimum(instance.data)
                end
            else
                instance._data_type = "String"
                #instance.data = collect(string(item) for item in instance.data)

                #descriptive statistics/string
                #instance.top = function()
                #    Base.top(instance.data)
                #end

                instance.unique = function()
                    return Base.unique(instance.data)
                end
            end

            #descriptive statistics/common
            instance.describe = function(;return_value = true, header_name = true)
                describe_dict = Dict()
                if instance._data_type == "Numeric"
                    describe_dict["count"] = instance.count()
                    describe_dict["min"] = instance.min()
                    describe_dict["max"] = instance.max()
                    describe_dict["sum"] = instance.sum()
                    describe_dict["std"] = instance.std()
                    describe_dict["avg"] = instance.avg()
                else
                    describe_dict["unique"] = instance.unique()
                end

                describe_string = join(["$key : $value" for (key, value) in describe_dict], "\n")
                if header_name
                    describe_string = string(instance._col_name, "@", instance._parent_name, "\n", describe_string)
                else
                    describe_string = string("ColumnSet-", instance._col_num, "@", instance._parent_name, "\n", describe_string)
                end
                if return_value
                    println(describe_string)
                    return describe_dict
                else
                    return describe_string
                end
            end

            instance.count = function()
                return length(instance.data)
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::ColumnSet)
        println(io, string(instance.describe(return_value = false, header_name = false), "\n\n"))
    end

    mutable struct RowSet
        ###parameters
        #arguments
        data::DataStructures.OrderedDict
        #hidden parameters
        _parent_name::String
        _row_num::Int

        ###functions
        null::Function
        keys::Function

        #eval
        where::Function
        equation::Function

        #init
        function RowSet(data = DataStructures.OrderedDict(), row_num = 0, parent = nothing)
            instance = new(data) ; instance._parent_name, instance._row_num = parent.name, row_num
            instance.null = function()
                result = []
                for (key, value) in instance.data
                    if value == nothing
                        Base.append!(result, key)
                    end
                end

                return result
            end

            instance.keys = function()
                result = Base.keys(instance.data)
                return result
            end

            #eval
            instance.where = function(where_string; drop_null::Bool = true)
                result = nothing
                key_list = instance.keys()
                eval(parse(string(join(key_list, ", "), " = $(values(instance.data))")))
                if eval(parse(where_string))
                    result = collect(values(instance.data))
                else
                    if !drop_null
                        result = [nothing for key in key_list]
                    end
                end

                return result
            end

            instance.equation = function(equation_string)
                result = nothing
                key_list = in(strip(split(equation_string, "=")[1]), instance.keys()) ? instance.keys() : append!(nstance.keys(), strip(split(equation_string, "=")[1]))
                eval(parse(string(join(instance.keys(), ", "), " = $(values(instance.data))")))

                #try
                eval(parse(equation_string))
                result = [eval(parse(key)) for key in key_list]
                #catch
                #    result = [nothing for key in key_list]
                #end

                return result
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::RowSet)
        print_string = string("RowSet-", instance._row_num, "@", instance._parent_name, "\n")
        print_string = string(
            print_string,
            string("keys : ", instance.keys(), "\n")
        )

        println(io, print_string[1:end - 1])
    end
end


text:
using Documenter
using MixedModelsPermutations

makedocs(
    sitename = "MixedModelsPermutations",
    doctest = true,
    pages = [
        "index.md",
    ],
)

deploydocs(repo = "github.com/palday/MixedModelsPermutations.jl.git", push_preview=true, devbranch="main")


text:
using JLD
function CreatePreprocessItems()
    traintest = load("dsge.jld", "traintest")
    thetas = traintest[:,1:9]
    Zs = traintest[:,10:end]
    mZs = mean(Zs,1)
    sZs = std(Zs,1)
    Zs = (Zs .- mZs) ./ sZs
    # get least squares fit
    x = [ones(size(Zs,1)) Zs]
    beta = x\thetas
    # least squares errors
    errors = thetas - x*beta
    # for standardization
    sErrors = std(errors,1)
    return mZs, sZs, beta, sErrors
end

function Preprocess(data)
    mZs, sZs, beta, sErrors = CreatePreprocessItems()
    Y_orig = data[:,1:9]
    Zs = data[:,10:end]
    Zs = (Zs .- mZs) ./ sZs
    x = [ones(size(Zs,1)) Zs]
    # least squares errors
    errors = Y_orig - x*beta
    # for standardization
    errors = errors ./ sErrors
    Y = errors'
    X = Zs'
    return Y,X, Y_orig
end

# this function takes the net's inputs and outputs
# and generates the fit for the original unprocessed
# outputs
function PostProcess(inputs, outputs)
    mZs, sZs, beta, sErrors = CreatePreprocessItems()
    x = [ones(size(inputs,1)) inputs]
    fit = x*beta + sErrors.*outputs # back to original location and scale
    # keep the fits inside support
    fit[:,6] = fit[:,6].*(fit[:,6].> 0.0)
    fit[:,8] = fit[:,8].*(fit[:,8].> 0.0)
    return fit
end
traintest = load("dsge.jld", "traintest")
Y,X = Preprocess(traintest)
trainsize = 500000
YT = Y[:,trainsize+1:end]
XT = X[:,trainsize+1:end]
Y = Y[:,1:trainsize]
X = X[:,1:trainsize]
montecarlo = load("dsge.jld", "montecarlo")
YMC, XMC, Y_orig = Preprocess(montecarlo)


text:
function not_disappearing(data)
    idx=[]
    for s in 1:data.S
        flag = true
        for t in 1:data.T-1
            if data.counts[s,1,t]==0 && data.counts[s,1,t+1]>0
                flag = false
                break
            end
        end

        if flag
            push!(idx,s)
        end
    end
    return idx
end

function filter_counts(data, rounds::Vector{Int}, counts_threshold::Int)
    @assert maximum(rounds) <= data.T
    @assert length(rounds) <= data.T

    data_filtered = deepcopy(data)
    for t in rounds
        data_filtered = subdata(data, findall(data.counts[:,1,t] .>= counts_threshold))
    end
    return data_filtered
end

function add_pseudocounts(data, pc)
    data_new = deepcopy(data)
    data_new.counts .+= pc
    return data_new
end

export not_disappearing, filter_counts, add_pseudocounts


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

using Base.Broadcast

# custom broadcasting style for TSeries
struct MVTSeriesStyle{F<:Frequency} <: Broadcast.BroadcastStyle end
@inline frequencyof(::Type{<:MVTSeriesStyle{F}}) where {F<:Frequency} = F

@inline Base.Broadcast.BroadcastStyle(::Type{<:MVTSeries{F}}) where {F<:Frequency} = MVTSeriesStyle{F}()
@inline Base.Broadcast.BroadcastStyle(S1::MVTSeriesStyle, S2::MVTSeriesStyle) = mixed_freq_error(S1, S2)
@inline Base.Broadcast.BroadcastStyle(S::MVTSeriesStyle{F}, ::MVTSeriesStyle{F}) where {F<:Frequency} = S

# mixing with TSeries - we treat is as a vector column
@inline Base.Broadcast.BroadcastStyle(S1::MVTSeriesStyle, S2::TSeriesStyle) = mixed_freq_error(S1, S2)
@inline Base.Broadcast.BroadcastStyle(S::MVTSeriesStyle{F}, ::TSeriesStyle{F}) where {F<:Frequency} = S

# mixing with other things - we take over as long as we have 2 dimensons or less
@inline Base.Broadcast.BroadcastStyle(S::MVTSeriesStyle, AS::Broadcast.AbstractArrayStyle{N}) where {N} = N <= 2 ? S : AS

function Base.similar(bc::Broadcast.Broadcasted{<:MVTSeriesStyle}, ::Type{ElType}) where {ElType}
    if !isa(bc.axes, Nothing) && length(bc.axes) == 2
        similar(Vector, ElType, bc.axes)
    else
        throw(DimensionMismatch("Cannot create a MVTSeries with axes $(bc.axes)."))
    end
end

# disable extrusion
@inline Base.Broadcast.extrude(x::MVTSeries) = x
@inline Base.Broadcast.extrude(x::TSeries) = x

const MVTSeriesIndexType = Tuple{UnitRange{<:MIT},NTuple{N,Symbol}} where {N}

@inline mvts_get_index(x::Number, p::MIT, c::Symbol) = x
@inline mvts_get_index(x::MVTSeries, p::MIT, c::Symbol) = x[p, c]
@inline mvts_get_index(x::TSeries, p::MIT, c::Symbol) = x[p]
@inline mvts_get_index(x::Base.Broadcast.Broadcasted, p::MIT, c::Symbol) = x[p, c]

function Base.Broadcast.getindex(bc::Base.Broadcast.Broadcasted, p::MIT, c::Symbol)
    args = (mvts_get_index(arg, p, c) for arg in bc.args)
    return bc.f(args...)
end

function Base.copyto!(dest::MVTSeries, bc::Broadcast.Broadcasted{Nothing})
    bc′ = Base.Broadcast.preprocess(dest, bc)
    bcrng, bcvars = axes(bc)
    xrng, xvars = axes(dest)
    for (p, c) in Iterators.product(intersect(bcrng, xrng), intersect(bcvars, xvars))
        @inbounds dest[p, c] = bc′[p, c]
    end
    return dest
end

function do_instantiate(bc)
    shape = mvts_combine_axes(bc.args...)
    if bc.axes !== nothing
        shape = mvts_broadcast_shape(bc.axes, shape)
    end
    mvts_instantiate(bc, shape)
end

@inline function Base.Broadcast.instantiate(bc::Base.Broadcast.Broadcasted{S,A}) where {N,S<:Base.Broadcast.BroadcastStyle,A<:Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}}
    do_instantiate(bc)
end

@inline function Base.Broadcast.instantiate(bc::Base.Broadcast.Broadcasted{S,A}) where {N,S<:Base.Broadcast.AbstractArrayStyle{0},A<:Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}}
    do_instantiate(bc)
end

@inline function Base.Broadcast.instantiate(bc::Base.Broadcast.Broadcasted{S,Nothing}) where {S<:MVTSeriesStyle}
    do_instantiate(bc)
end

# combine_axes works recursively
@inline mvts_combine_axes(A) = axes(A)
@inline mvts_combine_axes(A, B) = mvts_broadcast_shape(axes(A), axes(B))
@inline mvts_combine_axes(A, B...) = mvts_broadcast_shape(axes(A), mvts_combine_axes(B...))

@inline mvts_broadcast_shape(shape::Tuple) = shape
@inline mvts_broadcast_shape(::Tuple{}, shape::Tuple{}) = shape
@inline mvts_broadcast_shape(::Tuple{}, shape::Tuple) = shape
@inline mvts_broadcast_shape(shape::Tuple, ::Tuple{}) = shape
@inline mvts_broadcast_shape(shape1::Tuple{A}, shape2::Tuple) where {A} = mvts_broadcast_shape(shape2, shape1)
function mvts_broadcast_shape(shape1::Tuple, shape2::Tuple)
    if length(shape1) > 2 || length(shape2) > 2
        throw(ArgumentError("broadcasting MVTSeries with ndims > 2."))
    end
    return (mit_common_axes(shape1[1], shape2[1]),
        length(shape2) == 1 ? shape1[2] : sym_common_axes(shape1[2], shape2[2]),
    )
end

@inline mit_common_axes(a::UnitRange{<:MIT}, b::Base.OneTo) =
    length(b) == 1 || length(b) == length(a) ? a :
    throw(DimensionMismatch("Cannot broadcast with $a and $b."))

# mit_common_axes is the same as in TSeries
@inline sym_common_axes(a::NTuple{N1,Symbol}, b::NTuple{N2,Symbol}) where {N1,N2} = a == b ? a : tuple(intersect(a, b)...)
@inline sym_common_axes(a::NTuple{N,Symbol}, b::Any) where {N} = length(a) == length(b) && first(b) == 1 ? a : throw(DimensionMismatch("Cannot broadcast with $(a) and $b."))
@inline sym_common_axes(a::Any, b::NTuple{N,Symbol}) where {N} = sym_common_axes(b, a)

function mvts_instantiate(bc::Base.Broadcast.Broadcasted{S}, shape) where {S<:Base.Broadcast.BroadcastStyle}
    args = map(bc.args) do arg
        mvts_check_axes(shape, arg)
    end
    return Base.Broadcast.Broadcasted{S}(bc.f, args, shape)
end

struct MVTSBroadcasted{Axes,BC}
    singleton::Tuple{Bool,Bool}
    shape::Axes
    bc::BC
    function MVTSBroadcasted(shape, bc)
        bcshape = axes(bc)
        singleton = length(bcshape) == 0 ? (true, true) :
                    length(bcshape) == 1 ? (length(bcshape[1]) == 1, true) :
                    (length(bcshape[1]) == 1, length(bcshape[2]) == 1)
        # bc1 = Base.Broadcast.instantiate(bc)
        new{typeof(shape),typeof(bc)}(singleton, shape, bc)
    end
end

@inline Base.eltype(x::MVTSBroadcasted) = Base.Broadcast.combine_eltypes(x.bc.f, x.bc.args)

function mvts_get_index(x::MVTSBroadcasted, p::MIT, c::Symbol)
    ip = x.singleton[1] ? 1 : Int(p - first(x.shape[1]) + 1)
    ic = x.singleton[2] ? 1 : indexin([c], collect(x.shape[2]))[1]
    ind = CartesianIndex(ip, ic)
    x.bc[ind]
end

# nested broadcasts are processed recursively
@inline mvts_check_axes(shape::MVTSeriesIndexType, bc::Base.Broadcast.Broadcasted{<:MVTSeriesStyle}) = mvts_instantiate(bc, shape)
@inline mvts_check_axes(shape::MVTSeriesIndexType, bc::Base.Broadcast.Broadcasted{<:TSeriesStyle}) = mvts_instantiate(bc, shape)
@inline mvts_check_axes(shape::MVTSeriesIndexType, bc::Base.Broadcast.Broadcasted) = MVTSBroadcasted(shape, bc)

# MVTSs are left alone. If axes are wrong it'll error when indexing later.
@inline mvts_check_axes(shape::MVTSeriesIndexType, x::MVTSeries) = x
# TSeries are also left alone
@inline mvts_check_axes(shape::MVTSeriesIndexType, x::TSeries) = x
# Leave numbers alone too
@inline mvts_check_axes(shape::MVTSeriesIndexType, x::Number) = x
# For Vector, we wrap it in a TSeries
@inline mvts_check_axes(shape::MVTSeriesIndexType, x::AbstractVector) = TSeries(shape[1], x)
# For Matrix, we wrap them in an MVTSeries with the same dimensions
@inline mvts_check_axes(shape::MVTSeriesIndexType, x::AbstractMatrix) = MVTSeries(first(shape[1]), shape[2], x)

function Base.axes(bc::Base.Broadcast.Broadcasted{<:MVTSeriesStyle})
    bc.axes === nothing ? mvts_combine_axes(bc.args...) : bc.axes
end

function Base.axes(bc::Base.Broadcast.Broadcasted{<:MVTSeriesStyle}, d::Integer)
    1 <= d <= 2 ? axes(bc)[d] : Base.OneTo(1)
end


text:
import Base: _range, _rangestyle

maybe_sval(val::Val{V}) where {V} = seek_static_val(eltype(V), val)
maybe_sval(val::SVal) = val
maybe_sval(::Nothing) = nothing
maybe_sval(val) = val

# floor(Int, (stop-start)/step)
for (ST,BT) in S2B
    for B in (ST,BT)
    for S in (SVal,BaseAny)
    for E in (ST,BT)
        if B == ST && S == SVal && E == ST
            @eval begin
                _scolon(::Base.Ordered, ::Any, start::$B, step::$S, stop::$E) =
                    StepSRange(start, step, stop)
                _scolon(::Base.Ordered, ::Base.ArithmeticRounds, start::$B, step::$S, stop::$E) =
                    StepSRangeLen(start, step, floor(Int, (stop-start)/step)+SOne)
                _scolon(::Any, ::Any, start::$B, step::$S, stop::$E) =
                    StepSRangeLen(start, step, floor(Int, (stop-start)/step)+SOne)

            end
        else
            @eval begin
                _scolon(::Base.Ordered, ::Any, start::$B, step::$S, stop::$E) =
                    StepMRange(start, step, stop)
                _scolon(::Base.Ordered, ::Base.ArithmeticRounds, start::$B, step::$S, stop::$E) =
                    StepMRangeLen(start, step, floor(Int, (stop-start)/step)+SOne)
                _scolon(::Any, ::Any, start::$B, step::$S, stop::$E) =
                    StepMRangeLen(start, step, floor(Int, (stop-start)/step)+SOne)
            end
        end
    end
    end
    end
end

_range(b::SReal, ::Nothing, ::Nothing, l::BaseInteger) = UnitMRange{eltype(b)}(b, b + l - SOne)
_range(b::BaseReal, ::Nothing, ::Nothing, l::SInteger) = UnitMRange{eltype(b)}(b, b + l - SOne)
_range(b::SReal, ::Nothing, ::Nothing, l::SInteger) = UnitSRange{eltype(b)}(b, b + l - SOne)

Base.OrderStyle(x::Type{<:SReal}) = Base.OrderStyle(eltype(x))
Base.ArithmeticStyle(x::Type{<:SReal}) = Base.ArithmeticStyle(eltype(x))

_rangestyle(::Base.Ordered, ::Base.ArithmeticWraps, a::SVal, step::SVal, len::SInteger) =
    StepSRange{eltype(a)}(a, step, ofeltype(a, a+step*(len-SOne)))
_rangestyle(::Any, ::Any, a::SVal, step::SVal, len::SInteger) =
    StepSRangeLen{eltype(a+SZero*step)}(a, step, len)

_rangestyle(::Base.Ordered, ::Base.ArithmeticWraps, a::SVal, step::Real, len::Integer) =
    StepMRange{eltype(a)}(a, step, ofeltype(eltype(a), a+step*(len-SOne)))
_rangestyle(::Any, ::Any, a::SVal, step::Real, len::Integer) =
    StepMRangeLen{eltype(a+SZero*step)}(a, step, len)

_rangestyle(::Base.Ordered, ::Base.ArithmeticWraps, a::Real, step::SReal, len::Integer) =
    StepMRange{eltype(a)}(a, step, ofeltype(a, a+step*(len-SOne)))
_rangestyle(::Any, ::Any, a::Real, step::SReal, len::Integer) =
    StepMRangeLen{typeof(a+SZero*step)}(a, step, len)

_rangestyle(::Base.Ordered, ::Base.ArithmeticWraps, a::Real, step::Real, len::SInteger) =
    StepMRange{eltype(a)}(a, step, ofeltype(a, a+step*(len-SOne)))
_rangestyle(::Any, ::Any, a::Real, step::Real, len::SInteger) =
    StepMRangeLen{typeof(a+SZero*step)}(a, step, len)

for B in (SNumber,BaseNumber)
for E in (SNumber,BaseNumber)
for L in (SInteger,BaseInteger)
    if B == BaseNumber && E == BaseNumber && L == BaseInteger
        continue
    elseif B == SNumber && E == SNumber && L == SInteger
        @eval begin
            _range(start::$B, ::Nothing, stop::$E, len::$L) = LinSRange{eltype(start)}(start, stop, len)
        end
    else
        @eval begin
            _range(start::$B, ::Nothing, stop::$E, len::$L) = LinMRange{eltype(start)}(start, stop, len)
        end
    end
end
end
end

for (ST,BT) in SI2BI
    for B in (ST,BT)
    for E in (ST,BT)
    for L in (SInteger,BaseInteger)
        B == BT && E == BT && L == BaseInteger && continue
        @eval begin
            _range(start::$B, ::Nothing, stop::$E, len::$L) =
                linspace(float($BT), start, stop, len)
        end
    end
    end
    end
end

for (SF,BF) in SF2BF
    for B in (SF,BF)
    for S in (SF,BF)
    for L in (SInteger,BaseInteger)
        B == BF && S == SF && L == BaseInteger && continue
        @eval begin
            function _range(start::$B, ::Nothing, stop::$S, len::$L)
                len < 2 && return linspace1($BF, start, stop, len)
                if start == stop
                    return srangehp($BF, start, zero($B), SZero, len, SOne)
                end
                # Attempt to find exact rational approximations
                start_n, start_d = rat(start)
                stop_n, stop_d = rat(stop)
                if start_d != 0 && stop_d != 0
                    den = lcm(start_d, stop_d)
                    m = maxintfloat($B, Int)
                    if den != 0 && abs(den*start) <= m && abs(den*stop) <= m
                        start_n = round(Int, den*start)
                        stop_n = round(Int, den*stop)
                        if $BF(start_n/den) == start && $BF(stop_n/den) == stop
                            return linspace($BF, start_n, stop_n, len, den)
                        end
                    end
                end
                linspace(start, stop, len)
            end

            function _range(a::$B, st::$S, ::Nothing, len::$L)
                start_n, start_d = rat(a)
                step_n, step_d = rat(st)
                if start_d != 0 && step_d != 0 &&
                   ofeltype(a, start_n/start_d) == a &&
                   ofeltype(a, step_n/step_d) == st
                    den = lcm(start_d, step_d)
                    m = maxintfloat(eltype(a), Int)
                    if abs(den*a) <= m && abs(den*st) <= m &&
                            rem(den, start_d) == 0 && rem(den, step_d) == 0
                        start_n = round(Int, den*a)
                        step_n = round(Int, den*st)
                        return sfloatrange(eltype(a), start_n, step_n, len, den)
                    end
                end
                srangehp(eltype(a), a, st, SZero, len, SOne)
            end
        end
    end
    end
    end
end

function sfloatrange(::Type{T}, b::Integer, s::Integer, l::Real, d::Integer) where T
    if l < 2 || s == 0
        return srangehp(T, (b, d), (s, d), SZero, SOne, l)
    end
    # index of smallest-magnitude value
    imin = clamp(round(Int, -b/s+SOne), SOne, ofeltype(Int, l))
    # Compute smallest-magnitude element to 2x precision
    ref_n = b+(imin-SOne)*s  # this shouldn't overflow, so don't check
    nb = Base.nbitslen(T, l, imin)
    srangehp(T, (ref_n, d), (s, d), nb, ofeltype(Int, l), imin)
end

function sfloatrange(b::AbstractFloat, s::AbstractFloat, l::Integer, d::AbstractFloat)
    T = promote_type(typeof(b), typeof(s), typeof(d))
    m = maxintfloat(T, Int)
    if abs(b) <= m && abs(s) <= m && abs(d) <= m
        ia, ist, idivisor = round(Int, b), round(Int, s), round(Int, d)
        if ia == b && ist == s && idivisor == d
            # We can return the high-precision range
            return sfloatrange(T, ia, ist, ofeltype(Int, l), idivisor)
        end
    end
    # Fallback (misses the opportunity to set offset different from 1,
    # but otherwise this is still high-precision)
    srangehp(T, (b, d), (s, d), Base.nbitslen(T, l, SOne), ofeltype(Int, l), SOne)
end

for B in (SInteger,BaseInteger)
for E in (SInteger,BaseInteger)
for L in (SInteger,BaseInteger)
    if B == SInteger && E == SInteger && L == SInteger
        eval(:(linspace(::Type{T}, start::$B, stop::$E, len::$L) where T = LinSRange{T}(start, stop, len)))
    else
        eval(:(linspace(::Type{T}, start::$B, stop::$E, len::$L) where T = LinMRange{T}(start, stop, len)))
    end
    # this had to be included to avoid ambiguity
    eval(:(linspace(::Type{T}, start::$B, stop::$E, len::$L) where T<:Union{Float16, Float32, Float64} = linspace(T, start, stop, len, SOne)))
end
end
end

for (S,B) in ((AbstractFloat16, Float16),(AbstractFloat32, Float32),(AbstractFloat64,Float64))
    @eval begin
        @inline function linspace(b::$S, e::$S, l::Integer)
            (isfinite(b) && isfinite(e)) || throw(ArgumentError("start and stop must be finite, got $b and $e"))
            # Find the index that returns the smallest-magnitude element
            Δ, Δfac = l - b, one(l)
            if !isfinite(Δ)   # handle overflow for large endpoints
                Δ, Δfac = e/l - b/l, ofeltype(Int64, l)
            end
            tmin = -(b/Δ)/Δfac            # t such that (1-t)*start + t*stop == 0
            imin = round(Int, tmin*(l-one(l))+one(l))  # index approximately corresponding to t
            if 1 < imin < l
                # The smallest-magnitude element is in the interior
                t = (imin-one(l)/(l-one(l)))
                ref = (one(t)-t)*b + t*e
                step = imin-one(imin) < l-imin ? (ref-b)/(imin-one(imin)) : (e-ref)/(l-imin)
            elseif imin <= 1
                imin = SOne
                ref = b
                step = (Δ/(l-one(l)))*Δfac
            else
                imin = ofeltype(Int64, l)
                ref = e
                step = (Δ/(l-one(l)))*Δfac
            end
            if l == 2 && !isfinite(step)
                # For very large endpoints where step overflows, exploit the
                # split-representation to handle the overflow
                return srangehp($B, b, (-b, e), SZero, SInt(2))
            end
            # 2x calculations to get high precision endpoint matching while also
            # preventing overflow in ref_hi+(i-offset)*step_hi
            m, k = prevfloat(floatmax(b)), max(imin-one(imin), l-imin)
            step_hi_pre = clamp(step, max(-(m+ref)/k, (-m+ref)/k), min((m-ref)/k, (m+ref)/k))
            nb = Base.nbitslen($B, l, imin)
            step_hi = truncbits(step_hi_pre, nb)
            x1_hi, x1_lo = Base.add12((one(imin)-imin)*step_hi, ref)
            x2_hi, x2_lo = Base.add12((l-imin)*step_hi, ref)
            a, c = (b - x1_hi) - x1_lo, (e - x2_hi) - x2_lo
            step_lo = (c - a)/(l - one(l))
            ref_lo = a - (one(imin) - imin)*step_lo
            srangehp($B, (ref, ref_lo), (step_hi, step_lo), SZero, l, imin)
        end
    end
end

# range for rational numbers, start = start_n/den, stop = stop_n/den
# Note this returns a StepRangeLen
function linspace(::Type{T}, b::Integer, e::Integer, l::Integer, d::Integer
                 ) where {T<:Union{Float16, Float32, Float64}}
    l < 2 && return linspace1(T, b/d, e/d, l)
    b == e && return srangehp(T, (b, d), (zero(b), d), SZero, l)
    tmin = -b/(float(e) - float(b))
    imin = round(Int, tmin*(l-one(l))+one(l))
    imin = clamp(imin, SOne, ofeltype(Int, l))
    ref_num = ofeltype(Int128, l-imin) * b + ofeltype(Int128, imin-one(imin)) * e
    ref_denom = ofeltype(Int128, l-one(l)) * d
    srangehp(T, (ref_num, ref_denom), (ofeltype(Int128, e) - ofeltype(Int128, b), ref_denom),
             Base.nbitslen(T, l, imin), ofeltype(Int64, l), imin)
end

# For len < 2
function linspace1(::Type{T}, b::B, e::E, l::Integer) where {B,E,T<:Union{Float16, Float32, Float64}}
    l >= SZero || throw(ArgumentError("srange($(values(b)), stop=$(values(e)), length=$(values(l))): negative length"))
    if l <= 1
        l == 1 && (b == e || throw(ArgumentError("srange($(b), stop=$(e), length=$(l)): endpoints differ")))
        # Ensure that first(r)==start and last(r)==stop even for len==0
        # The output type must be consistent with steprangelen_hp
        if T<:Union{Float32,Float16}
            return StaticStepRangeLen{T}(f64(b), f64(b) - f64(e), l, SOne)
        else
            return StaticStepRangeLen(tpval(b, zero(T)), tpval(b, -e), l, SOne)
        end
    end
    throw(ArgumentError("should only be called for len < 2, got $l"))
end

SITuple = Tuple{<:Integer,<:Integer}

srangehp(::Type{Float64}, b::SITuple, s::SITuple, nb::Integer, l::Integer, f::Integer) =
    StaticStepRangeLen(tpval(Float64, b), tpval(Float64, s, nb), l, f)

srangehp(::Type{T}, b::SITuple, s::SITuple, nb::Integer, l::Integer, f::Integer) where T =
    StaticStepRangeLen{T}(b[1]/b[2], s[1]/s[2], ofeltype(Int64, l), f)

srangehp(::Type{Float64}, b::F_or_FF, s::F_or_FF, nb::Integer, l::Integer, f::Integer) =
    StaticStepRangeLen(tp64(b), Base.twiceprecision(tp64(s), nb), ofeltype(Int64, l), f)

srangehp(::Type{T}, b::F_or_FF, s::F_or_FF,
         nb::Integer, l::Integer, f::Integer) where {T<:Union{Float16,Float32}} =
    StaticStepRangeLen{T}(f64(b), f64(s), ofeltype(Int64, l), f)