text:
using LoopVectorization: forbroadcast

_extract(::Type{StaticInt{N}}) where {N} = N::Int
_extract(_) = nothing

abstract type AbstractStrideStyle{S,N} <: Base.Broadcast.AbstractArrayStyle{N} end
struct LinearStyle{S,N,R} <: AbstractStrideStyle{S,N} end
struct CartesianStyle{S,N} <: AbstractStrideStyle{S,N} end
@generated function Base.BroadcastStyle(
  ::Type{A},
) where {S,D,T,N,C,B,R,A<:AbstractStrideArray{S,D,T,N,C,B,R}}
  t = Expr(:curly, :Tuple)
  for s ∈ S.parameters
    _s = _extract(s)
    _s === nothing ? push!(t.args, Int) : push!(t.args, _s)
  end
  all(D) ? :(LinearStyle{$t,$N,$R}()) : :(CartesianStyle{$S,$N}())
end
# Base.BroadcastStyle(::Type{A}) where {S,T,N,X,SN,XN,A<:AbstractStrideArray{S,T,N,X,SN,XN,true}} = CartesianStyle{S,N}()
# function reverse_simplevec(S, N = length(S))
#     Srev = Expr(:curly, :Tuple)
#     for n ∈ 1:N
#         push!(Srev.args, S.parameters[N + 1 - n])
#     end
#     if N == 1
#         N += 1
#         insert!(Srev.args, 2, 1)
#     end
#     Srev, N
# # end
# @generated function Base.BroadcastStyle(::Type{Adjoint{T,A}}) where {S,T,N,A<:AbstractStrideArray{S,T,N}}
#     Srev, Nrev = reverse_simplevec(S, N)
#     Expr(:call, Expr(:curly, :CartesianStyle, Srev, Nrev))
# end
# @generated function Base.BroadcastStyle(::Type{Transpose{T,A}}) where {S,T,N,A<:AbstractStrideArray{S,T,N}}
#     Srev, Nrev = reverse_simplevec(S, N)
#     Expr(:call, Expr(:curly, :CartesianStyle, Srev, Nrev))
# end

const StrideArrayProduct = Union{
  LoopVectorization.Product{<:AbstractStrideArray},
  LoopVectorization.Product{<:Any,<:AbstractStrideArray},
  LoopVectorization.Product{<:AbstractStrideArray,<:AbstractStrideArray},
  # LoopVectorization.Product{Adjoint{<:Any,<:AbstractFixedSizeArray}},
  # LoopVectorization.Product{Transpose{<:Any,<:AbstractFixedSizeArray}},
  # LoopVectorization.Product{<:Any,Adjoint{<:Any,<:AbstractFixedSizeArray}},
  # LoopVectorization.Product{<:Any,Transpose{<:Any,<:AbstractFixedSizeArray}}
}

@generated function Base.BroadcastStyle(
  ::Type{P},
) where {
  SA,
  A<:AbstractStrideArray{SA},
  SB,
  B<:AbstractStrideArray{SB},
  P<:LoopVectorization.Product{A,B},
}
  t = Expr(:curly, :Tuple)
  M = _extract(SA.parameters[1])
  D = length(SB.parameters)
  D ∈ (1, 2) ||
    throw(ArgumentError("In A*B, B should be a vector or matrix, but ndims(B) == $D."))
  M === nothing ? push!(t.args, :Int) : push!(t.args, M)
  D == 1 && return :(CartesianStyle{$t,1}())
  N = _extract(SB.parameters[2])
  N === nothing ? push!(t.args, :Int) : push!(t.args, N)
  :(CartesianStyle{$t,2}())
end

@generated Base.BroadcastStyle(
  a::CartesianStyle{S,N1},
  ::Base.Broadcast.DefaultArrayStyle{N2},
) where {S,N1,N2} = N2 > N1 ? Base.Broadcast.Unknown() : :a
@generated Base.BroadcastStyle(
  ::LinearStyle{S,N1},
  ::Base.Broadcast.DefaultArrayStyle{N2},
) where {S,N1,N2} = N2 > N1 ? Base.Broadcast.Unknown() : CartesianStyle{S,N1}()
Base.BroadcastStyle(a::CartesianStyle{S,N}, ::AbstractStrideStyle{S,N}) where {S,N} = a
Base.BroadcastStyle(a::LinearStyle{S,N,R}, ::LinearStyle{S,N,R}) where {S,N,R} = a # ranks match
Base.BroadcastStyle(::LinearStyle{S,N}, ::LinearStyle{S,N}) where {S,N} =
  CartesianStyle{S,N}() # ranks don't match
@generated function Base.BroadcastStyle(
  ::AbstractStrideStyle{S1,N1},
  ::AbstractStrideStyle{S2,N2},
) where {S1,S2,N1,N2}
  N2 > N1 && return :(Base.Broadcast.Unknown())
  S = Expr(:curly, :Tuple)
  for n ∈ 1:N2#min(N1,N2)
    _s1 = S1.parameters[n]
    _s2 = S2.parameters[n]
    s1 = (_s1 isa Int ? _s1 : -1)::Int
    s2 = (_s2 isa Int ? _s2 : -1)::Int
    if s1 == s2
      push!(S.args, s1)
    elseif s2 == 1
      push!(S.args, s1)
    elseif s1 == 1
      push!(S.args, s2)
    elseif s2 == -1
      push!(S.args, s1)
    elseif s1 == -1
      push!(S.args, s2)
    else
      throw("Mismatched sizes: $S1, $S2.")
    end
  end
  for n ∈ N2+1:N1
    push!(S.args, S1.parameters[n])
  end
  Expr(:call, Expr(:curly, :CartesianStyle, S, max(N1, N2)))
end

# function Base.similar(
# ::Base.Broadcast.Broadcasted{FS}, ::Type{T}
# ) where {S,T<:Union{VectorizationBase.FloatingTypes,StrideArrays.VectorizationBase.IntTypes,StrideArrays.VectorizationBase.UIntTypes},N,FS<:AbstractStrideStyle{S,N}}
@generated function to_tuple(::Type{S}, s) where {N,S<:Tuple{Vararg{Any,N}}}
  t = Expr(:tuple)
  Sp = S.parameters
  for i = 1:length(Sp)
    l = _extract(Sp[i])
    if l === nothing
      push!(t.args, Expr(:ref, :s, i))
    else
      push!(t.args, StaticInt(l))
    end
  end
  t
end
function Base.similar(
  bc::Base.Broadcast.Broadcasted{FS},
  ::Type{T},
) where {S,T,N,FS<:AbstractStrideStyle{S,N}}
  StrideArray{T}(undef, to_tuple(S, size(bc)))
end
@generated function to_tuple(::Type{S}) where {N,S<:Tuple{Vararg{StaticInt,N}}}
  map(StaticInt, known(S))
end
function Base.similar(
  bc::Base.Broadcast.Broadcasted{FS},
  ::Type{T},
) where {S<:Tuple{Vararg{StaticInt}},T,N,FS<:AbstractStrideStyle{S,N}}
  StrideArray{T}(undef, to_tuple(S))
end
# Base.similar(sp::StackPointer, ::Base.Broadcast.Broadcasted{FS}, ::Type{T}) where {S,T,FS<:AbstractStrideStyle{S,T},N} = PtrArray{S,T}(sp)

function add_single_element_array!(
  ls::LoopVectorization.LoopSet,
  destname::Symbol,
  bcname::Symbol,
  elementbytes,
)
  LoopVectorization.pushprepreamble!(
    ls,
    Expr(:(=), Symbol("##", destname), Expr(:call, :first, bcname)),
  )
  LoopVectorization.add_constant!(ls, destname, elementbytes)
end
function add_fs_array!(
  ls::LoopVectorization.LoopSet,
  destname::Symbol,
  bcname::Symbol,
  loopsyms::Vector{Symbol},
  indexes,
  S,
  R,
  C,
  elementbytes,
)
  ref = Symbol[]
  # aref = LoopVectorization.ArrayReference(bcname, ref)
  vptrbc = LoopVectorization.vptr(bcname)
  LoopVectorization.add_vptr!(ls, bcname, vptrbc, true) #TODO: is this necessary?
  offset = 0
  Rnew = Int[]
  for (i, n) ∈ enumerate(indexes)
    _s = _extract(S.parameters[i])
    s = (_s === nothing ? -1 : _s)::Int
    # r = R[i]
    # (isone(n) & (stride != 1)) && pushfirst!(ref, LoopVectorization.DISCONTIGUOUS)
    if s == 1
      offset += 1
      bco = bcname
      bcname = Symbol(:_, bcname)
      v = Expr(:call, :view, Expr(:call, :parent, bco))
      foreach(_ -> push!(v.args, :(:)), 1:i-offset)
      push!(v.args, :(One()))
      foreach(_ -> push!(v.args, :(:)), i+1:length(indexes))
      LoopVectorization.pushprepreamble!(ls, Expr(:(=), bcname, v))
    else
      push!(Rnew, R[i])
      push!(ref, loopsyms[n])
    end
  end
  if iszero(length(ref))
    return add_single_element_array!(ls, destname, bcname, elementbytes)
  end
  bctemp = Symbol(:_, bcname)
  mref = LoopVectorization.ArrayReferenceMeta(
    LoopVectorization.ArrayReference(bctemp, ref),
    fill(true, length(ref)),
    vptrbc,
  )
  sp = sort_indices!(mref, Rnew, C)
  if sp === nothing
    LoopVectorization.pushprepreamble!(ls, Expr(:(=), bctemp, bcname))
  else
    ssp = Expr(:tuple)
    append!(ssp.args, sp)
    ssp = Expr(:call, Expr(:curly, :StaticInt, ssp))
    LoopVectorization.pushprepreamble!(
      ls,
      Expr(:(=), bctemp, Expr(:call, :permutedims, bcname, ssp)),
    )
  end
  loadop =
    LoopVectorization.add_simple_load!(ls, destname, mref, mref.ref.indices, elementbytes)
  LoopVectorization.doaddref!(ls, loadop)
end

# function add_broadcast_adjoint_array!(
#     ls::LoopVectorization.LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{A}, indices
# ) where {S, T, N, A <: AbstractStrideArray{S,T,N}}
#     if N == 1
#         if _extract(first(S.parameters)) == 1
#             add_single_element_array!(ls, destname, bcname, sizeof(T))
#         else
#             ref = LoopVectorization.ArrayReference(bcname, Symbol[loopsyms[2]])
#             LoopVectorization.add_load!( ls, destname, ref, sizeof(T) )
#         end
#     else
#         add_fs_array!(ls, destname, bcname, loopsyms, indices, S, sizeof(T))
#     end
# end

function sort_indices!(ar, R, C)
  any(i -> R[i-1] ≥ R[i], 2:length(R)) || return nothing
  li = ar.loopedindex
  inds = LoopVectorization.getindices(ar)
  offsets = ar.ref.offsets
  sp = rank_to_sortperm(R)
  lib = copy(li)
  indsb = copy(inds)
  offsetsb = copy(offsets)
  for i ∈ eachindex(li, inds)
    li[i] = lib[sp[i]]
    inds[i] = indsb[sp[i]]
    offsets[i] = offsetsb[sp[i]]
  end
  C > 0 || pushfirst!(inds, LoopVectorization.DISCONTIGUOUS)
  sp
end

function LoopVectorization.add_broadcast!(
  ls::LoopVectorization.LoopSet,
  destname::Symbol,
  bcname::Symbol,
  loopsyms::Vector{Symbol},
  ::Type{A},
  elementbytes::Int = 8,
) where {S,D,T,N,C,B,R,A<:AbstractStrideArray{S,D,T,N,C,B,R}}
  NN = min(N, length(loopsyms))
  op = add_fs_array!(ls, destname, bcname, loopsyms, Base.OneTo(NN), S, R, C, sizeof(T))
  op
end

_tuple_type_len(_) = nothing
function _tuple_type_len(::Type{S}) where {N,S<:Tuple{Vararg{StaticInt,N}}}
  L = 1
  for n = 1:N
    L *= _extract(S.parameters[n])::Int
  end
  L
end

function linear_indexing_broadcast_quote(S, D, T, N, C, B, R, UNROLL, @nospecialize(BC))
  # we have an N dimensional loop.
  # need to construct the LoopSet
  loopsyms = [gensym(:n)]
  ls = LoopVectorization.LoopSet(:StrideArrays)
  (inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads) = UNROLL
  LoopVectorization.set_hw!(ls, rs, rc, cls, l1, l2, l3)
  ls.vector_width = W
  ls.isbroadcast = isbroadcast
  itersym = first(loopsyms)
  L = _tuple_type_len(S)
  if L === nothing
    Lsym = gensym(:L)
    Rsym = gensym(:R)
    LoopVectorization.pushprepreamble!(
      ls,
      Expr(:(=), Lsym, Expr(:call, :static_length, :dest)),
    )
    LoopVectorization.pushprepreamble!(
      ls,
      Expr(:(=), Rsym, Expr(:call, :(:), :(One()), Lsym)),
    )
    LoopVectorization.add_loop!(
      ls,
      LoopVectorization.Loop(itersym, 1, Lsym, 1, Rsym, Lsym),
      itersym,
    )
  else
    LoopVectorization.add_loop!(
      ls,
      LoopVectorization.Loop(itersym, 1, L, 1, Symbol(""), Symbol("")),
      itersym,
    )
  end
  elementbytes = sizeof(T)
  LoopVectorization.add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
  storeop = LoopVectorization.add_simple_store!(
    ls,
    :dest,
    LoopVectorization.ArrayReference(:dest, loopsyms),
    elementbytes,
  )
  LoopVectorization.doaddref!(ls, storeop)
  resize!(ls.loop_order, LoopVectorization.num_loops(ls)) # num_loops may be greater than N, eg Product
  # fallback in case `check_args` fails
  fallback = :(copyto!(
    dest,
    Base.Broadcast.instantiate(
      Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{$N}}(
        bc.f,
        bc.args,
        axes(dest),
      ),
    ),
  ))
  Expr(
    :block,
    Expr(:meta, :inline),
    LoopVectorization.setup_call(
      ls,
      fallback,
      LineNumberNode(0),
      inline,
      false,
      u₁,
      u₂,
      v,
      threads % Int,
      0,
    ),
    :dest,
  )
  # ls
end
function cartesian_indexing_broadcast_quote(S, D, T, N, C, B, R, UNROLL, @nospecialize(BC))
  # we have an N dimensional loop.
  # need to construct the LoopSet
  loopsyms = [gensym(:n) for _ ∈ 1:N]
  ls = LoopVectorization.LoopSet(:StrideArrays)
  (inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, threads) = UNROLL
  LoopVectorization.set_hw!(ls, rs, rc, cls, l1, l2, l3)
  ls.vector_width = W
  ls.isbroadcast = isbroadcast
  destref = LoopVectorization.ArrayReference(:_dest, copy(loopsyms))
  destmref = LoopVectorization.ArrayReferenceMeta(
    destref,
    fill(true, length(LoopVectorization.getindices(destref))),
  )
  sp = sort_indices!(destmref, R, C)
  Sknown = known(S)
  for n ∈ 1:N
    itersym = loopsyms[n]#isnothing(sp) ? n : sp[n]]
    Sₙ = (Sknown[n])::Union{Int,Nothing}
    if Sₙ === nothing
      Sₙsym = gensym(:Sₙ)
      Rₙsym = gensym(:Rₙ)
      LoopVectorization.pushprepreamble!(
        ls,
        Expr(:(=), Sₙsym, Expr(:call, :size, :dest, static(n))),
      )
      LoopVectorization.pushprepreamble!(
        ls,
        Expr(:(=), Rₙsym, Expr(:call, :(:), :(One()), Sₙsym)),
      )
      LoopVectorization.add_loop!(
        ls,
        LoopVectorization.Loop(itersym, 1, Sₙsym, 1, Rₙsym, Sₙsym),
        itersym,
      )
    else#TODO: handle offsets
      LoopVectorization.add_loop!(
        ls,
        LoopVectorization.Loop(itersym, 1, Sₙ::Int, 1, Symbol(""), Symbol("")),
        itersym,
      )
    end
  end
  elementbytes = sizeof(T)
  # destadj = length(X.parameters) > 1 && last(X.parameters)::Int == 1
  # if destadj
  #     destsym = :dest′
  #     LoopVectorization.pushprepreamble!(ls, Expr(:(=), destsym, Expr(:call, Expr(:(.), :LinearAlgebra, QuoteNode(:Transpose)), :dest)))
  # else
  # end
  LoopVectorization.add_broadcast!(ls, :destination, :bc, loopsyms, BC, elementbytes)
  if isnothing(sp)
    LoopVectorization.pushprepreamble!(ls, Expr(:(=), :_dest, :dest))
  else
    ssp = Expr(:tuple)
    append!(ssp.args, sp)
    ssp = Expr(:call, Expr(:curly, :StaticInt, ssp))
    LoopVectorization.pushprepreamble!(
      ls,
      Expr(:(=), :_dest, Expr(:call, :permutedims, :dest, ssp)),
    )
  end
  storeop = LoopVectorization.add_simple_store!(ls, :destination, destmref, sizeof(T))
  LoopVectorization.doaddref!(ls, storeop)
  # destref = if destadj
  #     ref = LoopVectorization.ArrayReference(:dest′, reverse(loopsyms))
  # else
  #     if first(X.parameters)::Int != 1
  #         pushfirst!(LoopVectorization.getindices(ref), Symbol("##DISCONTIGUOUSSUBARRAY##"))
  #     end
  #     ref
  # end
  resize!(ls.loop_order, LoopVectorization.num_loops(ls)) # num_loops may be greater than N, eg Product
  # return ls
  # fallback in case `check_args` fails
  fallback = :(copyto!(
    dest,
    Base.Broadcast.instantiate(
      Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{$N}}(
        bc.f,
        bc.args,
        axes(dest),
      ),
    ),
  ))
  Expr(
    :block,
    Expr(:meta, :inline),
    LoopVectorization.setup_call(
      ls,
      fallback,
      LineNumberNode(0),
      inline,
      false,
      u₁,
      u₂,
      v,
      threads % Int,
      0,
    ),
    :dest,
  )
end
# function Base.Broadcast.materialize!(
@generated function _materialize!(
  dest::AbstractStrideArray{S,D,T,N,C,B,R},
  bc::BC,
  ::Val{UNROLL},
) where {S,D,T,N,C,B,R,LS,FS<:LinearStyle{LS,N,R},BC<:Base.Broadcast.Broadcasted{FS},UNROLL}
  if (LS.parameters...,) === known(S)
    linear_indexing_broadcast_quote(S, D, T, N, C, B, R, UNROLL, BC)
  else
    cartesian_indexing_broadcast_quote(S, D, T, N, C, B, R, UNROLL, BC)
  end
end
@generated function _materialize!(
  # function _materialize!(
  dest::AbstractStrideArray{S,D,T,N,C,B,R},
  bc::BC,
  ::Val{UNROLL},
) where {S,D,T,N,C,B,R,BC<:Union{Base.Broadcast.Broadcasted,StrideArrayProduct},UNROLL}
  cartesian_indexing_broadcast_quote(S, D, T, N, C, B, R, UNROLL, BC)
end

@inline function Base.Broadcast.materialize!(
  dest::AbstractStrideArray,
  bc::BC,
) where {BC<:Union{Base.Broadcast.Broadcasted,StrideArrayProduct}}
  _materialize!(
    dest,
    bc,
    LoopVectorization.avx_config_val(
      Val((true, zero(Int8), zero(Int8), zero(Int8), true, -1 % UInt)),
      pick_vector_width(eltype(dest)),
    ),
  )
end

# @generated function Base.Broadcast.materialize!(
#     dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::BC
# ) where {S, T <: Union{Float32,Float64}, N, X, A <: AbstractStrideArray{S,T,N,X}, BC <: Union{Base.Broadcast.Broadcasted,StrideArrayProduct}}
#     # we have an N dimensional loop.
#     # need to construct the LoopSet
#     loopsyms = [gensym(:n) for n ∈ 1:N]
#     ls = LoopVectorization.LoopSet(:StrideArrays)
#     LoopVectorization.pushprepreamble!(ls, Expr(:(=), :dest, Expr(:call, :parent, :dest′)))
#     for (n,itersym) ∈ enumerate(loopsyms)
#         LoopVectorization.add_loop!(ls, LoopVectorization.Loop(itersym, 1, (S.parameters[n])::Int))
#     end
#     elementbytes = sizeof(T)
#     LoopVectorization.add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
#     LoopVectorization.add_simple_store!(ls, :dest, LoopVectorization.ArrayReference(:dest, reverse(loopsyms)), elementbytes)
#     resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
#     q = LoopVectorization.lower(ls)
#     push!(q.args, :dest′)
#     pushfirst!(q.args, Expr(:meta,:inline))
#     q
#     # ls
# end

@inline function Base.Broadcast.materialize!(
  dest::AbstractStrideArray{S,D,T,N,C,B,R,X,O},
  bc::Base.Broadcast.Broadcasted{
    Base.Broadcast.DefaultArrayStyle{0},
    Nothing,
    typeof(identity),
    Tuple{T2},
  },
) where {S,D,T,N,C,B,R,X,O,T2<:Number}
  fill!(dest, first(bc.args))
end

@inline function Base.Broadcast.materialize(
  bc::Base.Broadcast.Broadcasted{S},
) where {S<:AbstractStrideStyle}
  ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
  Base.Broadcast.materialize!(similar(bc, ElType), bc)
end

@inline LoopVectorization.vmaterialize(
  bc::Base.Broadcast.Broadcasted{<:AbstractStrideStyle},
) = Base.Broadcast.materialize(bc)
@inline LoopVectorization.vmaterialize!(
  dest,
  bc::Base.Broadcast.Broadcasted{<:AbstractStrideStyle},
) = Base.Broadcast.materialize!(dest, bc)

@inline LoopVectorization.vmaterialize(bc::StrideArrayProduct) =
  Base.Broadcast.materialize(bc)
@inline LoopVectorization.vmaterialize!(dest, bc::StrideArrayProduct) =
  Base.Broadcast.materialize!(dest, bc)

@inline Base.:(+)(A::AbstractStrideArray, B::AbstractStrideArray) = A .+ B
@inline Base.:(-)(A::AbstractStrideArray, B::AbstractStrideArray) = A .- B

@inline Base.unaliascopy(A::AbstractStrideArray) = A


text:
using Rankers.Spectator

region = "kr"

#=
featured = Spectator.featured(region)
@info sort(featured.gameList, by = x -> x.gameId)[1].participants

version = Spectator.version(region)
@info :version version
=#

gameId = 4036891599
meta = Spectator.getGameMetaData(region, gameId)
@info :meta meta

#=
gameId = 4036891599
stats = Spectator.endOfGameStats(region, gameId)
@info :stats stats
=#

ext: jl


text:
#!/usr/bin/env julia

### images.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Commentary:
#
# This file produces images used in documentation.
#
### Code:

using StableRNGs

rng = StableRNG(42)

using Measurements, LombScargle, Plots

sz = (1200, 800)

# First periodograms
ntimes = 1001
t = range(0.01, stop=10pi, length=ntimes)
t += step(t)*rand(rng, ntimes)
s = @. sinpi(t) + 1.5cospi(2t) + rand(rng)
p = lombscargle(t, s)
plot(freqpower(p)..., size = sz, xlim = (0.0, 2.0), xlabel = "Frequency",
     ylabel = "Lomb–Scargle power")
savefig(joinpath(@__DIR__, "freq-periodogram.png"))
plot(periodpower(p)..., size = sz, xlim = (0.5, 2.5), xlabel = "Period",
     ylabel = "Lomb–Scargle power")
savefig(joinpath(@__DIR__, "period-periodogram.png"))

# signal with uncertainties
ntimes = 1001
t = range(0.01, stop=10pi, length=ntimes)
s = @. sinpi(2t)
errors = rand(rng, 0.1:1e-3:4.0, ntimes)
p = lombscargle(t, s, errors, maximum_frequency=1.5)
plot(freqpower(p)..., size = sz, xlim = (0.25, 1.5), xlabel = "Frequency",
     ylabel = "Lomb–Scargle power")
savefig(joinpath(@__DIR__, "freq-uncertainties.png"))
plot(periodpower(p)..., size = sz, xlim = (0.75, 2.0), xlabel = "Period",
     ylabel = "Lomb–Scargle power")
savefig(joinpath(@__DIR__, "period-uncertainties.png"))


text:
module GraphLaplacians

using ChainRulesCore: @non_differentiable
using SparseArrays: SparseMatrixCSC
using LinearAlgebra
using Graphs

include("adjmat.jl")
include("graphs.jl")

@non_differentiable adjacency_matrix(x...)
@non_differentiable degrees(x...)
@non_differentiable normalized_laplacian(x...)

end


text:
"""
Embedds all produced media in one big html file
"""
function generate_preview(path=joinpath(@__DIR__, "preview.html"); media_root=basedir("recorded"))
    open(path, "w") do io
        for file in readdir(media_root)
            media = joinpath(media_root, file)
            println(io, "
 $file
")
            if isdir(media)
                medias = joinpath.(media, readdir(media))
                embed_media(io, medias)
            else
                embed_media(io, media)
            end
        end
    end
end

function get_medias(path, result=[])
    if isdir(path)
        foreach(x-> get_medias(x, result), readdir(path, join=true))
    else
        file, ext = splitext(path)
        if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG", ".gif", ".pdf", ".svg")
            push!(result, load(path))
        elseif ext == ".mp4"
            append!(result, get_frames(path))
        else
            error("Unknown media extension: $ext with path: $path")
        end
    end
    return result
end

using ImageShow

function generate_test_summary(path, recorded_root, refimages_root, scores)
    open(path, "w") do io
        scores_sorted = sort!(collect(scores), by=last, rev=true)
        for (filename, score) in scores_sorted
            media_ref = get_medias(joinpath(refimages_root, filename))
            media_recorded = get_medias(joinpath(recorded_root, filename))
            println(io, "
 $filename : $(round(score, digits=4))
")
            for (ref, rec) in zip(media_ref, media_recorded)
                diff = (ref .- rec) ./ 2.0

                diff2 = (ref .* 0.5) .+ (rec .* 0.5)
                ctx = IOContext(io, :thumbnailsize=> size(diff), :thumbnail => false)
                show(ctx, "text/html", [diff, diff2])
            end
        end
    end
end

function generate_test_summary(path, recorded_root)
    open(path, "w") do io
        for filename in readdir(recorded_root)
            media_recorded = joinpath(recorded_root, filename)
            println(io, "
 $filename
")
            println(io, """

                    $(embed_media(media_recorded))

            """)
        end
    end
end

function tourl(path)
    if Sys.iswindows()
        # There might be a nicer way?
        # Anyways, this seems to be needed on windows
        if !startswith(path, "http")
            path = "file:///" * replace(path, "\\" => "/")
        end
    end
    return repr(path)
end

"""
    embed_image(path::AbstractString)
Returns the html to embed an image
"""
function embed_image(path::AbstractString, alt = "")
    if splitext(path)[2] == "pdf"
        return """

        """
    end
    """

    """
end

"""
    embed_video(path::AbstractString)

Generates a html formatted string for embedding video into Documenter Markdown files
(since `Documenter.jl` doesn't support directly embedding mp4's using ![]() syntax).
"""
function embed_video(path::AbstractString, alt = "")
    """
    <video controls autoplay loop muted>
      <source src=$(tourl(path)) type="video/mp4">
      Your browser does not support mp4. Please use a modern browser like Chrome or Firefox.
    </video>
    """
end

"""
Embeds the most common media types as html
"""
function embed_media(path::String, alt = "")
    file, ext = splitext(path)
    if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG", ".gif", ".pdf", ".svg")
        return embed_image(path, alt)
    elseif ext == ".mp4"
        return embed_video(path, alt)
    elseif isdir(path)
        return sprint(io-> embed_media(io, readdir(path)))
    else
        error("Unknown media extension: $ext with path: $path")
    end
end

embed_media(io::IO, path) = println(io, embed_media(path))

"""
Embeds a vector of media files as HTML
"""
function embed_media(io::IO, paths::AbstractVector{<: AbstractString}, caption = "")
    for path in paths
        println(io, """

                $(embed_media(path, caption))

        """)
    end
end


text:
mutable struct ProjectsPanel <: GtkScrolledWindow

    handle::Ptr{Gtk.GObject}
    list::GtkTreeStore
    treeview::GtkTreeView
    add_button::GtkToolButton
    main_window::MainWindow

    function ProjectsPanel(main_window::MainWindow)

        tv,list,cols = projectspanel_treeview()

        vbox = GtkBox(:v)
        add_button = GtkToolButton("gtk-add")
        rm_button = GtkToolButton("gtk-remove")
        bbox = GtkButtonBox(false)
        push!(bbox,add_button,rm_button)

        sc = GtkScrolledWindow()
        set_gtk_property!(sc,:vexpand,true)
        set_gtk_property!(sc,:hscrollbar_policy,1)
        push!(sc,tv)

        push!(vbox,bbox,sc)
        set_gtk_property!(bbox,:layout_style, Gtk.GConstants.GtkButtonBoxStyle.GTK_BUTTONBOX_START)

        t = new(vbox.handle,list,tv,add_button,main_window)

        signal_connect(projectspanel_treeview_clicked_cb,tv, "button-press-event",
                       Cint, (Ptr{Gtk.GdkEvent},), false,t)
        signal_connect(projectspanel_treeview_keypress_cb,tv, "key-press-event",
                       Cint, (Ptr{Gtk.GdkEvent},), false,t)

        signal_connect(projectspanel_add_button_clicked_cb, add_button, "clicked", Nothing, (), false, t)
        signal_connect(projectspanel_rm_button_clicked_cb, rm_button, "clicked", Nothing, (), false, t)

        Gtk.gobject_move_ref(t,vbox)
    end
end

function projectspanel_treeview()

    list = GtkTreeStore(String, Int)
    tv = GtkTreeView(GtkTreeModel(list))

    cell = GtkCellRendererText()
    cols = GtkTreeViewColumn("Name", cell, Dict([("text",0),("weight",1)]))#the second collumn controls the font weight
    push!(tv,cols)

    return tv,list,cols
end

function update!(w::ProjectsPanel)

    n = readdir(joinpath(HOMEDIR,"config","projects"))
    n = filter(x -> extension(x) == ".json", n)
    n = map(x->splitext(x)[1], n)

    empty!(w.list)
    for i = 1:length(n)
        weight = n[i] == w.main_window.project.name ? 800 : 400 #bold the select project
        push!(w.list,(n[i],weight))
    end

end

function add_project(w::ProjectsPanel)
    ok, name = input_dialog("Project name","project",(("Cancel",0),("Ok",1)),w.main_window)
    ok == 0 && return

    w.main_window.project.name = name
    save(w.main_window.project)
    save(w.main_window)
end

function remove_project(w::ProjectsPanel)

    v = selected(w.treeview, w.list)
    v == nothing && return

    name = v[1]
    name == "default" && return #can't delete default

    rm(joinpath(HOMEDIR,"config","projects","$(name).json"))

    if name == project.name #when current, load default
        load(project,"default", w.main_window.editor, w.main_window; dosave=false)#we shouldn't resave the project we just deleted
    end
end

@guarded nothing function projectspanel_add_button_clicked_cb(widgetptr::Ptr, user_data)
    w = user_data
    add_project(w)
    update!(w)
    return nothing
end

@guarded nothing function projectspanel_rm_button_clicked_cb(widgetptr::Ptr, user_data)
    w = user_data
    remove_project(w)
    update!(w)
    return nothing
end

function on_path_change(w::ProjectsPanel, path)
end

function load(w::ProjectsPanel,treeview::GtkTreeView,list::GtkTreeStore)

    v = selected(treeview, list)
    v == nothing && return

    main_window = w.main_window
    editor = main_window.editor

    load(project,v[1], editor, main_window)
    update!(w)
end

function load(project::Project,name::String, editor::Editor, main_window::MainWindow; dosave=true)

    dosave && save(main_window.project)

    project.name = name
    load(project)
    cd(project.path)
    empty!(editor)
    load_tabs(editor,project)
    save(main_window)
    on_path_change(main_window)
end

@guarded (PROPAGATE) function projectspanel_treeview_clicked_cb(widgetptr::Ptr, eventptr::Ptr, projectspanel)

    treeview = convert(GtkTreeView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    list = projectspanel.list

    if event.button == 3

    else
        if event.event_type == Gtk.GdkEventType.DOUBLE_BUTTON_PRESS
            load(projectspanel,treeview,list)
        end
    end
    return PROPAGATE
end

@guarded (PROPAGATE) function projectspanel_treeview_keypress_cb(widgetptr::Ptr, eventptr::Ptr, projectspanel)
    treeview = convert(GtkTreeView, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    list = projectspanel.list
    if event.keyval == GdkKeySyms.Return
        load(projectspanel,treeview,list)
    end
    return PROPAGATE
end


text:
using StructuralEquationModels, Test, FiniteDiff
import LinearAlgebra: diagind
# import StructuralEquationModels as SEM
include("helper.jl")

dat = example_data("holzinger_swineford")
dat_missing = example_data("holzinger_swineford_missing")
solution_lav = example_data("holzinger_swineford_solution")

dat_g1 = dat[dat.school .== "Pasteur", :]
dat_g2 = dat[dat.school .== "Grant-White", :]

dat_miss_g1 = dat_missing[dat_missing.school .== "Pasteur", :]
dat_miss_g2 = dat_missing[dat_missing.school .== "Grant-White", :]

############################################################################
### specification - RAMMatrices
############################################################################

x = Symbol.(:x, 1:36)
#x = [x...]

F = zeros(9, 12)
F[diagind(F)] .= 1.0

A = Matrix{Any}(zeros(12, 12))
A[1, 10] = 1.0; A[4, 11] = 1.0; A[7, 12] = 1.0
A[2:3, 10] .= [x...][16:17]; A[5:6, 11] .= [x...][18:19]; A[8:9, 12] .= [x...][20:21];

# group 1
S1 = Matrix{Any}(zeros(12, 12))
S1[diagind(S1)] .= [x...][1:12]
S1[10, 11] = x[13]; S1[11, 10] = x[13]
S1[10, 12] = x[14]; S1[12, 10] = x[14]
S1[12, 11] = x[15]; S1[11, 12] = x[15]

# group 2
S2 = Matrix{Any}(zeros(12, 12))
S2[diagind(S2)] .= [x...][22:33]
S2[10, 11] = x[34]; S2[11, 10] = x[34]
S2[10, 12] = x[35]; S2[12, 10] = x[35]
S2[12, 11] = x[36]; S2[11, 12] = x[36]

specification_g1 = RAMMatrices(;
    A = A,
    S = S1,
    F = F,
    parameters = x,
    colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed])

specification_g2 = RAMMatrices(;
    A = A,
    S = S2,
    F = F,
    parameters = x,
    colnames = [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :visual, :textual, :speed])

partable = EnsembleParameterTable(
    specification_g1,
    specification_g2;
    groups = [:Pasteur, :Grant_White]
    )

specification_miss_g1 = nothing
specification_miss_g2 = nothing

start_test = [fill(1.0, 9); fill(0.05, 3); fill(0.01, 3); fill(0.5, 6); fill(1.0, 9); fill(0.05, 3); fill(0.01, 3)]
semdiff = SemDiffOptim
@testset "RAMMatrices | constructor | Optim" begin include("build_models.jl") end

############################################################################
### specification - Graph
############################################################################

using StenoGraphs

# w.o. meanstructure -------------------------------------------------------

latent_vars = [:visual, :textual, :speed]
observed_vars = Symbol.(:x, 1:9)

graph = @StenoGraph begin
    # measurement model
    visual  → fixed(1.0, 1.0)*x1 + label(:λ₂, :λ₂)*x2 + label(:λ₃, :λ₃)*x3
    textual → fixed(1.0, 1.0)*x4 + label(:λ₅, :λ₅)*x5 + label(:λ₆, :λ₆)*x6
    speed   → fixed(1.0, 1.0)*x7 + label(:λ₈, :λ₈)*x8 + label(:λ₉, :λ₉)*x9
    # variances and covariances
    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars)   ⇔ _(latent_vars)
end

partable = EnsembleParameterTable(;
    graph = graph,
    observed_vars = observed_vars,
    latent_vars = latent_vars,
    groups = [:Pasteur, :Grant_White])

specification = RAMMatrices(partable)

specification_g1 = specification[:Pasteur]
specification_g2 = specification[:Grant_White]

# w. meanstructure (fiml) ----------------------------------------------------

latent_vars = [:visual, :textual, :speed]
observed_vars = Symbol.(:x, 1:9)

graph = @StenoGraph begin
    # measurement model
    visual  → fixed(1.0, 1.0)*x1 + label(:λ₂, :λ₂)*x2 + label(:λ₃, :λ₃)*x3
    textual → fixed(1.0, 1.0)*x4 + label(:λ₅, :λ₅)*x5 + label(:λ₆, :λ₆)*x6
    speed   → fixed(1.0, 1.0)*x7 + label(:λ₈, :λ₈)*x8 + label(:λ₉, :λ₉)*x9
    # variances and covariances
    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars)   ⇔ _(latent_vars)

    Symbol("1") → _(observed_vars)
end

partable_miss = EnsembleParameterTable(;
    graph = graph,
    observed_vars = observed_vars,
    latent_vars = latent_vars,
    groups = [:Pasteur, :Grant_White])

specification_miss = RAMMatrices(partable_miss)

specification_miss_g1 = specification_miss[:Pasteur]
specification_miss_g2 = specification_miss[:Grant_White]

start_test = [fill(0.5, 6); fill(1.0, 9); 0.05; 0.01; 0.01; 0.05; 0.01; 0.05; fill(1.0, 9); 0.05; 0.01; 0.01; 0.05; 0.01; 0.05]
semdiff = SemDiffOptim
@testset "Graph → Partable → RAMMatrices | constructor | Optim" begin include("build_models.jl") end


text:
using Documenter, SpringerQuarantineBooks

DocMeta.setdocmeta!(SpringerQuarantineBooks,
                    :DocTestSetup,
                    :(using SpringerQuarantineBooks;),
                    recursive = true)
makedocs(sitename = "SpringerQuarantineBooks",
         modules = [SpringerQuarantineBooks],
         pages = [
             "Home" => "index.md",
             "API" => "api.md"
             ]
         )

deploydocs(repo = "github.com/Nosferican/SpringerQuarantineBooks.jl.git", push_preview = true)


text:
using Plots

function calc_midpoints(edges::AbstractVector)::Vector{Float64}
    return Float64[0.5 * (edges[i] + edges[i+1]) for i ∈ 1:length(edges)-1]
end

function bin_data(data::AbstractVector, bins::Integer)
    low, high = extrema(data)
    edges = collect(range(low, stop = high, length = bins + 1))
    midpoints = calc_midpoints(edges)
    buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges))) - 1 for x ∈ data]
    counts = zeros(Int, length(midpoints))
    for b in buckets
        counts[b] += 1
    end
    return edges, midpoints, buckets, counts
end

rectangle(x, y, w, h) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

function myhistogram(data; bins = 10, norm = true, alpha = 0.5)
    edges, midpoints, buckets, _counts = bin_data(data, bins)
    if norm
        tot = length(data)
        counts = _counts ./ tot
    else
        counts = _counts
    end
    p = plot()
    for i = 1:length(midpoints)
        emin, emax = edges[i], edges[i+1]
        c = Float64(counts[i])
        r = rectangle(emin, zero(c), emax - emin, c)
        plot!(p, r, legend = false, color = :blue, alpha = alpha)
    end
    p = plot(p)
    return p
end

function main()
    data = randn(1000)
    dmin, dmax = extrema(data)
    bins = 10
    norm = true
    p = myhistogram(data, bins = bins, norm = norm, alpha = 1)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# List samples after alignment. Looks for sampleIDs by findind files called "SAMPLE_ID.bam".
# syn  - only used if path is a synapseID
# path - a local path or a synapseID
# runName - name of the run. Only used if path is a SynapseID, where it assumes runName is a subfolder. Otherwise path is used directly.
# returns paths, names (where paths are local paths or Synapse IDs)
function listaligned(syn, path::AbstractString, runName::AbstractString)
    if syn != nothing && SynapseClient.utils.is_synapse_id(path) != nothing
        folder = getchildbyname(syn, path, runName)
        isempty(folder) && error("Could not find \"$runName\" in \"$path\".")
    else
        folder = path
    end
    paths, names = listfiles(syn, folder)

    # filter by extension
    # mask = map( x->ismatch(r".bam$",x), names )
    mask = match.(r".bam$",names) .!= nothing
    paths[mask], names[mask]
end

# get sample names and references used for all samples in logFile
function referencefromlog(syn, logFile::AbstractString)
    logFileText = readlines(localpath(syn, logFile)) # load entire text file
    map!(rstrip, logFileText, logFileText)

    # find lines of type
    # Assigned reference "RefName" to sample "SampleName".

    pattern = r"^Assigned reference \"([^\"]+)\" to sample \"([^\"]+)\".$"
    matches = map( x->match(pattern,x), logFileText )
    matches = matches[map!(x->x!=nothing, falses(length(matches)), matches)] # remove non-matches

    logRefs = [splitext(m.captures[1])[1] for m in matches] # reference name without file ending (.fasta)
    logSamples = [m.captures[2] for m in matches]           # sample names

    logSamples, logRefs
end

# get references used for given logFile and list of samples
function referencefromlog(syn, logFile::AbstractString, sampleNames::AbstractArray)
    logSamples, logRefs = referencefromlog(syn, logFile)

    references = Vector{String}(undef,length(sampleNames))
    for (i,name) in enumerate(sampleNames)
        ind = findfirst(isequal(name), logSamples)
        references[i] = ind!=0 ? logRefs[ind] : ""
    end
    references
end

# get references used for given logFile and a single sample
function referencefromlog(syn, logFile::AbstractString, sampleNames::AbstractString)
    referencefromlog(syn, logFile, [sampleNames])[1]
end

function wrongreferencefromlog(syn, logFile::AbstractString)
    logFileText = readlines(localpath(syn, logFile)) # load entire text file
    map!(rstrip, logFileText, logFileText)

    # find lines of type
    # WARNING: SampleName is closer to reference "wrongRef.fasta" than to reference "RefName.fasta".

    pattern = r"^WARNING: ([^\"]+) is closer to reference \"([^\"]+)\" than to reference \"([^\"]+)\"\.$"
    matches = map( x->match(pattern,x), logFileText )
    matches = matches[map!(x->x!=nothing, falses(length(matches)), matches)] # remove non-matches

    logSamples = String[m.captures[1] for m in matches]     # sample names
    logWrongRefs = [splitext(m.captures[2])[1] for m in matches] # reference name without file ending (.fasta)
    logOrigRefs  = [splitext(m.captures[3])[1] for m in matches] # reference name without file ending (.fasta)

    logSamples, logWrongRefs, logOrigRefs
end


text:
abstract type NN end

struct Network <: NN

    id::Inf64
    name::string
    obj::Any

    Synapse[i] #synapses#  are always between the axon of one neuron and the sum[Op: operation] of another.
    #Q1. is this the best representation ?
    #Q2.how about a network of meta-conns(you know, with amongst synapses ) ?

    axon #Q3. how does this make the synapse better ? [ isa significant axon asif it hasa `distict features` ]

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
const CARFILE = joinpath(@__DIR__, "..", "icons", "racing_car_top_view.svg")
const PEDFILE = joinpath(@__DIR__, "..", "icons", "walking_person.svg")

function render_fancy_car(
    ctx           :: CairoContext,
    x             :: Real, # x-pos of the center of the vehicle
    y             :: Real, # y-pos of the center of the vehicle
    yaw           :: Real, # heading angle [rad]
    length        :: Real, # vehicle length
    width         :: Real, # vehicle width
    color_fill    :: Colorant
    )

    # # renders a car from the top view sport car svg
    # # (x,y) are in meters and yaw is the radians, counter-clockwise from pos x axis

    save(ctx)
    r = open(CARFILE, "r") do carfile
        car_svg = read(carfile, String)
        car_svg = replace(car_svg, "#CAAC00"=>"#"*hex(color_fill))
        Rsvg.handle_new_from_data(car_svg)
    end
    d = Rsvg.handle_get_dimensions(r)

    # scaling factor
    xdir, ydir = length/d.width, width/d.height

    translate(ctx, x, y)
    scale(ctx, xdir, ydir)
    rotate(ctx, yaw)
    translate(ctx, -d.width/2, -d.height/2)

    Rsvg.handle_render_cairo(ctx, r)

    restore(ctx)

end

function render_fancy_pedestrian(
    ctx           :: CairoContext,
    x             :: Real, # x-pos of the center of the vehicle
    y             :: Real, # y-pos of the center of the vehicle
    yaw           :: Real, # heading angle [rad]
    length        :: Real, # vehicle length
    width         :: Real, # vehicle width
    color_fill    :: Colorant
    )
    # renders a car from the top view sport car svg
    # (x,y) are in meters and yaw is the radians, counter-clockwise from pos x axis

    save(ctx)
    r = open(PEDFILE, "r") do pedfile
        ped_svg = read(pedfile, String)
        ped_svg = replace(ped_svg, "#C90000"=>"#"*hex(color_fill))
        Rsvg.handle_new_from_data(ped_svg)
    end
    d = Rsvg.handle_get_dimensions(r)

    # scaling factor
    xdir, ydir = length/d.width, width/d.height

    translate(ctx, x, y)
    scale(ctx, xdir, ydir)
    rotate(ctx, yaw + pi/2)
    translate(ctx, -d.width/2, -d.height/2)

    Rsvg.handle_render_cairo(ctx, r)

    restore(ctx)
end


text:
""" Predator prey model with logistic function from http://www.math.lsa.umich.edu/~rauch/256/F2Lab5.pdf """
function foxes_and_rabbits(oct::Bool = true)
    m = Model(SOLVER_SILENT)
    t = 20
    r = 0.2
    x1 = 0.6
    y1 = 0.5
    @variable(m, x[1:t] >= 0.01) # Note: Ipopt solution does not converge with an upper bound!!
    @variable(m, dx[1:t-1])
    @variable(m, y[1:t] >= 0.01)
    @variable(m, dy[1:t-1])
    @constraint(m, x[1] == x1)
    @constraint(m, y[1] == y1)
    @constraint(m, [i=2:t], x[i] == x[i-1] + dx[i-1])
    @constraint(m, [i=2:t], y[i] == y[i-1] + dy[i-1])

    # GlobalModel representation
    set_upper_bound.(x, 1)
    set_upper_bound.(dx, 1)
    set_lower_bound.(dx, -1)
    set_upper_bound.(y, 1)
    set_upper_bound.(dy, 1)
    set_lower_bound.(dy, -1)

    if !oct
        @NLconstraint(m, [i=1:t-1], dx[i] == x[i]*(1-x[i]) - x[i]*y[i]/(x[i]+1/5))
        @NLconstraint(m, [i=1:t-1], dy[i] == r*y[i]*(1-y[i]/x[i]))
        return m
    else
        gm = GlobalModel(model = m, name = "foxes_rabbits")
        add_nonlinear_constraint(gm, :((x, y, dx) -> dx[1] - (x[1]*(1-x[1]) -x[1]*y[1]/(x[1]+$(r)))),
                                vars = [x[1], y[1], dx[1]], equality=true)
        add_nonlinear_constraint(gm, :((x, y, dy) -> dy[1] - $(r)*y[1]*(1-y[1]/x[1])),
                                vars = [x[1], y[1], dy[1]], equality=true)
        for i = 2:t-1
            add_linked_constraint(gm, gm.bbls[1], [x[i], y[i], dx[i]])
            add_linked_constraint(gm, gm.bbls[2], [x[i], y[i], dy[i]])
        end
        return gm
    end
end


text:
# ArgMacros
# constants.jl

const USAGE_SYMBOL = Symbol("@helpusage")
const DESCRIPTION_SYMBOL = Symbol("@helpdescription")
const EPILOG_SYMBOL = Symbol("@helpepilog")
const ARGHELP_SYMBOL = Symbol("@arghelp")

const ARGUMENT_REQUIRED_SYMBOL = Symbol("@argumentrequired")
const ARGUMENT_DEFAULT_SYMBOL = Symbol("@argumentdefault")
const ARGUMENT_OPTIONAL_SYMBOL = Symbol("@argumentoptional")
const ARGUMENT_FLAG_SYMBOL = Symbol("@argumentflag")
const ARGUMENT_COUNT_SYMBOL = Symbol("@argumentcount")
const FLAGGED_SYMBOLS = [ARGUMENT_REQUIRED_SYMBOL, ARGUMENT_DEFAULT_SYMBOL, ARGUMENT_OPTIONAL_SYMBOL,
                         ARGUMENT_FLAG_SYMBOL, ARGUMENT_COUNT_SYMBOL]

const POSITIONAL_REQUIRED_SYMBOL = Symbol("@positionalrequired")
const POSITIONAL_DEFAULT_SYMBOL = Symbol("@positionaldefault")
const POSITIONAL_OPTIONAL_SYMBOL = Symbol("@positionaloptional")
const POSITIONAL_OPTIONAL_SYMBOLS = [POSITIONAL_DEFAULT_SYMBOL, POSITIONAL_OPTIONAL_SYMBOL]


text:
module FastTransforms

using AbstractFFTs, DSP, FastGaussQuadrature, FFTW, Libdl, LinearAlgebra, SpecialFunctions, ToeplitzMatrices

import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
             inv, size, view

import Base.GMP: Limb
import Base.MPFR: BigFloat, _BigFloat

import AbstractFFTs: Plan

import DSP: conv

import FastGaussQuadrature: unweightedgausshermite

import FFTW: dct, dct!, idct, idct!,
             plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
             plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct,
			 plan_bfft, plan_bfft!, plan_brfft, fftwNumber

import LinearAlgebra: mul!, lmul!, ldiv!

export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
       lag2lag, jac2ultra, ultra2jac, jac2cheb,
       cheb2jac, ultra2cheb, cheb2ultra,
       sph2fourier, sphv2fourier, disk2cxf, tri2cheb, tet2cheb,
       fourier2sph, fourier2sphv, cxf2disk, cheb2tri, cheb2tet

export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
       plan_lag2lag, plan_jac2ultra, plan_ultra2jac, plan_jac2cheb,
       plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra,
       plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis,
       plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
       plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,
       plan_tri2cheb, plan_tri_synthesis, plan_tri_analysis,
       plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis

include("libfasttransforms.jl")

export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3, plan_inufft1, plan_inufft2
export nufft, nufft1, nufft2, nufft3, inufft1, inufft2

include("nufft.jl")
include("inufft.jl")

export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!,
       paduapoints, plan_paduatransform!, plan_ipaduatransform!

include("PaduaTransform.jl")

export plan_chebyshevtransform, plan_ichebyshevtransform, plan_chebyshevtransform!, plan_ichebyshevtransform!,
            chebyshevtransform, ichebyshevtransform, chebyshevpoints,
            plan_chebyshevutransform, plan_ichebyshevutransform, plan_chebyshevutransform!, plan_ichebyshevutransform!,
            chebyshevutransform, ichebyshevutransform

include("chebyshevtransform.jl")

export plan_clenshawcurtis, clenshawcurtisnodes, clenshawcurtisweights
export plan_fejer1, fejernodes1, fejerweights1,
       plan_fejer2, fejernodes2, fejerweights2

include("clenshawcurtis.jl")
include("fejer.jl")

export weightedhermitetransform, iweightedhermitetransform

include("hermite.jl")

include("fftBigFloat.jl")

export gaunt

include("gaunt.jl")

export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
       sphvones, sphvzeros, sphvrand, sphvrandn,
       diskones, diskzeros, diskrand, diskrandn,
       triones, trizeros, trirand, trirandn, trievaluate,
       tetones, tetzeros, tetrand, tetrandn

lgamma(x) = logabsgamma(x)[1]

include("specialfunctions.jl")

end # module


text:
col_1_to_30 = [ -2.5   0.5           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0       -2.0       1.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
1.0  -3.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0       -2.0       1.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0          -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.5
0.0   0.0           0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        2.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       2.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         2.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   2.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   2.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   2.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     2.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   2.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.5
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        2.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
-0.0  -0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       2.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         2.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   2.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   2.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   2.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     2.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   2.0    0.0
0.0   0.0           1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -6.0        6.66667  -2.66667     0.0   0.0   0.0   0.0     0.0   0.0    1.0
0.0   0.0           0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.666667  -2.0       0.666667    0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.666667  -2.66667  -4.0         6.0  -4.0   2.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0       -14.0  20.0  -8.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0        -2.0  -2.0   2.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         2.0  -8.0  -8.0  18.0   -12.0   6.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0  -3.5     2.5  -1.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0  -0.25   -2.0   0.25   0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.25   -1.0  -2.75   0.0
0.0   6.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0   -6.0
0.0   2.498e-16     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        2.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    1.0
-0.0  -1.66533e-15   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       2.0         0.0   0.0   0.0   0.0     0.0   0.0   -4.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         2.0   0.0   0.0   0.0     0.0   0.0    9.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   2.0   0.0   0.0     0.0   0.0    3.747e-16
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   2.0   0.0     0.0   0.0   -2.498e-15
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   2.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     2.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   2.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0
0.0   0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0       0.0         0.0   0.0   0.0   0.0     0.0   0.0    0.0]

col_31_to_40 = [   0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.5           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.5            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.5            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.5            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.5           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.5     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.5    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.5     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.5           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.5            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.5            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.5            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.5           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.5     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.5    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.5     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
1.0           0.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           1.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            1.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            1.0            0.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            1.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           1.0     0.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     1.0    0.0     0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.0     0.0    1.0     0.0   1.5
-4.0           1.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
-3.0          -1.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
10.0          -9.0            0.0            0.0            0.0           0.0     0.0    0.0     0.0   0.0
-18.0          27.0          -12.0          -12.0            3.0           0.0     0.0    0.0     0.0   0.0
-7.49401e-16   1.1241e-15     3.0           -3.0           -3.0           0.0     0.0    0.0     0.0   0.0
4.996e-15    -7.49401e-15  -12.0           30.0          -21.0           0.0     0.0    0.0     0.0   0.0
0.0           0.0            1.125         -2.25           3.375        -4.125  -1.5    0.375   0.0   0.0
0.0           0.0            4.68375e-17   -9.36751e-17    1.40513e-16   0.375  -3.0   -0.375   0.0   0.0
0.0           0.0           -3.1225e-16     6.245e-16     -9.36751e-16  -1.5     3.75  -5.25    0.0   0.0
0.0           0.0            0.0            0.0            0.0           0.25   -0.5    0.75   -2.5   2.0
0.0           0.0            0.0            0.0            0.0           0.5    -1.0    1.5     1.0  -2.0]

B_data = [col_1_to_30 col_31_to_40]


text:
using Test
using Statistics
using DelimitedFiles

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SurfaceFlow
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, vars_state
using ClimateMachine.ArtifactWrappers

@testset "V Catchment Maxwell River Model" begin
    tv_dataset = ArtifactWrapper(
        @__DIR__,
        isempty(get(ENV, "CI", "")),
        "tiltedv",
        ArtifactFile[ArtifactFile(
            url = "https://caltech.box.com/shared/static/qi2gftjw2vu2j66b0tyfef427xxj3ug7.csv",
            filename = "TiltedVOutput.csv",
        ),],
    )
    tv_dataset_path = get_data_folder(tv_dataset)

    ClimateMachine.init()
    FT = Float64

    soil_water_model = PrescribedWaterModel()
    soil_heat_model = PrescribedTemperatureModel()
    soil_param_functions = nothing

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)

    function x_slope(x, y)
        MFT = eltype(x)
        if x < MFT(800)
            MFT(-0.05)
        elseif x <= MFT(820)
            MFT(0)
        else
            MFT(0.05)
        end
    end

    function y_slope(x, y)
        MFT = eltype(x)
        MFT(-0.02)
    end

    function channel_mannings(x, y)
        MFT = eltype(x)
        return x >= MFT(800) && x <= MFT(820) ? MFT(2.5 * 60 * 10^-3) :
               MFT(2.5 * 60 * 10^-4)
    end

    m_surface = OverlandFlowModel(x_slope, y_slope; mannings = channel_mannings)

    bc = LandDomainBC(
        miny_bc = LandComponentBC(
            surface = Dirichlet((aux, t) -> eltype(aux)(0)),
        ),
        minx_bc = LandComponentBC(
            surface = Dirichlet((aux, t) -> eltype(aux)(0)),
        ),
        maxx_bc = LandComponentBC(
            surface = Dirichlet((aux, t) -> eltype(aux)(0)),
        ),
    )

    function init_land_model!(land, state, aux, localgeo, time)
        state.surface.height = eltype(state)(0)
    end

    # units in m / s
    precip(x, y, t) = t < eltype(t)(90 * 60) ? eltype(t)(3e-6) : eltype(t)(0.0)

    sources = (Precip{FT}(precip),)

    m = LandModel(
        param_set,
        m_soil;
        surface = m_surface,
        boundary_conditions = bc,
        source = sources,
        init_state_prognostic = init_land_model!,
    )

    N_poly_hori = 1
    N_poly_vert = 1
    xres = FT(20)
    yres = FT(20)
    zres = FT(1)
    # Specify the domain boundaries.
    zmax = FT(1)
    zmin = FT(0)
    xmax = FT(1620)
    ymax = FT(1000)

    driver_config = ClimateMachine.MultiColumnLandModel(
        "LandModel",
        (N_poly_hori, N_poly_vert),
        (xres, yres, zres),
        xmax,
        ymax,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = RusanovNumericalFlux(),
    )

    t0 = FT(0)
    timeend = FT(180 * 60)
    dt = FT(0.5)

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    mygrid = solver_config.dg.grid
    Q = solver_config.Q

    h_index = varsindex(vars_state(m, Prognostic(), FT), :surface, :height)
    n_outputs = 60

    every_x_simulation_time = ceil(Int, timeend / n_outputs)

    dons = Dict([k => Dict() for k in 1:n_outputs]...)

    iostep = [1]
    callback = GenericCallbacks.EveryXSimulationTime(
        every_x_simulation_time,
    ) do (init = false)
        t = ODESolvers.gettime(solver_config.solver)
        h = Q[:, h_index, :]
        all_vars = Dict{String, Array}("t" => [t], "h" => h)
        dons[iostep[1]] = all_vars
        iostep[1] += 1
        return
    end

    ClimateMachine.invoke!(solver_config; user_callbacks = (callback,))

    aux = solver_config.dg.state_auxiliary
    x = Array(aux[:, 1, :])
    y = Array(aux[:, 2, :])
    z = Array(aux[:, 3, :])
    # Get points at outlet (y = ymax)
    mask2 = (Float64.(y .== 1000.0)) .== 1
    n_outputs = length(dons)
    function compute_Q(h, xv)
        height = max.(h, 0.0)
        v = calculate_velocity(m_surface, xv, 1000.0, height)# put in y = 1000.0
        speed = sqrt(v[1]^2.0 + v[2]^2.0 + v[3]^2.0)
        Q_outlet = speed .* height .* 60.0 # multiply by 60 so it is per minute, not per second
        return Q_outlet
    end
    # We divide by 4 because we have 4 nodal points with the same value at each x (z = 0, 1)
    # Multiply by xres because the solution at each point roughly represents that the solution for that range in x
    Q =
        [
            sum(compute_Q.(Array(dons[k]["h"])[:][mask2[:]], x[mask2[:]]))
            for k in 1:n_outputs
        ] ./ 4.0 .* xres
    data = joinpath(tv_dataset_path, "TiltedVOutput.csv")
    ds_tv = readdlm(data, ',')
    error = sqrt(mean(Q .- ds_tv))
    @test error < 5e-7
end


text:
"""
    visualize_solution(results::Dict{String, Any}, data::Dict{String, Any}; gate_sequence = false)

Given dictionaries of results and data, and assuming that the optimization model had a feasible solution,
this function aids in visualizing the optimal circuit decomposition.
"""
function visualize_solution(results::Dict{String, Any}, data::Dict{String, Any}; gate_sequence = false)

    if results["primal_status"] != MOI.FEASIBLE_POINT

        if results["termination_status"] != MOI.TIME_LIMIT
            Memento.warn(_LOGGER, "Infeasible primal status. Gate decomposition may be inaccurate")
        else
            Memento.warn(_LOGGER, "Optimizer hits time limit with an infeasible primal status. Gate decomposition may be inaccurate")
        end

        return
    else
        gates_sol, gates_sol_compressed = QCO.get_postprocessed_circuit(results, data)
    end

    if !isempty(gates_sol_compressed)

        printstyled("\n","=============================================================================","\n"; color = :cyan)

        printstyled("Quantum Circuit Model Data","\n"; color = :red)

        printstyled("\n","  ","Number of qubits: ", data["num_qubits"], "\n"; color = :cyan)

        printstyled("  ","Total number of elementary gates (after presolve): ",size(data["gates_real"])[3],"\n"; color = :cyan)

        printstyled("  ","Maximum depth of decomposition: ", data["maximum_depth"],"\n"; color = :cyan)

        printstyled("  ","Input elementary gates: ", data["elementary_gates"],"\n"; color = :cyan)

        if "discretization" in keys(data)
            for i in keys(data["discretization"])
                printstyled("    ","$i discretization: ", ceil.(rad2deg.(data["discretization"][i]), digits = 1),"\n"; color = :cyan)
            end
        end

        # printstyled("  ","Input target gate: ", data["target_gate"]["type"],"\n"; color = :cyan)

        printstyled("  ","Type of decomposition: ", data["decomposition_type"],"\n"; color = :cyan)

        printstyled("  ","MIP optimizer: ", results["optimizer"],"\n"; color = :cyan)

        printstyled("\n","Optimal Circuit Decomposition","\n","\n"; color = :red)

        print("  ")

        for i=1:length(gates_sol_compressed)

            if i != length(gates_sol_compressed)
                printstyled(gates_sol_compressed[i], " * "; color = :cyan)
            else
                if data["decomposition_type"] == "exact"
                    printstyled(gates_sol_compressed[i], " = ", "Target gate","\n"; color = :cyan)
                elseif data["decomposition_type"] == "approximate"
                    printstyled(gates_sol_compressed[i], " ≈ ", "Target gate","\n"; color = :cyan)
                end
            end

        end

        if data["decomposition_type"] == "approximate"
            printstyled("  ","||Decomposition error||₂: ", round(LA.norm(results["solution"]["slack_var"]), digits = 10),"\n"; color = :cyan)
        end

        if data["objective"] == "minimize_depth"

            if length(data["identity_idx"]) >= 1 && !(results["termination_status"] == MOI.TIME_LIMIT)
                printstyled("  ","Minimum optimal depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
            else
                printstyled("  ","Decomposition depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
            end

        elseif data["objective"] == "minimize_cnot"

            if !isempty(data["cnot_idx"])

                if data["decomposition_type"] == "exact"
                    printstyled("  ","Minimum number of CNOT gates: ", round(results["objective"], digits = 6),"\n"; color = :cyan)

                elseif data["decomposition_type"] == "approximate"
                    printstyled("  ","Minimum number of CNOT gates: ", round((results["objective"] - results["objective_slack_penalty"]*LA.norm(results["solution"]["slack_var"])^2), digits = 6),"\n"; color = :cyan)
                end

            elseif !isempty(data["identity_idx"])

                printstyled("  ","Objective function changed to minimize_depth (CNOT gate not found in input)", "\n"; color = :cyan)

                if !(results["termination_status"] == MOI.TIME_LIMIT)
                    printstyled("  ","Minimum optimal depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
                else
                    printstyled("  ","Decomposition depth: ", length(gates_sol_compressed),"\n"; color = :cyan)
                end

            else

                printstyled("  ","Objective function changed to feasiblity (CNOT and Identity gates not found in input)", "\n"; color = :cyan)
                printstyled("  ","Compressed depth for a feasbile decomposition: ", length(gates_sol_compressed),"\n"; color = :cyan)

            end

        end

        printstyled("  ","Optimizer run time: ", ceil(results["solve_time"], digits=2)," sec.","\n"; color = :cyan)

        if results["termination_status"] == MOI.TIME_LIMIT
            printstyled("  ","Termination status: MOI.TIME_LIMIT", "\n"; color = :cyan)
            printstyled("  ","Decomposition may not be optimal", "\n"; color = :cyan)
        end

        printstyled("=============================================================================","\n"; color = :cyan)

    else
        Memento.warn(_LOGGER, "Valid integral feasible solutions could not be found to visualize the solution")
    end

    if gate_sequence
        return gates_sol
    end

end

function get_postprocessed_circuit(results::Dict{String, Any}, data::Dict{String, Any})

    gates_sol = Array{String,1}()
    id_sequence = Array{Int64,1}()

    for d = 1:data["maximum_depth"]
        id = findall(isone.(round.(abs.(results["solution"]["z_onoff_var"][:,d]), digits=3)))[1]
        push!(id_sequence, id)

        gate_id = data["gates_dict"]["$id"]

        if !("Identity" in gate_id["type"])

            s1 = gate_id["type"][1]

            if occursin(kron_symbol, s1)
                push!(gates_sol, s1)

            elseif !(QCO._parse_gate_string(s1, type = true) in union(QCO.ONE_QUBIT_GATES_ANGLE_PARAMETERS, QCO.TWO_QUBIT_GATES_ANGLE_PARAMETERS))
                push!(gates_sol, s1)

            else

                s2 = String[]
                for i_qu = 1:data["num_qubits"]
                    if gate_id["qubit_loc"] == "qubit_$i_qu"
                        s2 = "$i_qu"
                    end
                end

                if "angle" in keys(gate_id)

                    if length(keys(gate_id["angle"])) == 1
                        θ = round(rad2deg(gate_id["angle"]), digits = 3)
                        s3 = "$(θ)"
                        push!(gates_sol, string(s1,"(", s3, ")"))
                    elseif length(keys(gate_id["angle"])) == 3
                        θ = round(rad2deg(gate_id["angle"]["θ"]), digits = 3)
                        ϕ = round(rad2deg(gate_id["angle"]["ϕ"]), digits = 3)
                        λ = round(rad2deg(gate_id["angle"]["λ"]), digits = 3)
                        s3 = string("(","$(θ)",",","$(ϕ)", ",","$(λ)",")")
                        push!(gates_sol, string(s1, s3))
                    end

                end
            end
        end
    end

    QCO.validate_circuit_decomposition(data, id_sequence)

    gates_sol_compressed = QCO.get_depth_compressed_circuit(data["num_qubits"], gates_sol)

    return gates_sol, gates_sol_compressed
end

"""
    validate_circuit_decomposition(data::Dict{String, Any}, id_sequence::Array{Int64,1})

This function validates the circuit decomposition if it is indeed exact with respect to the specified target gate.
"""
function validate_circuit_decomposition(data::Dict{String, Any}, id_sequence::Array{Int64,1})

    M_sol = Array{Complex{Float64},2}(Matrix(LA.I, 2^(data["num_qubits"]), 2^(data["num_qubits"])))

    for i in id_sequence
        M_sol *= data["gates_dict"]["$i"]["matrix"]
    end

    # This tolerance is very important for the final feasiblity check
    if data["are_gates_real"]
        target_gate = real(data["target_gate"])
    else
        target_gate = QCO.real_to_complex_gate(data["target_gate"])
    end

    if (data["decomposition_type"] == "exact") && (!isapprox(M_sol, target_gate, atol = 1E-4))
        Memento.error(_LOGGER, "Decomposition is not valid: Problem may be infeasible")
    end

end

"""
    get_depth_compressed_circuit(num_qubits::Int64, gates_sol::Array{String,1})

Given the number of qubits and the sequence of gates from the solution, this function returns a
decomposition of gates after compressing adjacent pair of gates represented on two separate qubits.
For example, gates H1 and H2 appearing in a sequence will be compressed to H1xH2 (kron(H1,H2)).
This functionality is currently supported only for two qubit circuits and gates without angle parameters.
"""
function get_depth_compressed_circuit(num_qubits::Int64, gates_sol::Array{String,1})
    # This part of the code may be hacky. This needs to be updated once the input format gets cleaned up for elementary gates with U and R gates.

    if (length(gates_sol) == 1) || (num_qubits > 2)
        return gates_sol
    end

    gates_sol_compressed = String[]

    angle_param_gate = false
    for i=1:length(gates_sol)
        if !occursin(kron_symbol, gates_sol[i])
            gates_sol_type = QCO._parse_gate_string(gates_sol[i], type = true)

            if gates_sol_type in union(QCO.ONE_QUBIT_GATES_ANGLE_PARAMETERS, QCO.TWO_QUBIT_GATES_ANGLE_PARAMETERS)
                angle_param_gate = true
                break
            end
        end
    end

    if !angle_param_gate

        status = false

        for i=1:(length(gates_sol))
            if i <= length(gates_sol) - 1
                if status
                    status = false
                    continue
                else
                    gate_i = QCO.is_multi_qubit_gate(gates_sol[i])
                    gate_iplus1 = QCO.is_multi_qubit_gate(gates_sol[i+1])

                    if !(gate_i) && !(gate_iplus1)
                        if (occursin('1', gates_sol[i]) && occursin('2', gates_sol[i+1])) || (occursin('2', gates_sol[i]) && occursin('1', gates_sol[i+1]))
                            if occursin('1', gates_sol[i])
                                gate_string = string(gates_sol[i],"x",gates_sol[i+1])
                            else
                                gate_string = string(gates_sol[i+1],"x",gates_sol[i])
                            end
                            push!(gates_sol_compressed, gate_string)
                            status = true
                            continue
                        else
                            push!(gates_sol_compressed, gates_sol[i])

                        end
                    else
                        push!(gates_sol_compressed, gates_sol[i])

                    end
                end
            else
                if !status
                    push!(gates_sol_compressed, gates_sol[i])
                end
            end
        end

    else
        return gates_sol
    end

    if isempty(gates_sol_compressed)
        Memento.error(_LOGGER, "Compressed gates solution is empty")
    end

    return gates_sol_compressed
end


text:
logistic(x) = 1 / (1 + exp(-x))

function autoencoder_cost(We1, We2, Wd, b1, b2, x)
    firstLayer = logistic.(We1 * x .+ b1)
    encodedInput = logistic.(We2 * firstLayer .+ b2)
    reconstructedInput = logistic.(Wd * encodedInput)
    cost = sum((reconstructedInput .- x) .^ 2.0f0)
    return cost
end

function mlp1(w1, w2, w3, x1)
    xx2 = w1 * x1
    x2 = log.(1.0f0 .+ exp.(xx2))
    xx3 = w2 * x2
    x3 = log.(1.0f0 .+ exp.(xx3))
    x4 = w3 * x3
    sum(1.0f0 ./ (1.0f0 .+ exp.(-x4)))
end


text:
using HPack
using Base.Test
import HttpCommon: Headers

include("codes.jl")
include("musts.jl")


text:
abstract type Guard{T} end

# ******************************************************************************
# T
# ******************************************************************************
struct GuardT{T} <: Guard{T}
    T :: AbstractVector{T}
    W :: AbstractVector{T}
end

function Guard(grid::GridT{T}, field) where T
    (; arch, wu, t, w, tguard, wguard) = grid

    Tguard = guard_window_both(t, tguard)
    Tguard = adapt_array(arch, Tguard)

    Wguard = @. exp(-((w * wu)^2 / wguard^2)^20)
    Wguard = adapt_array(arch, Wguard)

    return GuardT{T}(Tguard, Wguard)
end

function apply_guard_real_domain!(E::AbstractVector, guard::GuardT)
    @. E = E * guard.T
    return nothing
end

function apply_guard_spec_domain!(E::AbstractVector, guard::GuardT)
    @. E = E * guard.W
    return nothing
end

# ******************************************************************************
# RT
# ******************************************************************************
struct GuardRT{T} <: Guard{T}
    R :: AbstractVector{T}
    T :: AbstractVector{T}
    K :: AbstractMatrix{T}
    W :: AbstractVector{T}
end

function Guard(grid::GridRT{T}, field) where T
    (; arch, Nr, Nt, wu, ku, r, t, k, w) = grid
    (; rguard, tguard, kguard, wguard) = grid

    Rguard = guard_window_right(r, rguard)
    Rguard = adapt_array(arch, Rguard)

    Tguard = guard_window_both(t, tguard)
    Tguard = adapt_array(arch, Tguard)

    Kguard = zeros((Nr, Nt))
    for it=1:Nt
        kmax = k_func(w[it] * wu) * sind(kguard)
        if kmax != 0
            for ir=1:Nr
                Kguard[ir, it] = exp(-((k[ir] * ku)^2 / kmax^2)^20)
            end
        end
    end
    Kguard = adapt_array(arch, Kguard)

    Wguard = @. exp(-((w * wu)^2 / wguard^2)^20)
    Wguard = adapt_array(arch, Wguard)

    return GuardRT{T}(Rguard, Tguard, Kguard, Wguard)
end

function apply_guard_real_domain!(E::AbstractMatrix, guard::GuardRT)
    guard_mul!(E, guard.R; dim=1)
    guard_mul!(E, guard.T; dim=2)
    return nothing
end

function apply_guard_spec_domain!(E::AbstractMatrix, guard::GuardRT)
    @. E = E * guard.K
    guard_mul!(E, guard.W; dim=2)
    return nothing
end

# ******************************************************************************
# Tools
# ******************************************************************************
function guard_window_left(x::AbstractArray, width; p::Int=6)
    if width >= (x[end] - x[1])
        error("Guard width is larger or equal than the grid size.")
    end
    N = length(x)
    if width == 0
        guard = ones(N)
    else
        xloc1 = x[1]
        xloc2 = x[1] + width
        gauss1 = zeros(N)
        gauss2 = ones(N)
        for i=1:N
            if x[i] >= xloc1
                gauss1[i] = 1 - exp(-((x[i] - xloc1) / (width / 2))^p)
            end
            if x[i] <= xloc2
                gauss2[i] = exp(-((x[i] - xloc2) / (width / 2))^p)
            end
        end
        guard = @. (gauss1 + gauss2) / 2
    end
    return guard
end

function guard_window_right(x::AbstractArray, width; p::Int=6)
    if width >= (x[end] - x[1])
        error("Guard width is larger or equal than the grid size.")
    end
    N = length(x)
    if width == 0
        guard = ones(N)
    else
        xloc1 = x[end] - width
        xloc2 = x[end]
        gauss1 = ones(N)
        gauss2 = zeros(N)
        for i=1:N
            if x[i] >= xloc1
                gauss1[i] = exp(-((x[i] - xloc1) / (width / 2))^p)
            end
            if x[i] <= xloc2
                gauss2[i] = 1 - exp(-((x[i] - xloc2) / (width / 2))^p)
            end
        end
        guard = @. (gauss1 + gauss2) / 2
    end
    return guard
end

function guard_window_both(x::AbstractArray, width; p::Int=6)
    if width >= (x[end] - x[1]) / 2
        error("Guard width is larger or equal than the grid size.")
    end
    lguard = guard_window_left(x, width; p=p)
    rguard = guard_window_right(x, width; p=p)
    return @. lguard + rguard - 1
end

# ------------------------------------------------------------------------------
function guard_mul!(E::Matrix, A::Vector; dim::Int=1)
    ci = CartesianIndices(size(E))
    for ici=1:length(ci)
        idim = ci[ici][dim]
        E[ici] = E[ici] * A[idim]
    end
    return nothing
end

function guard_mul!(E::CuMatrix, A::CuVector; dim::Int=1)
    N = length(E)
    ckernel = @cuda launch=false guard_mul_kernel(E, A, dim)
    config = launch_configuration(ckernel.fun)
    threads = min(N, config.threads)
    blocks = cld(N, threads)
    ckernel(E, A, dim; threads=threads, blocks=blocks)
    return nothing
end

function guard_mul_kernel(E, A, dim)
    id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    ci = CartesianIndices(size(E))
    for ici=id:stride:length(ci)
        idim = ci[ici][dim]
        E[ici] = E[ici] * A[idim]
    end
    return nothing
end


text:
"""
All stochastic iterated integrals are written in Stratonovich sense as indicated
by the J.
"""

abstract type AbstractJ end

"""
Kloeden, P. E., Platen, E., & Wright, I. W., The approximation of multiple stochastic integrals.
Stochastic analysis and applications, 10(4), pp. 431-441 (1992).
DOI: 10.1080/07362999208809281
"""
abstract type AbstractKPWJ <: AbstractJ end

struct KPWJ_oop <: AbstractKPWJ end

mutable struct KPWJ_iip{WikAType, WikJType, randType} <: AbstractKPWJ
  WikA::WikAType
  WikJ::WikJType
  μ::randType
  ζ::randType
  η::randType
  kronprod1::WikAType
  kronprod2::WikAType
end

function KPWJ_iip(ΔW)
  WikJ = false .* vec(ΔW) .* vec(ΔW)'
  WikA = false .* vec(ΔW) .* vec(ΔW)'

  μ = zero(ΔW)
  ζ = zero(ΔW)
  η = zero(ΔW)

  kronprod1 = zero(WikA)
  kronprod2 = zero(WikA)

  KPWJ_iip(WikA,WikJ,μ,ζ,η,kronprod1,kronprod2)
end

function get_iterated_I(dt, dW, dZ, Wik::KPWJ_oop, p=nothing, C=1, γ=1//1)
  m  = length(dW)

  WikJ = vec(dW) .* vec(dW)'
  WikA = zero(WikJ)

  p==nothing && (p = Int(floor(C*dt^(1//1-2//1*γ)) + 1))

  # Eq.(20)
  ρp = (π^2)/6
  for k in 1:p
    ρp -= (1/k^2)
  end
  ρp = ρp/(2*π^2)

  if dZ==nothing
    μ = randn(eltype(dW),m)
  else
    μ = convert(eltype(dW),sqrt(1/dt))*dZ[1:m]
  end

  for k in 1:p
    if dZ==nothing
      ζ = randn(eltype(dW),m)
      η = randn(eltype(dW),m)
    else
      ζ = convert(eltype(dW),sqrt(1/dt))*(dZ[Int(m*(2*k-1)+1):Int(m*2*k)])
      η = convert(eltype(dW),sqrt(1/dt))*(dZ[Int(m*2*k+1):Int(m*(2*k+1))])
    end

    kronprod1 = vec(ζ) .* vec(η+convert(eltype(dW),sqrt(2/dt))*dW)'
    kronprod2 = vec(η+convert(eltype(dW),sqrt(2/dt))*dW) .* vec(ζ)'

    WikA += (kronprod1 - kronprod2)/k
  end

  WikA *=  dt/(2*π)

  kronprod1 = vec(μ) .* vec(dW)'
  kronprod2 = vec(dW).* vec(μ)'

  WikJ = WikJ/2 + WikA + sqrt(dt*ρp)*(kronprod1-kronprod2)  # Ito: (UniformScaling(dt))/2
  return WikJ
end

function get_iterated_I!(dt, dW, dZ, Wik::KPWJ_iip, p=nothing, C=1, γ=1//1)
  m  = length(dW)
  @unpack WikA,WikJ,μ,ζ,η,kronprod1,kronprod2 = Wik

  mul!(WikJ,vec(dW),vec(dW)')
  fill!(WikA, zero(eltype(WikJ)))

  # Below Eq.(26): truncation
  p==nothing && (p = Int(floor(C*dt^(1//1-2//1*γ)) + 1))
  #@show p, dZ
  # Eq.(20)
  ρp = (π^2)/6
  for k in 1:p
    ρp -= (1/k^2)
  end
  ρp = ρp/(2*π^2)

  if dZ==nothing
    randn!(μ)
  else
    μ .= convert(eltype(dW),sqrt(1/dt))*@view(dZ[1:m])
  end

  # Eq. (21)
  for k in 1:p
    if dZ==nothing
      randn!(ζ)
      randn!(η)
    else
      ζ .= convert(eltype(dW),sqrt(1/dt))*@view(dZ[Int(m*(2*k-1)+1):Int(m*2*k)])
      η .= convert(eltype(dW),sqrt(1/dt))*@view(dZ[Int(m*2*k+1):Int(m*(2*k+1))])
    end
    mul!(kronprod1,vec(ζ),vec(η+convert(eltype(dW),sqrt(2/dt))*dW)')
    mul!(kronprod2,vec(η+convert(eltype(dW),sqrt(2/dt))*dW),vec(ζ)')

    @.. WikA = WikA + (kronprod1 - kronprod2)/k
  end

  @.. WikA = WikA*dt/(2*π) # check if factor of 2 is correct
  mul!(kronprod1,vec(μ),vec(dW)')
  mul!(kronprod2,vec(dW),vec(μ)')
  #@show WikA, sqrt(dt*ρp)*(kronprod1-kronprod2) , WikA + sqrt(dt*ρp)*(kronprod1-kronprod2)
  #@show sqrt(dt*ρp)
  @.. WikJ = WikJ/2 + WikA + sqrt(dt*ρp)*(kronprod1-kronprod2) # Ito: -UniformScaling(dt))/2
  return nothing
end

# Wiktorssson Approximation

abstract type AbstractWikJ <: AbstractJ end
abstract type AbstractWikJDiagonal <: AbstractWikJ end
abstract type AbstractWikJCommute <: AbstractWikJ end
abstract type AbstractWikJGeneral <: AbstractWikJ end

struct WikJDiagonal_oop <: AbstractWikJDiagonal end

mutable struct WikJDiagonal_iip{WikJType} <: AbstractWikJDiagonal
  WikJ::WikJType
  WikJDiagonal_iip(ΔW) = new{typeof(ΔW)}(false .* ΔW .* ΔW)
end

struct WikJCommute_oop <: AbstractWikJCommute end

mutable struct WikJCommute_iip{WikJType} <: AbstractWikJCommute
  WikJ::WikJType
  function WikJCommute_iip(ΔW)
    WikJ = false .* vec(ΔW) .* vec(ΔW)'
    new{typeof(WikJ)}(WikJ)
  end
end

struct WikJGeneral_oop <: AbstractWikJGeneral
  m_seq::Matrix{Int}
  function WikJGeneral_oop(ΔW)
    m = length(ΔW)
    M = div(m*(m-1),2)
    m_seq = Matrix{Int}(undef, M, 2)
    k = 1
    for i in 1:length(ΔW)
      for j in i+1:length(ΔW)
        m_seq[k,1] = i
        m_seq[k,2] = j
        k += 1
      end
    end
    new(m_seq)
  end
end

mutable struct WikJGeneral_iip{rateNoiseElTypeNoUnits, WikJType} <: AbstractWikJGeneral
  WikJ::WikJType
  WikJ2::WikJType
  WikJ3::WikJType
  m_seq::Matrix{Int}
  vec_ζ::Vector{rateNoiseElTypeNoUnits}
  vec_η::Vector{rateNoiseElTypeNoUnits}
  Gp₁::Vector{rateNoiseElTypeNoUnits}
  Gp₂::Vector{rateNoiseElTypeNoUnits}
  Aᵢ::Vector{rateNoiseElTypeNoUnits}
end

function WikJGeneral_iip(ΔW)
  WikJ = false .* vec(ΔW) .* vec(ΔW)'
  WikJ2 = false .* vec(ΔW) .* vec(ΔW)'
  WikJ3 = false .* vec(ΔW) .* vec(ΔW)'
  m = length(ΔW)
  M = div(m*(m-1),2)
  m_seq = Array{Int}(undef, M, 2)
  k = 1
  for i in 1:length(ΔW)
    for j in i+1:length(ΔW)
      m_seq[k,1] = i
      m_seq[k,2] = j
      k += 1
    end
  end
  vec_ζ = zero(ΔW)
  vec_η = zero(ΔW)
  Gp₁ = false .* Array{eltype(ΔW)}(undef, M)
  Gp₂ = false .* Array{eltype(ΔW)}(undef, M)
  Aᵢ = zero(ΔW)
  WikJGeneral_iip{eltype(ΔW), typeof(WikJ)}(WikJ, WikJ2, WikJ3, m_seq, vec_ζ, vec_η, Gp₁, Gp₂, Aᵢ)
end

function get_WikJ(ΔW,prob,alg)
  if isinplace(prob)
    if typeof(ΔW) <: Number || is_diagonal_noise(prob)
      return WikJDiagonal_iip(ΔW)
    elseif alg.ii_approx isa IICommutative
      return WikJCommute_iip(ΔW)
    else
      return KPWJ_iip(ΔW) # WikJGeneral_iip(ΔW)
    end
  else
    if typeof(ΔW) <: Number || is_diagonal_noise(prob)
      return WikJDiagonal_oop()
    elseif alg.ii_approx isa IICommutative
      return WikJCommute_oop()
    else
      return KPWJ_oop() #WikJGeneral_oop(ΔW)
    end
  end
end

function get_WikJ(ΔW,prob,alg::RKMilCommute)
  if isinplace(prob)
    if typeof(ΔW) <: Number || is_diagonal_noise(prob)
      return WikJDiagonal_iip(ΔW)
    else
      return WikJCommute_iip(ΔW)
    end
  else
    if typeof(ΔW) <: Number || is_diagonal_noise(prob)
      return WikJDiagonal_oop()
    else
      return WikJCommute_oop()
    end
  end
end

function get_iterated_I(dt, dW, dZ, Wik::WikJDiagonal_oop, p=nothing, C=1, γ=1//1)
  WikJ = 1//2 .* dW .* dW
  WikJ
end

function get_iterated_I!(dt, dW, dZ, Wik::WikJDiagonal_iip, p=nothing, C=1, γ=1//1)
  @unpack WikJ = Wik
  if typeof(dW) <: Number
    Wik.WikJ = 1//2 .* dW .^ 2
  else
    @.. WikJ = 1//2*dW^2
  end
  return nothing
end

function get_iterated_I(dt, dW, dZ, Wik::WikJCommute_oop, p=nothing, C=1, γ=1//1)
  WikJ = 1//2 .* vec(dW) .* vec(dW)'
  WikJ
end

function get_iterated_I!(dt, dW, dZ, Wik::WikJCommute_iip, p=nothing, C=1, γ=1//1)
  @unpack WikJ = Wik
  mul!(WikJ,vec(dW),vec(dW)')
  @.. WikJ *= 1//2
  return nothing
end

"""

    get_iterated_I(dt, dW, Wik::WikJGeneral_oop, C=1)

This function calculates WikJ, a mxm Array for a m dimensional general noise problem, which is a approximation
to the second order iterated integrals.

WikJDiagonal and WikJCommute use the properties of respective noises to simplify the calculations.
While the calculation for General Noise case is taken from section 4 of [SDELab: A Package for solving stochastic differential
equations in MATLAB](https://doi.org/10.1016/j.cam.2006.05.037) and SDELAB2(https://github.com/tonyshardlow/SDELAB2)
which is the Implementation of SDELab in Julia.
```math
    𝒜ᵖ = (Iₘ² - Pₘ)Kₘᵀ Δt/(2π) √(𝑎ₚ) √(Σ∞) Gp₁
```

```math
    √(Σ∞) = (Σ∞ + 2αIₘ)/(√2 * (1 + α))
```

let the combined operators be,
```math
    F = Kₘ(Iₘ² - Pₘ)(Iₘ ⨂ W(Δt)W(Δt)ᵀ)(Iₘ² - Pₘ)Kₘᵀ
```

```math
    Σ∞ = 2Iₘ + (2/Δt)F
```

See the paper for further details of specific operators.
Here we've only shown in which order these are implemented in this code.

From above we can see:

```math
    Δt/(2π) √(𝑎ₚ) √(Σ∞) Gp₁ = Δt/π √(𝑎ₚ) (√(Σ∞)/2 Gp₁)
```

let,
```math
    Oper2(Gp₁) = (√(Σ∞)/2 Gp₁) = (Iₘ/√2 + F/(√2 * (1 + α) * Δt))(Gp₁)
    Oper2(Gp₁) = Iₘ*Gp₁/√2 + F(Gp₁/(√2*(1+α)*Δt))
    Oper2(Gp₁) = Gp₁/√2 + F(Gp₁/(√2*(1+α)*Δt))
```

we have,

    Gp₂ = Gp₁/(sqrt(2)*(1+α)*dt)

```math
    𝒜ᵖ = (Iₘ² - Pₘ)Kₘᵀ Δt/π √(𝑎ₚ) Oper2(Gp₁/√2 + F(Gp₂))
    𝒜ᵖ = √(𝑎ₚ)*Δt/π * (Iₘ² - Pₘ)Kₘᵀ Oper2(Gp₁/√2 + F(Gp₂))
```
In the code we have

```math
    WikJ2 = (Iₘ² - Pₘ)Kₘᵀ(Oper2(Gp₁))
```

"""
function get_iterated_I(dt, dW, dZ, Wik::WikJGeneral_oop, p=nothing, C=1, γ=1//1)
  @unpack m_seq = Wik
  m = length(dW)
  M = div(m*(m-1),2)
  sum_dW² = dot(dW,dW)

  WikJ = dW*dW'
  eltWikJ = eltype(WikJ)
  WikJ2 =zeros(eltWikJ, size(WikJ))
  Gp₁ = randn(eltWikJ,M)
  α = convert(eltWikJ,sqrt(1 + sum_dW²/dt))
  Gp₂ = Gp₁/convert(eltWikJ, sqrt(2)*(1+α)*dt)
  #operator (Iₘ² - Pₘ)Kₘᵀ
  for i in 1:M
    WikJ2[m_seq[i,1], m_seq[i,2]] = Gp₂[i]
    WikJ2[m_seq[i,2], m_seq[i,1]] = -Gp₂[i]
  end

  #operator (Iₘ X W*Wᵀ)
  WikJ2 = WikJ*WikJ2

  #operator Kₘ(Iₘ² - Pₘ)
  WikJ2 = WikJ2 - WikJ2'
  for i in 1:M
    Gp₂[i] = WikJ2[m_seq[i,1], m_seq[i,2]]
  end
  Gp₁ = Gp₁/convert(eltWikJ,sqrt(2))+Gp₂

  #operator (Iₘ² - Pₘ)Kₘᵀ
  for i in 1:M
    WikJ2[m_seq[i,1], m_seq[i,2]] = Gp₁[i]
    WikJ2[m_seq[i,2], m_seq[i,1]] = -Gp₁[i]
  end

  WikJ *= 1//2
  𝑎ₚ = (π^2)/6
  p==nothing && (p = Int(floor((1/(π))*sqrt(M/(12*dt*C))*sqrt(m + 4*sum_dW²/dt) + 1)))
  Aᵢ = false .* vec(dW)   # Aᵢ is vector of aᵢ₀
  for r in 1:p
    𝑎ₚ -= (1/r^2)
    var = convert(eltWikJ, sqrt(dt/(2*π*r)))
    vec_ζ = randn(eltWikJ,m)*var
    vec_η = randn(eltWikJ,m)*var
    WikJ += (vec_ζ*vec_η' - vec_η*vec_ζ')
    Aᵢ +=  convert(eltWikJ, 2/(sqrt(π*r)))*vec_ζ
  end

  WikJ -= 1//2*(dW*Aᵢ' - Aᵢ*dW')
  WikJ -= convert(eltWikJ, sqrt(𝑎ₚ)*dt/π)*WikJ2 #-= to agree with operator version from Wiktorssson
  WikJ
end

"""

    get_iterated_I!(dt, dW, Wik::WikJGeneral_iip, C=1)

This function calculates WikJ, a mxm Array for a m dimensional general noise problem, which is a approximation
to the second order iterated integrals.

WikJDiagonal and WikJCommute use the properties of respective noises to simplify the calculations.
While the calculation for General Noise case is taken from section 4 of [SDELab: A Package for solving stochastic differential
equations in MATLAB](https://doi.org/10.1016/j.cam.2006.05.037) and SDELAB2(https://github.com/tonyshardlow/SDELAB2)
which is the Implementation of SDELab in Julia.
```math
    𝒜ᵖ = (Iₘ² - Pₘ)Kₘᵀ Δt/(2π) √(𝑎ₚ) √(Σ∞) Gp₁
```

```math
    √(Σ∞) = (Σ∞ + 2αIₘ)/(√2 * (1 + α))
```

let the combined operators be,
```math
    F = Kₘ(Iₘ² - Pₘ)(Iₘ ⨂ W(Δt)W(Δt)ᵀ)(Iₘ² - Pₘ)Kₘᵀ
```

```math
    Σ∞ = 2Iₘ + (2/Δt)F
```

See the paper for further details of specific operators.
Here we've only shown in which order these are implemented in this code.

From above we can see:

```math
    Δt/(2π) √(𝑎ₚ) √(Σ∞) Gp₁ = Δt/π √(𝑎ₚ) (√(Σ∞)/2 Gp₁)
```

let,
```math
    Oper2(Gp₁) = (√(Σ∞)/2 Gp₁) = (Iₘ/√2 + F/(√2 * (1 + α) * Δt))(Gp₁)
    Oper2(Gp₁) = Iₘ*Gp₁/√2 + F(Gp₁/(√2*(1+α)*Δt))
    Oper2(Gp₁) = Gp₁/√2 + F(Gp₁/(√2*(1+α)*Δt))
```

we have,

    Gp₂ = Gp₁/(sqrt(2)*(1+α)*dt)

```math
    𝒜ᵖ = (Iₘ² - Pₘ)Kₘᵀ Δt/π √(𝑎ₚ) Oper2(Gp₁/√2 + F(Gp₂))
    𝒜ᵖ = √(𝑎ₚ)*Δt/π * (Iₘ² - Pₘ)Kₘᵀ Oper2(Gp₁/√2 + F(Gp₂))
```
In the code we have

```math
    WikJ2 = (Iₘ² - Pₘ)Kₘᵀ(Oper2(Gp₁))
```

"""
function get_iterated_I!(dt, dW, dZ, Wik::WikJGeneral_iip, p=nothing, C=1, γ=1//1)
  @unpack WikJ, WikJ2, WikJ3, m_seq, vec_ζ, vec_η, Gp₁, Gp₂, Aᵢ = Wik

  m = length(dW)
  M = div(m*(m-1),2)

  sum_dW² = dot(dW,dW) #zero(eltype(dW))
  # mul!(sum_dW²,dW', dW)

  Gp₁ .= randn(M)
  α = sqrt(1 + sum_dW²/dt)
  @.. Gp₂ = Gp₁/(sqrt(2)*(1+α)*dt)

  #operator (Iₘ² - Pₘ)Kₘᵀ
  for i in 1:M
    WikJ2[m_seq[i,1], m_seq[i,2]] = Gp₂[i]
    WikJ2[m_seq[i,2], m_seq[i,1]] = -Gp₂[i]
  end

  #operator (Iₘ X W*Wᵀ)
  mul!(WikJ,dW,dW')
  mul!(WikJ3,WikJ,WikJ2)

  #operator Kₘ(Iₘ² - Pₘ)
  @.. WikJ2 = WikJ3 - WikJ3'
  for i in 1:M
    Gp₂[i] = WikJ2[m_seq[i,1], m_seq[i,2]]
  end
  @.. Gp₁ = Gp₁/sqrt(2) + Gp₂

  #operator (Iₘ² - Pₘ)Kₘᵀ
  for i in 1:M
    WikJ2[m_seq[i,1], m_seq[i,2]] = Gp₁[i]
    WikJ2[m_seq[i,2], m_seq[i,1]] = -Gp₁[i]
  end

  @.. WikJ *= 1//2
  𝑎ₚ = (π^2)/6
  p==nothing && (p = Int(floor((1/(π))*sqrt(M/(12*dt*C))*sqrt(m + 4*sum_dW²/dt) + 1)))
  Aᵢ .= false .* vec(dW)    # Aᵢ is vector of aᵢ₀
  for r in 1:p
    𝑎ₚ -= (1/r^2)
    var = sqrt(dt/(2*π*r))
    vec_ζ .= randn(m) .* var
    vec_η .= randn(m) .* var
    mul!(WikJ3, vec_ζ, vec_η')
    @.. WikJ += WikJ3 - WikJ3'
    @.. Aᵢ += (2/sqrt(π*r))*vec_ζ
  end
  mul!(WikJ3, dW, Aᵢ')
  @.. WikJ = WikJ - 1//2*(WikJ3 - WikJ3')
  @.. WikJ2 *= -(sqrt(𝑎ₚ)*dt/π) # minus sign needed to agree with operator version
  #@show WikJ, (sqrt(𝑎ₚ)*dt/π)*WikJ2
  @.. WikJ = WikJ + WikJ2
  return nothing
end


text:
module TestResampling

#  using Revise
using Test
using MLJ
using MLJBase
import MLJModels
import Random.seed!
seed!(1234)

include("foobarmodel.jl")

@test CV(nfolds=6) == CV(nfolds=6)
@test CV(nfolds=5) != CV(nfolds=6)
@test MLJ.train_test_pairs(CV(), 1:10) !=
     MLJ.train_test_pairs(CV(shuffle=true), 1:10)
@test MLJ.train_test_pairs(Holdout(), 1:10) !=
     MLJ.train_test_pairs(Holdout(shuffle=true), 1:10)

@testset "checking measure/model compatibility" begin
    model = ConstantRegressor()
    y = rand(4)
    override=false
    @test MLJ._check_measure(:junk, :junk, :junk, :junk, true) == nothing
    @test_throws ArgumentError MLJ._check_measure(model, rms, y, predict, override)
    @test MLJ._check_measure(model, rms, y, predict_mean, override) == nothing
    @test MLJ._check_measure(model, rms, y, predict_median, override) == nothing
    y=categorical(collect("abc"))
    @test_throws(ArgumentError,
                 MLJ._check_measure(model, rms, y, predict_median, override))
    model = ConstantClassifier()
    @test_throws(ArgumentError,
                 MLJ._check_measure(model, misclassification_rate, y, predict, override))
    @test MLJ._check_measure(model, misclassification_rate, y,
                            predict_mode, override) == nothing
    model = MLJModels.DeterministicConstantClassifier()
    @test_throws ArgumentError MLJ._check_measure(model, cross_entropy, y,
                            predict, override)
end

@testset "folds specified" begin
    x1 = ones(10)
    x2 = ones(10)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    my_rms(yhat, y) = sqrt(mean((yhat -y).^2))
    my_mav(yhat, y) = abs.(yhat - y)
    MLJ.reports_each_observation(::typeof(my_mav)) = true

    resampling = [(3:10, 1:2),
                  ([1, 2, 5, 6, 7, 8, 9, 10], 3:4),
                  ([1, 2, 3, 4, 7, 8, 9, 10], 5:6),
                  ([1, 2, 3, 4, 5, 6, 9, 10], 7:8),
                  (1:8, 9:10)]

    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)

    # check detection of incompatible measure (cross_entropy):
    @test_throws ArgumentError evaluate!(mach, resampling=resampling,
                                         measure=[cross_entropy, rmslp1])
    result = evaluate!(mach, resampling=resampling,
                       measure=[my_rms, my_mav, rmslp1])
    v = [1/2, 3/4, 1/2, 3/4, 1/2]
    @test result.per_fold[1] ≈ v
    @test result.per_fold[2] ≈ v
    @test result.per_fold[3][1] ≈ abs(log(2) - log(2.5))
    @test ismissing(result.per_observation[1])
    @test result.per_observation[2][1] ≈ [1/2, 1/2]
    @test result.per_observation[2][2] ≈ [3/4, 3/4]
    @test result.measurement[1] ≈ mean(v)
    @test result.measurement[2] ≈ mean(v)
end

@testset "holdout" begin
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0]

    @test MLJBase.show_as_constructed(Holdout)
    holdout = Holdout(fraction_train=0.75)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    result = evaluate!(mach, resampling=holdout,
                       measure=[rms, rmslp1])
    result = evaluate!(mach, verbosity=0, resampling=holdout)
    result.measurement[1] ≈ 2/3

    # test direct evaluation of a model + data:
    result = evaluate(model, X, y, verbosity=0, resampling=holdout, measure=rms)
    @test result.measurement[1] ≈ 2/3

    X = (x=rand(100),)
    y = rand(100)
    mach = machine(model, X, y)
    evaluate!(mach, verbosity=0, Holdout(shuffle=true, rng=123))
    e1 = evaluate!(mach, verbosity=0, Holdout(shuffle=true)).measurement[1]
    @test e1 != evaluate!(mach, verbosity=0, Holdout()).measurement[1]

end

@testset "cv" begin
    x1 = ones(10)
    x2 = ones(10)
    X = (x1=x1, x2=x2)
    y = [1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0]

    @test MLJBase.show_as_constructed(CV)
    cv=CV(nfolds=5)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    result = evaluate!(mach, resampling=cv, measure=[rms, rmslp1])
    @test result.per_fold[1] ≈ [1/2, 3/4, 1/2, 3/4, 1/2]

    shuffled =  evaluate!(mach, resampling=CV(shuffle=true)) # using rms default
    @test shuffled.measurement[1] != result.measurement[1]
end

@testset "weights" begin

    # cv:
    x1 = ones(4)
    x2 = ones(4)
    X = (x1=x1, x2=x2)
    y = [1.0, 2.0, 3.0, 1.0]
    w = 1:4
    cv=CV(nfolds=2)
    model = MLJModels.DeterministicConstantRegressor()
    mach = machine(model, X, y)
    e = evaluate!(mach, resampling=cv, measure=l1,
                  weights=w, verbosity=0).measurement[1]
    @test e ≈ (1/3 + 13/14)/2

end

@testset "resampler as machine" begin
    N = 50
    X = (x1=rand(N), x2=rand(N), x3=rand(N))
    y = X.x1 -2X.x2 + 0.05*rand(N)

    ridge_model = FooBarRegressor(lambda=20.0)
    holdout = Holdout(fraction_train=0.75)
    resampler = Resampler(resampling=holdout, model=ridge_model, measure=mav)
    resampling_machine = machine(resampler, X, y)
    @test_logs((:info, r"^Training"), fit!(resampling_machine))
    e1=evaluate(resampling_machine).measurement[1]
    mach = machine(ridge_model, X, y)
    @test e1 ≈  evaluate!(mach, resampling=holdout,
                          measure=mav, verbosity=0).measurement[1]
    ridge_model.lambda=1.0
    fit!(resampling_machine, verbosity=2)
    e2=evaluate(resampling_machine).measurement[1]
    @test e1 != e2
    resampler.weights = rand(N)
    fit!(resampling_machine, verbosity=0)
    e3=evaluate(resampling_machine).measurement[1]
    @test e3 != e2

    @test MLJBase.package_name(Resampler) == "MLJ"
    @test MLJBase.is_wrapper(Resampler)
    rnd = randn(5)
    @test evaluate(resampler, rnd) === rnd
end

struct DummyResamplingStrategy <: MLJ.ResamplingStrategy end

@testset "custom strategy using resampling depending on X, y" begin
    function MLJ.train_test_pairs(resampling::DummyResamplingStrategy,
                              rows, X, y)
        train = filter(rows) do j
            y[j] == y[1]
        end
        test = setdiff(rows, train)
        return [(train, test),]
    end

    X = (x = rand(8), )
    y = categorical([:x, :y, :x, :x, :y, :x, :x, :y])
    @test MLJ.train_test_pairs(DummyResamplingStrategy(), 2:6, X, y) ==
        [([3, 4, 6], [2, 5]),]

    e = evaluate(ConstantClassifier(), X, y,
                 measure=misclassification_rate,
                 resampling=DummyResamplingStrategy(),
                 operation=predict_mode)
    @test e.measurement[1] ≈ 1.0
end

end
true