text:
module DataFrames

##############################################################################
##
## Dependencies
##
##############################################################################

using Base.Intrinsics
using Reexport
@reexport using StatsBase
@reexport using DataArrays
using GZip
using SortingAlgorithms

##############################################################################
##
## Extend methods in Base by default
##
##############################################################################

import Base: Sort, Order
import Base.Sort: sort, sort!, sortperm, Algorithm, defalg, issorted
import Base.Order: Ordering, By, Lt, Perm, Forward, lt, ord
import Base.AsyncStream

##############################################################################
##
## Exported methods and types
##
##############################################################################

export @~,
       AbstractDataFrame,
       array,
       by,
       cbind,
       coefnames,
       colwise,
       combine,
       complete_cases,
       complete_cases!,
       DataFrame,
       DataFrameRow,
       deleterows!,
       describe,
       eachcol,
       eachrow,
       eltypes,
       flipud!,
       flipud,
       Formula,
       gl,
       GroupApplied,
       groupby,
       GroupedDataFrame,
       interaction_design_matrix,
       loaddf,
       model_response,
       ModelFrame,
       ModelMatrix,
       names!,
       ncol,
       nonunique,
       nrow,
       nullable!,
       order,
       pool,
       pool!,
       printtable,
       rbind,
       readtable,
       rename!,
       rename,
       save,
       showcols,
       stack,
       SubDataFrame,
       subset,
       unique,
       unique!,
       unstack,
       writetable,
       xtab,
       xtabs,
       stackdf,
       melt,
       meltdf,
       pivottable,
       read_rda,
       vecbind,
       push!

##############################################################################
##
## Load files
##
##############################################################################

include(joinpath("other", "utils.jl"))
include(joinpath("other", "index.jl"))
include(joinpath("abstractdataframe", "abstractdataframe.jl"))
include(joinpath("dataframe", "dataframe.jl"))
include(joinpath("subdataframe", "subdataframe.jl"))
include(joinpath("dataframerow", "dataframerow.jl"))
include(joinpath("dataframe", "sort.jl"))
include(joinpath("dataframe", "iteration.jl"))
include(joinpath("dataframe", "show.jl"))
include(joinpath("dataframe", "join.jl"))
include(joinpath("groupeddataframe", "grouping.jl"))
include(joinpath("dataframe", "reshape.jl"))
include(joinpath("statsmodels", "formula.jl"))
include(joinpath("statsmodels", "statsmodel.jl"))
include(joinpath("dataframe", "io.jl"))
include("RDA.jl")
include("deprecated.jl")

end # module DataFrames


text:
using DemultiplexFASTQ
processargs()

ext: jl


text:
using IntensityScans
using IntensityScans.ThorlabsLTStage
using IntensityScans.TcpInstruments
using IntensityScans.Unitful

using Test

include("emulate_scan.jl")

@testset "IntensityScans.jl" begin
    include("test_scans.jl")
    include("test_metrics.jl")
    include("test_plotting.jl")
    include("test_loading_scans.jl")
end


text:
export tesselation

"""
`expander(P::HPolygon)` returns a list of polygons formed by reflecting `P`
across each of its sides.
"""
function expander(P::HPolygon)
    slist = sides(P)
    return [ reflect_across(P,S) for S in slist ]
end

"""
`tesselation(n,k,deep)`: Tesselate the hyperbolic plane by regular `n`-gons in which
each vertex is a corner of `k` polygons. `deep` controls how many layers.
The center of the first `k`-gon is placed at the origin.

May also be called `tesselation(n,k,deep,true)` in which case
a vertex of the first `k`-gon is placed at the origin and the tesselation
is seeded by copies of this first polygon around the origin.
"""
function tesselation(n::Int, k::Int, deep::Int=2,  vertex_centered::Bool=false)
    theta = 2pi/k
    P = equiangular(n,theta)
    outlist = HContainer()
    todo = HContainer()

    if vertex_centered
        v = P.plist[1]
        f = move2zero(v)
        P = f(P)
        r = rotation(theta)
        for j=1:k
            add_object!(todo,P)
            add_object!(outlist,P)
            P = r(P)
        end
    else
        add_object!(todo, P)
        add_object!(outlist,P)
    end

    for j=1:deep
        new_todo = HContainer()
        for X in todo.objs
            add_object!(outlist,X)
            for Y in expander(X)
                if !in(Y,outlist)
                    add_object!(new_todo,Y)
                end
            end
        end
        todo = new_todo
    end

    return outlist
end


text:
using EzXML

function getxml(; from="", until="")
    baseuri = "http://export.arxiv.org/oai2?verb=ListRecords"
    uri = "$baseuri&metadataPrefix=arXiv"
    if !isempty(from)
        @assert !isempty(until)
        uri = "$uri&from=$from&until=$until"
    end

    while true
        xml = readxml(download(uri))
        nodes = find(xml, "//resumptionToken")
        isempty(nodes) && break
        token = nodecontent(nodes[1])
        uri = "$baseuri&resumptionToken=$token"
    end
end

function readxml(path::String)
    lines = open(readlines, path)
    @assert startswith(lines[2],"<OAI-PMH")
    lines[2] = "<OAI-PMH>"
    for i = 1:length(lines)
        startswith(lines[i]," <arXiv") && (lines[i] = " <arXiv>")
    end
    parsexml(join(lines,"\n"))
end


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# test/streaminfo.jl: test stream information functions

@testset "StreamInfo" begin

info = StreamInfo(name="testname",
                  type="testtype",
                  channel_count=2,
                  nominal_srate=10,
                  channel_format=Float64,
                  source_id="testid")

@test name(info) == "testname"
@test type(info) == "testtype"
@test channel_count(info) == 2
@test nominal_srate(info) == 10
@test channel_format(info) == Float64
@test source_id(info) == "testid"

@test version(info) == protocol_version()

channels = append_child(desc(info), "channels")
for label in ["C3", "C4", "Cz", "FPz", "POz", "CPz", "O1", "O2"]
  ch = append_child(channels, "channel")
  append_child_value(ch, "label", label)
  append_child_value(ch, "unit", "microvolts")
  append_child_value(ch, "type", "EEG")
end
append_child_value(desc(info), "manufacturer", "SCCN")
cap = append_child(desc(info), "cap")
append_child_value(cap, "name", "EasyCap")
append_child_value(cap, "size", "54")
append_child_value(cap, "labelscheme", "10-20")

end

# Force cleanup
GC.gc()


text:
# functional definitions

import PDESolver: createFunctional, _evalFunctional, _evalFunctionalDeriv_m,
                 _evalFunctionalDeriv_q
import ODLCommonTools: getParallelData, setupFunctional

#------------------------------------------------------------------------------
# Lift and drag

@doc """
###EulerEquationMod.BoundaryForceData

Composite data type for storing data pertaining to the boundaryForce. It holds
lift and drag values

"""->

mutable struct BoundaryForceData{Topt, fname} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int,1}

  # factors to multiply x, y, z momenta by, determines if lift or drag is
  # calculated
  facx::Topt
  facy::Topt
  facz::Topt

  # things needed for the calculation
  qg::Vector{Topt}
  euler_flux::Vector{Topt}
end

"""
  Constructor for BoundaryForceData{Topt, :lift}
"""
function LiftForceDataConstructor(::Type{Topt}, mesh, sbp, eqn, opts, bcnums) where Topt

  facx = 0
  facy = 0
  facz = 0

  qg = zeros(Topt, mesh.numDofPerNode)
  euler_flux = zeros(Topt, mesh.numDofPerNode)

  return BoundaryForceData{Topt, :lift}(bcnums, facx, facy, facz, qg,
                                        euler_flux)
end

"""
  Constructor for BoundaryForceData{Topt, :drag}
"""
function DragForceDataConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                             bcnums) where Topt

  facx = 0
  facy = 0
  facz = 0

  qg = zeros(Topt, mesh.numDofPerNode)
  euler_flux = zeros(Topt, mesh.numDofPerNode)

  return BoundaryForceData{Topt, :drag}(bcnums, facx, facy, facz, qg,
                                        euler_flux)
end

function setupFunctional(mesh::AbstractMesh, sbp, eqn::AbstractSolutionData,
                         opts::Dict, func::BoundaryForceData{Topt, :lift}) where {Topt}

  if mesh.dim == 2
    func.facx = -sin(eqn.params.aoa)
    func.facy =  cos(eqn.params.aoa)
  else
    func.facx = -sin(eqn.params.aoa)
    func.facy =  0
    func.facz =  cos(eqn.params.aoa)
  end

  return nothing
end

function setupFunctional(mesh::AbstractMesh, sbp, eqn::AbstractSolutionData,
                         opts::Dict, func::BoundaryForceData{Topt, :drag}) where {Topt}

  if mesh.dim == 2
    func.facx = cos(eqn.params.aoa)
    func.facy = sin(eqn.params.aoa)
  else
    func.facx = cos(eqn.params.aoa)
    func.facy = 0
    func.facz = sin(eqn.params.aoa)
  end

  return nothing
end

"""
  Functional for computing lift coefficient.  Uses the lift functional to
  compute the force and then divides by the (non-dimensional) dynamic pressure
  0.5*rho_free*Ma^2.  Note that this assumes the chord length (in 2d) is 1
"""
mutable struct LiftCoefficient{Topt} <: AbstractBoundaryFunctional{Topt}
  func::BoundaryForceData{Topt, :lift}
  bcnums::Array{Int, 1}
end

"""
  Constructor for LiftCoefficient functional
"""
function LiftCoefficientConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                    bcnums) where Topt

  lift = LiftForceDataConstructor(Topt, mesh, sbp, eqn, opts, bcnums)

  return LiftCoefficient{Topt}(lift, bcnums)
end

"""
  Functional for computing drag coefficient.  See [`LiftCoefficient`](@ref)
  for the non-dimensionalization.
"""
mutable struct DragCoefficient{Topt} <: AbstractBoundaryFunctional{Topt}
  func::BoundaryForceData{Topt, :drag}
  bcnums::Array{Int, 1}
end

"""
  Constructor for DragCoefficient functional
"""
function DragCoefficientConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                    bcnums) where Topt

  drag = DragForceDataConstructor(Topt, mesh, sbp, eqn, opts, bcnums)

  return DragCoefficient{Topt}(drag, bcnums)
end

"""
  Typedef for functionals that compute aerodynamic coefficients (lift, drag)
"""
const AeroCoefficients{Topt} = Union{LiftCoefficient{Topt}, DragCoefficient{Topt}}

function setupFunctional(mesh::AbstractMesh, sbp, eqn::AbstractSolutionData,
                         opts::Dict, func::AeroCoefficients)

  setupFunctional(mesh, sbp, eqn, opts, func.func)

end

#------------------------------------------------------------------------------
# other functionals

"""
  Type for computing the mass flow rate over a boundary (integral rho*u dot n
  dGamma)
"""
mutable struct MassFlowData{Topt} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int, 1}
end

"""
  Constructor for MassFlowData.  This needs to have a different name from
  the type so it can be put in a dictionary
"""
function MassFlowDataConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                            bcnums) where Topt
  return MassFlowData{Topt}(bcnums)
end

"""
  Type for computing the entropy flux rate over a boundary (integral S * u_i dot n_i
  dGamma), where S is the entropy function (from the IR entropy variables).
"""
mutable struct EntropyFluxData{Topt} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int, 1}
end

"""
  Constructor for EntropyFlux.  This needs to have a different name from
  the type so it can be put in a dictionary
"""
function EntropyFluxConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                            bcnums) where Topt
  return EntropyFluxData{Topt}(bcnums)
end

"""
  Functional that computes function `w^T d(u)`, where `w` are the entropy
  variables and `d(u)` is the entropy stable dissipation computed by
  [`ELFPenaltyFaceIntegral`](@ref).  Note that this is an integral
  over interior faces and not boundary faces
"""
mutable struct EntropyDissipationData{Topt} <: EntropyPenaltyFunctional{Topt}
  func::ELFPenaltyFaceIntegral
  func_sparseface::FluxType
  func_sparseface_revq::FluxType_revq
  func_sparseface_revm::FluxType_revm
  do_face_term::Bool
  do_lps::Bool
  do_sc::Bool
end

"""
  Constructor for [`EntropyDissipationData`](@ref)

  This function takes `bcnums` as an argument for consistency with
  the boundary functional constructors, but doesn't use it.
"""
function EntropyDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  func = ELFPenaltyFaceIntegral(mesh, eqn)
  func_sparseface = LFPenalty()
  func_sparseface_revq = LFPenalty_revq()
  func_sparseface_revm = LFPenalty_revm()

  do_face_term = true
  do_lps = false
  do_sc = false

  return EntropyDissipationData{Topt}(func, func_sparseface,
          func_sparseface_revq, func_sparseface_revm,
          do_face_term, do_lps, do_sc)
end

function LPSDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                   bcnums) where Topt

  func = ELFPenaltyFaceIntegral(mesh, eqn)
  func_sparseface = ZeroFlux()
  func_sparseface_revq = ZeroFlux_revq()
  func_sparseface_revm = ZeroFlux_revm()

  do_face_term = false
  do_lps = true
  do_sc = false

  return EntropyDissipationData{Topt}(func, func_sparseface,
          func_sparseface_revq, func_sparseface_revm, do_face_term,
          do_lps, do_sc)
end

function SCDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                   bcnums) where Topt

  func = ELFPenaltyFaceIntegral(mesh, eqn)
  func_sparseface = ZeroFlux()
  func_sparseface_revq = ZeroFlux_revq()
  func_sparseface_revm = ZeroFlux_revm()

  do_face_term = false
  do_lps = false
  do_sc = true

  return EntropyDissipationData{Topt}(func, func_sparseface,
          func_sparseface_revq, func_sparseface_revm, do_face_term,
          do_lps, do_sc)
end

"""
  Returns the negative of [`EntropyDissipationData`](@ref).  That functional
  is always negative, so this one is always positive.
"""
mutable struct NegEntropyDissipationData{Topt} <: EntropyPenaltyFunctional{Topt}
  func::EntropyDissipationData{Topt}
end

"""
  Constructor for [`NegEntropyDissipationData`](@ref).  `bcnums` argument is
  unused.
"""
function NegEntropyDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  func = EntropyDissipationConstructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegEntropyDissipationData{Topt}(func)
end

function NegLPSDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  func = LPSDissipationConstructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegEntropyDissipationData{Topt}(func)
end

function NegSCDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  func = SCDissipationConstructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegEntropyDissipationData{Topt}(func)
end

"""
  Functional that computes function `w^T d(u)`, where `w` are the entropy
  variables and `d(u)` is the entropy stable dissipation computed by
  [`ELFPenaltyFaceIntegral`](@ref).  Note that this is an integral
  over interior faces and not boundary faces
"""
mutable struct EntropyJumpData{Topt} <: EntropyPenaltyFunctional{Topt}
  func::EntropyJumpPenaltyFaceIntegral
  do_face_term::Bool
  do_lps::Bool
  do_sc::Bool
end

"""
  Constructor for [`EntropyDissipationData`](@ref)

  This function takes `bcnums` as an argument for consistency with
  the boundary functional constructors, but doesn't use it.
"""
function EntropyJumpConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  func = EntropyJumpPenaltyFaceIntegral(mesh, eqn)
  do_face_term = true
  do_lps = false
  do_sc = false

  return EntropyJumpData{Topt}(func, do_face_term, do_lps, do_sc)
end

"""
  Computes entropy dissipation by computing

  -W^T * R_face + \int_\Gamma_interior (psi_L - psi_R)

  This works for flux functions that do not have the EC + dissipation form,
  and should be < 0 for any entropy stable scheme.
"""
mutable struct EntropyDissipation2Data{Topt} <: EntropyPenaltyFunctional{Topt}
  flux_functional::EntropyDissipationData{Topt}
  potential_flux::PotentialFlux
  potential_flux_revq::PotentialFlux_revq
  potential_flux_revm::PotentialFlux_revm
end

function EntropyDissipation2Constructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  flux_functional = EntropyDissipationConstructor(Topt, mesh, sbp, eqn, opts,
                                                  bcnums)

  # configure the functional to compute (wL - wR)^T F, where F is the flux
  # currently used by the discretization
  flux_functional.func_sparseface = eqn.flux_func
  flux_functional.func_sparseface_revq = eqn.flux_func_revq
  flux_functional.func_sparseface_revm = eqn.flux_func_revm

  potential_flux = FluxDict["PotentialFlux"]
  potential_flux_revq = FluxDict_revq["PotentialFlux"]
  potential_flux_revm = FluxDict_revm["PotentialFlux"]

  return EntropyDissipation2Data{Topt}(flux_functional, potential_flux,
                                       potential_flux_revq, potential_flux_revm)
end

mutable struct TotalEntropyDissipationData{Topt} <: EntropyPenaltyFunctional{Topt}
  bcnums::Vector{Int}
end

"""
  Computes the entropy dissipation of the entire residual (not just face
  terms), except for the boundaries sepcified by bcnums:

  w^T R - other boundary terms

"""
function TotalEntropyDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                       bcnums) where Topt

  return TotalEntropyDissipationData{Topt}(copy(bcnums))
end

"""
  Computes negative of EntropyDissipation2
"""
mutable struct NegEntropyDissipation2Data{Topt} <: EntropyPenaltyFunctional{Topt}
  func::EntropyDissipation2Data{Topt}
end

function NegEntropyDissipation2Constructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                           bcnums) where Topt
  func = EntropyDissipation2Constructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegEntropyDissipation2Data{Topt}(func)
end

mutable struct NegTotalEntropyDissipationData{Topt} <: EntropyPenaltyFunctional{Topt}
  func::TotalEntropyDissipationData{Topt}
end

function NegTotalEntropyDissipationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                           bcnums) where Topt
  func = TotalEntropyDissipationConstructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegTotalEntropyDissipationData{Topt}(func)
end

function getParallelData(obj::EntropyPenaltyFunctional)
  return PARALLEL_DATA_ELEMENT
end

function getParallelData(obj::Union{EntropyDissipationData, EntropyJumpData})
  if obj.do_face_term
    return PARALLEL_DATA_ELEMENT
  else
    return PARALLEL_DATA_NONE
  end
end

const NegEntropyDissipations = Union{NegEntropyDissipationData,
                                     NegEntropyDissipation2Data,
                                     NegTotalEntropyDissipationData}

function getParallelData(obj::NegEntropyDissipations)
  return getParallelData(obj.func)
end

"""
  Functional that computes -delta_w^T lambda_max Y^T Y delta_w for
  boundary faces.  delta_w is the difference between the (entropy variable)
  numerical solution and the boundary state.  This only works for boundary
  conditions that have [`getDirichletState`](@ref) defined.
"""
mutable struct BoundaryEntropyDiss{Topt} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int, 1}
end

function BoundaryEntropyDissConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                            bcnums) where Topt
  return BoundaryEntropyDiss{Topt}(bcnums)
end

"""
  Negative of [`BoundaryEntropyDiss`](@ref).  That one is always negative
  so this one is always positive.
"""
mutable struct NegBoundaryEntropyDiss{Topt} <: AbstractBoundaryFunctional{Topt}
  func::BoundaryEntropyDiss{Topt}
end

function NegBoundaryEntropyDissConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                            bcnums) where Topt
  obj = BoundaryEntropyDissConstructor(Topt, mesh, sbp, eqn, opts, bcnums)
  return NegBoundaryEntropyDiss{Topt}(obj)
end

# this is actually a volume integral, not a boundary integral, but it
# doesn't make a difference
"""
  Computes the function from the Zahr and Perssons paper "An Optimization-based
  approach for high-order accurate discretization of conservation laws with
  discontinuous solutions", ie.

  \\sum_k \\int || u - h_bar ||^2 dOmega_k
"""
mutable struct SolutionDeviation{Topt} <: AbstractBoundaryFunctional{Topt}
end

function SolutionDeviationConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                            bcnums) where Topt
  return SolutionDeviation{Topt}()
end

"""
  Creates a functional object.

**Arguments**

 * `mesh` : Abstract PUMI mesh
 * `sbp`  : Summation-by-parts operator
 * `eqn`  : Euler equation object
 * `opts` : Options dictionary
 * `functional_name`: the name of the functional (in [`FunctionalDict`](@ref)
 * `functional_bcs`: the boundary condition numbers the functional is
                     computed on.
"""
function createFunctional(mesh::AbstractMesh, sbp::AbstractOperator,
                  eqn::EulerData{Tsol}, opts,
                  functional_name::AbstractString,
                  functional_bcs::Vector{I}) where {Tsol, I<:Integer}

  func_constructor = FunctionalDict[functional_name]
  objective = func_constructor(Tsol, mesh, sbp, eqn, opts, functional_bcs)

  return objective
end

"""
  Maps functional names to their outer constructors.

  All outer constructors must have the signature

  MyTypeNameConstructor{Topt}(::Type{Topt}, mesh, sbp, eqn, opts, bcnums)

  where MyTypeName is the name of the type, bcnums are the
  boundary conditions that the functional is
  defined on (Array{Int, 1}),

  Currently only boundary functionals are supported.

  For non-boundary functionals
  the region numbers associated with the faces would also be needed.
  Consider:

                | -> edge 1
     face 1     |         face2
                |
                |

  The mesh edges lying on geometric edge 1 have two possible parent elements,
  one of face 1 and one on face 2.  `geom_regions` picks between them.  This
  effectively determines the direction of the normal vector.

  Note that the outer constructor name is the type name with the suffix "Constructor"
"""
global const FunctionalDict = Dict{String, Function}(
"lift" => LiftForceDataConstructor,
"drag" => DragForceDataConstructor,
"liftCoefficient" => LiftCoefficientConstructor,
"dragCoefficient" => DragCoefficientConstructor,
"massflow" => MassFlowDataConstructor,
"entropyflux" => EntropyFluxConstructor,
"entropydissipation" => EntropyDissipationConstructor,
"lpsdissipation" => LPSDissipationConstructor,
"scdissipation" => SCDissipationConstructor,
"entropydissipation2" => EntropyDissipation2Constructor,
"totalentropydissipation" => TotalEntropyDissipationConstructor,
"negentropydissipation" => NegEntropyDissipationConstructor,
"neglpsdissipation" => NegLPSDissipationConstructor,
"negscdissipation" => NegSCDissipationConstructor,
"negtotalentropydissipation" => NegTotalEntropyDissipationConstructor,
"negentropydissipation2" => NegEntropyDissipation2Constructor,
"entropyjump" => EntropyJumpConstructor,
"boundaryentropydiss" => BoundaryEntropyDissConstructor,
"negboundaryentropydiss" => NegBoundaryEntropyDissConstructor,
"solutiondeviation" => SolutionDeviationConstructor,
)


text:
module Utils

using AsyncFinalizers
using Test

function outlined(f)
    @noinline wrapper() = f()
    wrapper()
end

function check_executor()
    @assert !istaskdone(AsyncFinalizers.Internal.EXECUTOR[])
end

return_nothing(_...) = nothing

@noinline function collect_garbage(;
    object = () -> Ref(0),
    destruct = return_nothing,
    async_finalizer = return_nothing,
    npolls = 1_000_000_000,
    check_executor::Bool = true,
)
    check_executor && (@__MODULE__).check_executor()
    phase1 = Threads.Atomic{Int}(0)
    phase2 = Threads.Atomic{Int}(0)
    outlined() do
        local obj = object()
        AsyncFinalizers.onfinalize(obj) do obj
            phase1[] = 1
            local y = destruct(obj)
            function ()
                phase2[] = 1
                async_finalizer(y)
            end
        end
        return
    end
    GC.gc()

    n = 0
    while phase1[] == 0
        GC.safepoint()
        yield()
        if (n += 1) > npolls
            error("too many pollings. finalizer not called?")
        end
    end
    @test phase1[] == 1

    n = 0
    while phase2[] == 0
        GC.safepoint()
        yield()
        if (n += 1) > npolls
            error("too many pollings. async finalizer not called?")
        end
    end
    @test phase2[] == 1
end

end  # module


text:
export LSForcingModel, NoLSForcing, HadGEMVertical

abstract type LSForcingModel end

vars_state(::LSForcingModel, ::AbstractStateType, FT) = @vars()

function compute_gradient_argument!(
    lsforcing::LSForcingModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

function compute_gradient_flux!(
    lsforcing::LSForcingModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

"""
    NoLSForcing <: LSForcingModel
No large-scale forcing
"""
struct NoLSForcing <: LSForcingModel end

"""
    Container for GCM variables from HadGEM2-A forcing,
used in the AMIP experiments
"""
struct HadGEMVertical <: LSForcingModel end

vars_state(m::HadGEMVertical, ::Auxiliary, FT) = @vars(
    ta::FT,
    hus::FT,
    ua::FT,
    va::FT,
    tntΣhava::FT,
    Σtemp_tendency::FT,
    Σqt_tendency::FT,
    w_s::FT,
)

vars_state(::HadGEMVertical, ::Gradient, FT) = @vars(ta::FT, hus::FT)
vars_state(::HadGEMVertical, ::GradientFlux, FT) = @vars(∇ᵥta::FT, ∇ᵥhus::FT)

function compute_gradient_argument!(
    lsforcing::HadGEMVertical,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.lsforcing.ta = aux.lsforcing.ta
    transform.lsforcing.hus = aux.lsforcing.hus
end

function compute_gradient_flux!(
    lsforcing::HadGEMVertical,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.lsforcing.∇ᵥta = ∇transform.lsforcing.ta[3]
    diffusive.lsforcing.∇ᵥhus = ∇transform.lsforcing.hus[3]
end


text:
boxcoordtype(::Type{T}) where T<:AbstractFloat = T
boxcoordtype(::Type{T}) where T<:Real = Float64

# The overall domain of the problem, plus initial evaluation locations
# and the function value at the initial point. Every box within the
# world will hold a link to this structure
struct World{T<:Real,P,M}
    lower::Vector{T}      # lower bounds on parameters
    upper::Vector{T}      # upper bounds on parameters
    position::Vector{P}   # initial position data along each dimension
    meta::M               # metadata at baseposition (see below)

    function World{T,P,M}(lower, upper, position, meta) where {T,P,M}
        # Validate the inputs
        N = length(lower)
        length(lower) == length(upper) == length(position) == N ||
            throw(DimensionMismatch("lower, upper, and position must have length $N (got $(length(lower)), $(length(upper)), and $(length(position))"))
        for i = 1:N
            world_validate(position[i], lower[i], upper[i])
        end
        return new{T,P,M}(lower, upper, position, meta)
    end
end

"""
    world = World(lower::AbstractVector, upper::AbstractVector, position::AbstractVector, meta)

Return a structure representing the rectangular domain and
initialization of a CoordinateSplittingPTree.  `lower` and `upper` are
the bounds of the coordinates, `position` is a vector that contains
position information---in the simplest case this may be just the
coordinates of an initial point. `meta` is the metadata at `position`.

In general, `position[i]` can contain whatever information you require
about dimension `i`. Packages may specialize the behavior of `World`
for different types of information---see the source for details.
"""
function World(lower::AbstractVector, upper::AbstractVector, position::AbstractVector, meta)
    T = promote_type(eltype(lower), eltype(upper))
    for x in position
        T = promote_type(T, world_eltype(x))
    end
    Tb = boxcoordtype(T)
    positionb = world_ofeltype.(Tb, position)
    metaval = metagen(meta, baseposition(positionb))
    return World{Tb,eltype(positionb),typeof(metaval)}(lower, upper, positionb, metaval)
end

function World(position::AbstractVector{T}, meta) where T<:Real
    Tb = boxcoordtype(T)
    n = length(position)
    return World(fill(Tb(-Inf), n), fill(Tb(Inf), n), position, meta)
end

metagen(meta, x) = meta
metagen(f::Function, x) = f(x)

## API for extending World:
#  - `baseposition([T], position)` must return an actual position,
#    i.e., a vector of real-valued numbers that lie between `lower`
#    and `upper`. If `T` is supplied it must be used as the eltype of
#    the returned vector.
#  - `baseposition(position[i])` should return the position in dimension `i`.
#  - `world_eltype(position[i])` must return a suitable element type
#    to be used by `baseposition`
#  - `world_ofeltype(T, position[i])` must convert to world_eltype `T`
#  - `world_validate(position[i], lower[i], upper[i])` should throw an
#    error if `position[i]` is invalid, and otherwise return `nothing`.
#  - `world_newposition(x, position[i], lower[i], upper[i])` should
#    return a valid location different from `x`.

baseposition(x::Real) = x
baseposition(position::AbstractVector{<:Real}) = position
baseposition(::Type{T}, position::AbstractVector{<:Real}) where T<:Real =
    T[T(x) for x in position]
world_eltype(x::T) where T<:Real = T
world_ofeltype(::Type{T}, x::Real) where T<:Real = T(x)
function world_validate(x::Real, l::Real, u::Real)
    l <= x <= u || throw(ArgumentError("position $x is not within bounds [$l, $u]"))
    return nothing
end
function world_newposition(x, s::Real, l::Real, u::Real)
    if x != s
        return s
    end
    return x + 1 < u ? x + 1 :
           x - 1 >= l ? x - 1 : error("no valid value found. Consider using `partition_interval` to find a new point.")
end

# We're also going to "reserve" the behavior of World when position[i]
# is a 2-tuple (as used by QuadSplit.jl). Since tuples are Base types,
# doing it in this package avoids risk of type piracy.

baseposition(s::Tuple{Real,Real}) = s[1]
baseposition(::Type{T}, splits::AbstractVector{TT}) where {T,TT<:Tuple{Real,Real}} =
    T[s[1] for s in splits]
baseposition(splits::AbstractVector{TT}) where TT<:Tuple{Real,Real} =
    [s[1] for s in splits]

world_eltype(v::Tuple{Real,Real}) = promote_type(typeof(v[1]), typeof(v[2]))
world_ofeltype(::Type{T}, v::Tuple{Real,Real}) where T = (T(v[1]), T(v[2]))

function world_validate(s::Tuple{Real,Real}, l::Real, u::Real)
    s[1] != s[2] || throw(ArgumentError("split points must be distinct, got $s"))
    l <= s[1] <= u || throw(ArgumentError("position $(s[1]) is not within bounds [$l, $u]"))
    l <= s[2] <= u || throw(ArgumentError("position $(s[2]) is not within bounds [$l, $u]"))
    return nothing
end

world_newposition(x, s::Tuple{Real,Real}, l::Real, u::Real) = x == s[1] ? s[2] : s[1]

function World(lower::AbstractVector{T}, upper::AbstractVector{T}, splits::AbstractVector{Tuple{T,T}}, meta) where T<:Real
    Tb = boxcoordtype(T)
    metaval = metagen(meta, baseposition(Tb, splits))
    return World{Tb, Tuple{Tb,Tb}, typeof(metaval)}(lower, upper, splits, metaval)
end

function World(lower::AbstractVector{Tl}, upper::AbstractVector{Tu}, splits::AbstractVector{Tuple{T,T}}, meta) where {Tl<:Real,Tu<:Real,T<:Real}
    Tb = boxcoordtype(promote_type(Tl,Tu,T))
    metaval = metagen(meta, baseposition(Tb, splits))
    return World{Tb, Tuple{Tb,Tb}, typeof(metaval)}(lower, upper, splits, metaval)
end

function World(splits::AbstractVector{Tuple{T,T}}, meta) where T<:Real
    Tb = boxcoordtype(T)
    metaval = metagen(meta, baseposition(Tb, splits))
    n = length(splits)
    return World{Tb, Tuple{Tb,Tb}, typeof(metaval)}(fill(Tb(-Inf), n), fill(Tb(Inf), n), splits, metaval)
end

Base.ndims(world::World) = length(world.position)
baseposition(world::World{T}) where T = baseposition(T, world.position)

Base.in(position::AbstractVector, world::World) = all(t->((x, l, u) = t; l <= x <= u), zip(position, world.lower, world.upper))

## Box

# The tree structure was designed to be non-redundant: there is only
# one way to get any value of interest.  We could copy some data and
# reduce traversal of the tree, but at the risk of values going
# out-of-sync.
# The exception is the world structure, to which a reference is held
# by every box. Without introducing type-instability there may not be
# a way around this.

# For the types below,
# - p is the degree of splitting (the number of dimensions that we split simultaneously)
# - L = 2^p-1 is the number of new (non-self) children created when splitting a box
# - T is the type used for spatial positions
# - M is the type used for the metadata associated with each n-dimensional point
# - B is a parameter representing Box before that type exists (avoids circular definitions)

struct Children{M,B,L}
    children::NTuple{L,B} # the non-self child nodes
    metas::NTuple{L,M}    # metadata for the evaluation points corresponding to new children
end

"""
    struct Split
        dims       # a tuple like (3,5) holding the splitting dimensions
        xs         # a tuple like (1.2,-15.3), the new positions along the splitting dimensions
        self       # the child-box with the same position as the parent box
        others     # the other children and metadata of the parent box
    end
"""
struct Split{p,T,M,B,L}
    dims::NTuple{p,Int}
    xs::NTuple{p,T}
    self::B
    others::Children{M,B,L}
end
Split{p,T}(dims::NTuple{p,Integer},
           xs::NTuple{p,Any},
           self::B,
           others::Children{M,B,L}) where {T,p,M,B,L} =
   Split{p,T,M,B,L}(dims, xs, self, others)

boxtype(::Type{Split{p,T,M,B,L}}) where {p,T,M,B,L} = B
function Base.show(io::IO, split::Split)
    print(io, "Split(")
    show(io, split.self)
    print(io, " along $(split.dims) at $(split.xs))")
end

mutable struct Box{p,T,M,L,P,MB}
    world::World{T,P,M}      # the overall problem domain
    parent::Box{p,T,M,L,P,MB}   # the node above this one
    childindex::UInt8        # which of parent's children is this? (0=self)
    split::Split{p,T,M,Box{p,T,M,L,P,MB},L}  # undefined if this box is a leaf node
    metabox::MB

    function Box{p,T,M,L,P,MB}(world::World{T,P,M}, metabox) where {p,T,M,L,P,MB}
        p > 8 && error("degree must be less than or equal to 8 (change childindex::UInt8 for p>8)")
        root = new{p,T,M,L,P,MB}(world)
        root.parent = root
        root.childindex = 0
        root.metabox = metabox
        return root
    end

    function Box{p,T,M,L,P,MB}(parent::Box{p,T,M,L,P,MB}, splitdims::NTuple{p,Integer},
                               xs::NTuple{p,Real}, metas::NTuple{L,M}) where {p,T,M,L,P,MB}
        function box(parent::Box{p,T,M,L,P,MB}, childindex::Integer) where {p,T,M,L,P,MB}
            return new{p,T,M,L,P,MB}(parent.world, parent, childindex)
        end

        @noinline throw0(sd, mxd) = error("got split along dimension $sd, max allowed is $mxd")
        @noinline throw1(x, bb, sd) = error("out-of-bounds evaluation position $x along dimension $sd (bounds are $bb)")
        @noinline throw2(x, bb, u, sd) = error("cannot split at the upper edge $x along dimension $sd because the boxbounds $bb are not at the upper world edge $u")
        @noinline throw3(x, sd) = error("position $x is identical to the parent along $sd")

        @assert(isleaf(parent))
        # We allow fictive dimensions when needed to "round out"
        # splits. For example, if n is odd but p = 2, we allow a
        # dimension n+1 so that (n+1)/2 pairs can cover all dimensions.
        n = ndims(parent)
        maxdim = ceil(Int, n/p)*p
        for (splitdim,x) in zip(splitdims, xs)
            0 < splitdim <= maxdim || throw0(splitdim, maxdim)
            splitdim > n && continue
            bb = boxbounds(parent, splitdim)
            if (x < bb[1]) | (x > bb[2])
                throw1(x, bb, splitdim)
            elseif x == bb[2] && bb[2] != (u = parent.world.upper[splitdim])
                # To ensure that each point corresponds to a unique leaf
                # (needed for proper behavior of `find_leaf_at`), we
                # essentially define the boxbounds as the half-open
                # interval [lower, upper).  Therefore we have to prevent
                # splitting with an evaluation point at the upper edge of
                # box. (This forces the user to split the other box that
                # shares this edge.) The only exception is for boxes that
                # extend to the upper edge of the world bounds, for which
                # the bounds are defined as the closed interval [lower, upper].
                throw2(x, bb, u, splitdim)
            end
            x != position(parent, splitdim) || throw3(x, splitdim)
        end
        local self, others
        let parent = parent
            self, others = box(parent, 0), ntuple(i->box(parent, i), Val(L))
        end
        parent.split = Split{p,T}(splitdims, xs, self, Children(others, metas))
        return others
    end
end
Box{p}(world::World{T,P,M}, metabox=nothing) where {p,T,M,P} = Box{p,T,M,calcL(Val(p)),P,typeof(metabox)}(world, metabox)

"""
    root = Box{p}(world::World)

Create the root box of a CSp-tree for domain `world`.

    box = Box(parent::Box, splitdims::(Integer...), xs::(Real...), metas::(Any...))

Split a `parent` box along dimensions `splitdims`, where the length of
`splitdims` and `xs` is equal to `p`. The `2^p - 1` new evaluation
points have metadata `metas`, stored in "column-major" order
corresponding to positions

    (xs[1], xp[2], ...)  # (xp[1], xp[2], ...), the "self" box, is omitted
    (xp[1], xs[2], ...)
    (xs[1], xs[2], ...)

where `xp` is `position(parent)[splitdims]`.  For dimensions not
listed in `splitdims`, the positions of these children are all
identical to the parent evaluation point.
"""
Box(parent::Box{p,T,M,L,P,MB}, splitdims::NTuple{p,Integer}, xs::NTuple{p,Real}, metas::NTuple{L,Any}) where {p,T,M,L,P,MB} =
    Box{p,T,M,L,P,MB}(parent, splitdims, xs, metas)

Box(parent::Box{1,T,M,1}, splitdim::Integer, x::Real, meta) where {T,M} =
    Box(parent, (splitdim,), (x,), (meta,))

# You can supply fewer than p dimensions, in which case we fill in
# with fictive dimensions
function Box(parent::Box{p,T,M,L,P,MB}, splitdims::NTuple{k,Integer}, xs::NTuple{k,Real}, metas) where {p,T,M,k,L,P,MB}
    0 < k <= p || throw(DimensionMismatch("got $k dimensions, max allowed is $p"))
    nmeta = 2^k-1
    length(metas) == nmeta || error("got $(length(metas)) metadatas for $k split dimensions, need $nmeta")
    n = ndims(parent)
    parentmeta = meta(parent)
    local splitdimspad, xspad, metaspad
    let n = n, nmeta = nmeta, parentmeta = parentmeta   # julia issue 15276
        splitdimspad = ntuple(Val(p)) do d
            d <= k ? Int(splitdims[d]) : n + d - k
        end
        xspad = ntuple(Val(p)) do d
            d <= k ? T(xs[d]) : zero(T)
        end
        metaspad = ntuple(Val(L)) do d
            i = d & nmeta
            return i == 0 ? parentmeta : metas[i]
        end
    end
    return Box{p,T,M,L,P,MB}(parent, splitdimspad, xspad, metaspad)
end

"""
    isroot(box)

Return `true` if `box` is the root-node of a CSp-tree.
"""
isroot(box::Box) = box.parent == box
"""
    isleaf(box)

Return `true` if `box` has no children, and hence is a leaf-node of a CSp-tree.
"""
isleaf(box::Box) = !isdefined(box, :split)
isself(box::Box)  = !isroot(box) && box.parent.split.self == box
isother(box::Box) = !isroot(box) && box.parent.split.self != box
Base.eltype(::Type{B}) where B<:Box = B
Base.eltype(box::Box) = eltype(typeof(box))
Base.ndims(box::Box) = length(box.world.lower)
degree(::Type{B}) where B<:Box{p} where p = p
degree(box::Box) = degree(typeof(box))
boxcoordtype(::Type{B}) where B<:Box{p,T} where {p,T} = T
boxcoordtype(box::Box) = boxcoordtype(typeof(box))
maxchildren(::Type{Box{p,T,M,L,P,MB}}) where {p,T,M,L,P,MB} = L+1
maxchildren(::Type{B}) where B<:Box{p} where p = 1<<p  # for partial type like Box{2}
maxchildren(box::Box) = maxchildren(typeof(box))
isnonleaf(box) = !isleaf(box)

function Base.show(io::IO, box::Box)
    print(io, "Box")
    show(IOContext(io, :compact=>true), meta(box))
    print(io, "@", position(box))
end

# Compute 2^p-1 in a way that makes it inferrable
@generated function calcL(::Val{p}) where p
    return 1<<p - 1
end

## AbstractTrees interface

struct ChildrenIterator{B<:Box}
    box::B
end

AbstractTrees.children(box::Box) = ChildrenIterator(box)

function Base.iterate(iter::ChildrenIterator)
    (isleaf(iter.box) | maxchildren(iter.box) == 0) && return nothing
    split = iter.box.split
    return split.self, 1
end
function Base.iterate(iter::ChildrenIterator, s::Int)
    (isleaf(iter.box) | maxchildren(iter.box) == s) && return nothing
    split = iter.box.split
    return split.others.children[s], s+1
end
Base.length(iter::ChildrenIterator) = maxchildren(iter.box)
Base.keys(iter::ChildrenIterator) = isleaf(iter.box) ? (1:0) : (0:maxchildren(iter.box)-1)

AbstractTrees.printnode(io::IO, box::Box) = isleaf(box) ? print(io, box) : print(io, box.split.dims)

## Other iteration.

# The iterators in AbstractTrees are not sufficiently high-performance
# for our needs. Moreover we need customized visitation patterns.

Base.IteratorSize(::Type{<:Box}) = Base.SizeUnknown()

# Iterator API: support boxtype and getbox
boxtype(::Type{B}) where B<:Box = B
getbox(box::Box) = box

# state type for most iterators
struct VisitorState{B<:Box}
    box::B
    childindex::Int
end

abstract type CSpTreeIterator end

Base.IteratorSize(::Type{<:CSpTreeIterator}) = Base.SizeUnknown()

struct SplitIterator{B<:Box} <: CSpTreeIterator
    base::B
end

Base.eltype(::Type{SplitIterator{B}}) where B<:Box{p,T,M,L} where {p,T,M,L} =
    Split{p,T,M,B,L}

maxchildren(splits::SplitIterator) = maxchildren(splits.base)

# ClimbingState supports `splits(box)` and is designed to descend and then
# visit siblings as it ascends.
struct ClimbingState{B<:Box,BI,BS}
    box::B                  # current top-level node
    visited::Bool           # true if box.split has already been returned
    skipchildindex::UInt8   # index of box's branch from which we climbed
    childindex::Int         # index of box's branch we're currently exploring
    branchiter::BI          # iterator for the current branch
    branchstate::BS         # state for the current branch
end

# Replace the branch state
ClimbingState(cstate::ClimbingState{B,BI,BS}, vstate::BS) where {B,BI,BS} =
    ClimbingState(cstate.box, cstate.visited, cstate.skipchildindex,
                  cstate.childindex, cstate.branchiter, vstate)

struct ChainIterator{B<:Box} <: CSpTreeIterator
    base::B
end
Base.eltype(::Type{ChainIterator{B}}) where B<:Box{p,T,M,L} where {p,T,M,L} =
    Split{p,T,M,B,L}

## Incremental Gaussian Elimination

struct IGE{T}
    coefs::Matrix{T}
    rhs::Vector{T}
    rowtmp::Vector{T}
end

function IGE{T}(n::Integer) where T
    coefs = fill(zero(T), n, n)
    rhs = fill(zero(T), n)
    rowtmp = Vector{T}(undef, n)
    return IGE{T}(coefs, rhs, rowtmp)
end

## SymmetricArray

# A "utility type" (not really core to this package) to make it easy
# to work with higher-order polynomial models.
struct SymmetricArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::A
end
# SymmetricArray(A::AbstractArray{T,N}) where {T,N} =
#     SymmetricArray{T,N,typeof(A)}(A)

Base.size(S::SymmetricArray) = size(S.data)

function Base.getindex(S::SymmetricArray{T,N}, I::Vararg{Int,N}) where {T,N}
    J = CoordinateSplittingPTrees.tuplesort(I)
    return S.data[J...]
end

function Base.setindex!(S::SymmetricArray{T,N}, val, I::Vararg{Int,N}) where {T,N}
    J = tuplesort(I)
    S.data[J...] = val
    return val
end

tuplesort(dims::Tuple{})        = dims
tuplesort(dims::Tuple{Int})     = dims
tuplesort(dims::Tuple{Int,Int}) = dims[1] > dims[2] ? (dims[2], dims[1]) : dims
tuplesort(dims::NTuple{N,Int}) where N = (sort([dims...])...,)

# Convenience functions
mtrx(S::SymmetricArray{T,2}) where T = Symmetric(S.data, :L)
LinearAlgebra.eigen(S::SymmetricArray{T,2}) where T = LinearAlgebra.eigen(mtrx(S))

const SymmetricMatrix{T} = Union{Symmetric{T}, SymTridiagonal{T}}


text:
using ACE1, JuLIP, Test, IPFitting
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB, observation, vec_obs
using LinearAlgebra: qr, norm

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Cu)
calc = let r0=r0
   LennardJones(r0=r0) * C2Shift(2.5*r0)
end
data = generate_data(:Cu, 3, 0.33*r0, 50, calc)
rcut2 = cutoff(calc)

B2 = rpi_basis(species = :Cu, N = 2, maxdeg = 6,
   r0 = r0, rcut = cutoff(calc), rin = 0.5*r0,
   pin = 2)

# generate a database so that we generate the row indices
db = LsqDB("", B2, data)

## generate a LSQ system by hand
atlen = length(data[1])
len = 1 + length(observation(data[1], "F"))
Y_man = zeros(length(data) * len)
Ψ_man = zeros(length(data) * len, length(B2))
for (okey, evalb, w) in zip(["E", "F"], (energy, forces), (1.345/sqrt(atlen), 2.987))
   for d in data
      o = observation(d, okey)
      irows = IPFitting.DB.matrows(d, okey)
      leno = length(o)
      Y_man[irows] = w * o
      #for colidx = 1:length(B2)
      Ψ_man[irows, :] = w * vec_obs(okey, evalb(B2, d.at))
      #end
   end
end

## test that Ψ_man, F_man are equivalent to a system generated from a db

Ψ, Y = IPFitting.Lsq.get_lsq_system( db,
               Vref = OneBody(:Cu => 0.0),
               weights = Dict("default" => Dict("E" => 1.345, "F" => 2.987)),
            )

##
print_tf(@test Y ≈ Y_man)
print_tf(@test Ψ ≈ Ψ_man)
println()


text:
module Cuas
using ExternalBallistics
using HCubature
using ReferenceFrameRotations
# Write your package code here.
include("types.jl")
include("roundPerRound.jl")
include("generate.jl")
include("rotations.jl")
#include("phit.jl")
#using .phit

export pHit, RectErrorB, createTourelle, createProjectile, targetPos, muzzlePos, muzzleVel, wind, Tourelle,
Hull, Turret, Canon, Sight, Tank, SSHP

end


text:
using Test
using BSDiff
using ArgTools
using Pkg.Artifacts

import bsdiff_classic_jll
import bsdiff_endsley_jll
import zrl_jll

const test_data = artifact"test_data"
const FORMATS = sort!(collect(keys(BSDiff.FORMATS)))

@testset "BSDiff" begin
    @testset "API coverage" begin
        # create new and old files
        dir = mktempdir()
        old_file = joinpath(dir, "old")
        new_file = joinpath(dir, "new")
        write(old_file, "Goodbye, world.")
        write(new_file, "Hello, world!")

        # generate reference index
        index_file = joinpath(dir, "index")
        index_file = bsindex(old_file)

        @testset "bsindex" begin
            arg_readers(old_file) do old
                index_file′ = @arg_test old bsindex(old)
                @test read(index_file′) == read(index_file)

                arg_writers() do index_file′, index
                    @arg_test old index begin
                        @test index == bsindex(old, index)
                    end
                    @test read(index_file′) == read(index_file)
                end
            end
        end

        for format in (nothing, :classic, :endsley)
            fmt = format == nothing ? [] : [:format => format]

            # generate and test reference patch first
            patch_file = bsdiff(old_file, new_file; fmt...)
            @testset "reference patch" begin
                new_file′ = bspatch(old_file, patch_file; fmt...)
                @test read(new_file′) == read(new_file)
                rm(new_file′)
            end

            @testset "bsdiff" begin
                arg_readers(old_file) do old
                    arg_readers(new_file) do new
                        patch_file′ = @arg_test old new begin
                            bsdiff(old, new; fmt...)
                        end
                        @test read(patch_file′) == read(patch_file)
                        rm(patch_file′)

                        arg_readers(index_file) do index
                            patch_file′ = @arg_test old index new begin
                                bsdiff((old, index), new; fmt...)
                            end
                            @test read(patch_file′) == read(patch_file)
                            rm(patch_file′)
                        end

                        arg_writers() do patch_file′, patch
                            @arg_test old new patch begin
                                @test patch == bsdiff(old, new, patch; fmt...)
                            end
                            @test read(patch_file′) == read(patch_file)
                        end
                    end
                end
            end

            @testset "bspatch" begin
                arg_readers(old_file) do old
                    arg_readers(patch_file) do patch
                        # auto-detected format
                        new_file′ = @arg_test old patch begin
                            bspatch(old, patch)
                        end
                        @test read(new_file′) == read(new_file)
                        rm(new_file′)

                        arg_writers() do new_file′, new
                            @arg_test old new patch begin
                                @test new == bspatch(old, new, patch)
                            end
                            @test read(new_file′) == read(new_file)
                        end
                        format === nothing &&
                            return # tests below are the same

                        # explicit format
                        new_file′ = @arg_test old patch begin
                            bspatch(old, patch; fmt...) # explicit format
                        end
                        @test read(new_file′) == read(new_file)
                        rm(new_file′)

                        arg_writers() do new_file′, new
                            @arg_test old new patch begin
                                @test new == bspatch(old, new, patch; fmt...)
                            end
                            @test read(new_file′) == read(new_file)
                        end
                    end
                end
            end
        end
        rm(dir, recursive=true, force=true)
    end

    @testset "registry data" begin
        registry_data = joinpath(test_data, "registry")
        old = joinpath(registry_data, "before.tar")
        new = joinpath(registry_data, "after.tar")
        ref = joinpath(registry_data, "reference.diff")
        old_data = read(old)
        new_data = read(new)
        @testset "hi-level API (w/ timing)" for format in FORMATS
            println("[ raw data ]")
            @show format
            index = bsindex(old)
            patch = @time bsdiff((old, index), new, format = format)
            patch = @time bsdiff((old, index), new, format = format)
            patch = @time bsdiff((old, index), new, format = format)
            new′ = bspatch(old, patch)
            @test read(new) == read(new′)
            @show filesize(patch)
        end

        @testset "ZRL data (w/timing)" for format in FORMATS
            println("[ ZRL data ]")
            @show format
            old_zrl = tempname()
            new_zrl = tempname()
            zrl_jll.zrle() do zrle
                run(pipeline(`$zrle $old`, old_zrl))
                run(pipeline(`$zrle $new`, new_zrl))
            end
            index = bsindex(old_zrl)
            patch = @time bsdiff((old_zrl, index), new_zrl, format = format)
            patch = @time bsdiff((old_zrl, index), new_zrl, format = format)
            patch = @time bsdiff((old_zrl, index), new_zrl, format = format)
            new_zrl′ = bspatch(old_zrl, patch)
            @test read(new_zrl) == read(new_zrl′)
            @show filesize(patch)
            rm(old_zrl)
            rm(new_zrl)
        end

        @testset "low-level API" begin
            # test that diff is identical to reference diff
            index = BSDiff.generate_index(old_data)
            diff = sprint() do io
                patch = BSDiff.EndsleyPatch(io, length(new_data))
                BSDiff.generate_patch(patch, old_data, new_data, index)
            end |> codeunits
            @test read(ref) == diff
            # test that applying reference patch to old produces new
            new_data′ = open(ref) do io
                patch = BSDiff.EndsleyPatch(io, length(new_data))
                sprint() do new_io
                    BSDiff.apply_patch(patch, old_data, new_io)
                end |> codeunits
            end
            @test new_data == new_data′
        end

        if !Sys.iswindows() # bsdiff JLLs don't compile on Windows
            for (format, jll) in [
                    (:classic, bsdiff_classic_jll),
                    (:endsley, bsdiff_endsley_jll),
                ]
                @testset "high-level API" begin
                    # test that bspatch command accepts patches we generate
                    patch = bsdiff(old, new, format = format)
                    new′ = tempname()
                    jll.bspatch() do bspatch
                        run(`$bspatch $old $new′ $patch`)
                    end
                    @test new_data == read(new′)
                    rm(new′)
                    # test that we accept patches generated by bsdiff command
                    patch = tempname()
                    jll.bsdiff() do bsdiff
                        run(`$bsdiff $old $new $patch`)
                    end
                    new′ = bspatch(old, patch, format = format)
                    @test new_data == read(new′)
                    new′ = bspatch(old, patch) # format auto-detected
                    @test new_data == read(new′)
                end
            end
        end
    end
end


text:
"""
Implementation of local response normalization.
"""
module LRNs
using Flux
import Flux: params

export LRN

"""
    LRN{T,I}

Local response normalization layer. Input `i` is processed according to
out(x,y,f,b) = x / [c + α * sum( i(x-k÷2:x+k÷2, y-k÷2:y+k÷2, f-n÷2:f+n÷2, b)^2 )]^β

Todo: β is currently ignored (always set to 0.5)

# Arguments
- `c::T`: additive constant for normalization sum. Default 1.0
- `α::T`: multiplicative scaling constant for normalization sum. Default 1.0
- `β::T`: power scaling constant for normalization sum. Default 0.5
- `size::Tuple{I, I}`: size of normalization field in image space. Default (1, 1)
- `depth::I`: number of neighboring feature maps that go into normalization. Default 5.
"""
struct LRN{T,I,A}
    c::T
    α::T
    β::T
    size::Tuple{I, I}
    depth::I
    kernel::A
end

# convenience constructor that generates the appropriate summation kernel
function LRN(c::T, α::T, β::T, size::Tuple{I, I}, depth::I, features::I) where {T, I}
    # using convolution to do the sum, like https://github.com/FluxML/Flux.jl/pull/720
    # create kernel:
    kernel = zeros(T, size[1], size[2], features, features)
    for i in 1:features
        lower = max(i-depth÷2, 1)
        upper = min(i+depth÷2, features)
        kernel[:, :, lower:upper, i] .= 1
    end
    LRN{T, I, typeof(kernel)}(c, α, β, size, depth, kernel)
end

# convenience constructor with default arguments
function LRN(features; c=1.0, α=1.0, β=0.5, size=(1, 1), depth=5)
    return LRN(c, α, β, size, depth, features)
end # function LRN

"""
    (l::LRN)(i)

Applies a local response normalization layer according to:
out(x,y,f,b) = x ./ [c + α * sum( i(x, y, f-n÷2:f+n÷2, b)^2 )]^(β)

Todo: β is currently ignored (always set to 0.5)
"""
function (l::LRN)(x)
    norm = l.c .+ l.α .* conv(x .^ 2, l.kernel)
    norm = norm .^ l.β
    return x ./ norm
end # function (l::LRN)

Flux.params(l::LRN) = nothing
Flux.@functor LRN

end # module LRNs


text:
module ExceptionRequest

export StatusError

import ..HTTP
using ..Messages: iserror

export exceptionlayer

"""
    exceptionlayer(stream) -> HTTP.Response

Throw a `StatusError` if the request returns an error response status.
"""
function exceptionlayer(handler)
    return function(stream; status_exception::Bool=true, kw...)
        res = handler(stream; kw...)
        if status_exception && iserror(res)
            throw(StatusError(res.status, res.request.method, res.request.target, res))
        else
            return res
        end
    end
end

"""
    StatusError <: Exception

The `Response` has a `4xx`, `5xx` or unrecognised status code.

Fields:
 - `status::Int16`, the response status code.
 - `response` the [`HTTP.Response`](@ref)
"""
struct StatusError <: Exception
    status::Int16
    method::String
    target::String
    response::HTTP.Response
end

# for backwards compatibility
StatusError(status, response::HTTP.Response) = StatusError(status, response.request.method, response.request.target, response)

end # module ExceptionRequest


text:
"""
    spy(A; kwargs...)

Description
============

Plots the sparsity pattern for the given matrix `A`. This means
that a scatterplot that resembles the matrix is drawn, in which
only the pixel for non-zero elements of the matrix are set.

If the parameters `width` and `height` are not explicitly
specified, then the function will attempt to preserve the aspect
ratio of the matrix, while also attempting to fit the resulting
plot withing the bounding box specified by `maxwidth` and
`maxheight`

Usage
======

    spy(A; maxwidth = 70, maxheight = 40, title = "Sparsity Pattern", labels = true, border = :solid, margin = 3, padding = 1, color = :auto, out_stream::Union{Nothing,IO} = nothing,width = 0, height = 0, canvas = BrailleCanvas, grid = true)

Arguments
==========

- **`A`** : The matrix of interest for which non-zero elements
  should be drawn

- **`maxwidth`** : Maximum number of characters per row that
  should be used for plotting

- **`maxheight`** : Maximum number of character rows that should
  be used for plotting.

- **`title`** : Text to display on the top of the plot.

- **`labels`** : Boolean. Can be used to hide the labels by
  setting `labels = false`.

- **`border`** : The style of the bounding box of the plot.
  Supports `:corners`, `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.
- **`margin`** : Number of empty characters to the left of the
  whole plot.

- **`padding`** : Space of the left and right of the plot between
  the labels and the canvas.

- **`color`** : Color of the drawing.
  Can be any of `:green`, `:blue`, `:red`, `:yellow`, `:cyan`,
  `:magenta`, `:white`, `:normal`.

- **`width`** : Exact number of characters per row that should be
  used for plotting. `0` stands for automatic.

- **`height`** : Exact number of character rows that should be
  used for plotting. `0` stands for automatic.

- **`canvas`** : The type of canvas that should be used for drawing.

Returns
========

A plot object of type `Plot{T<:Canvas}`

Author(s)
==========

- Dominique Orban (Github: https://github.com/dpo)
- Christof Stocker (Github: https://github.com/Evizero)
- Jake Bolewski (Github: https://github.com/jakebolewski)

Examples
=========

```julia-repl
julia> spy(sprandn(50, 120, .05))
                           Sparsity Pattern
      ┌────────────────────────────────────────────────────────────┐
    1 │⢀⠀⠒⠀⢄⠂⠀⡀⠀⠀⠀⠰⠀⠐⠀⠀⠀⠀⠂⠀⠀⠀⠠⠰⠀⠀⠀⡀⠄⠀⠠⠀⠀⠀⠀⠀⡀⠀⠀⠀⡀⠐⠀⠐⠀⠀⡀⠀⠀⠀⠀⠀⠀⠂⡀⠀⡀⠠⠀⠀│ > 0
      │⠀⠀⠀⠀⠀⢀⢀⠀⠀⠀⠀⠐⠄⠀⢀⠂⠀⠠⠀⠉⠀⠀⠀⠀⢐⠀⠀⠀⠀⠀⠠⠐⠄⠄⠀⠀⠀⠐⠐⠈⢁⠀⠀⠀⠀⠀⠀⠀⠠⠀⠐⡐⠀⠄⠀⠀⠈⠀⠀⠠│ < 0
      │⡀⠀⠀⠀⡀⠊⠀⠀⠀⢀⢀⠀⠄⠀⠂⡠⠂⠐⠐⠀⠀⠄⠀⠀⠀⠄⠀⠀⠀⠀⡀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠠⠀⠂⠀⠠⠀⠀⠀⠆⠄⠈⠀⢀⠀⠀⠀⠐⠄⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⢠⠀⢀⠀⠀⠀⢀⠀⠄⠀⠀⠀⠂⠀⠀⠀⠇⠀⢀⠀⠀⠁⠁⠀⠄⠄⠀⠀⠁⠁⠀⡈⠀⠁⠄⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡐⢀⠀⠀⠀⠐⠂⠁⢀⠠⠀⠀⠀⠀⠀⠁⠀⠀⠀⠠⠀⠀⠀⠀⢀⠀⢄⠀⠀⠀⠐⠀⠀⠈⠀⠀⠐⠀⠄⠂⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠠⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠠⠈⠀⠀⠀⠀⢀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠐⠀⠉⠐⠀⠀⠁⠀⠀⠀⠀⠀⠀⠌⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠌⠀⠀⠄⠀⠀⡀⠀⠀⠀⠀⠀⢀⠈⠄⠀⡀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀│
      │⠂⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⡨⠀⠠⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠈⡠⠀⡀⠀⠀⠁⠀⠁⠠⠀⠀⠀⠀⠀⠀⠀⠄⡀⠀⠄⡠⠀⠀⠀⠀│
      │⠐⠀⠀⠀⠀⠀⠀⠃⠀⠀⠀⠀⠁⠀⡀⠁⠀⢂⠂⠀⠀⠀⠈⠁⠀⠀⠁⠁⠄⠀⠀⠀⠀⠀⡀⠀⠀⠀⠂⠀⠀⠀⠀⠐⠀⠂⠀⠂⠁⠀⠆⠈⠀⠀⠄⠄⠀⠀⠀⠀│
      │⠀⠀⠀⠠⠀⠁⠀⠀⠀⠀⠂⠀⠠⠀⠀⠌⠀⠀⠂⠄⠀⡀⠂⠀⠀⡀⠠⠁⠄⠀⠀⠀⠀⠀⠄⠈⠀⢂⠐⠀⢀⠀⠠⠀⠄⠁⡀⠀⠀⠐⠀⠀⡐⠀⠀⠀⠁⠀⠀⠔│
      │⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠁⠃⡀⡁⠀⠀⠀⢀⠀⠀⠂⠀⠀⠠⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠬⠀⡀│
      │⠐⠀⠀⠀⠀⠀⠀⠀⠢⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠈⠁⠁⠀⠀⠄⠀⠀⠐⠀⠄⠀⠠⠀⠁⡀⠁⠀⠀⠠⠡⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⡀│
   50 │⠀⠀⠀⠀⠁⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠂⠀⠀⠀⠀⠀⠀⠐⠀⠈⠂⠀⠀⠀⠀⠀⠀⠀⠀⠠⠈⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⠀⠀⠂⠁⠀⠀⠀⠂⠀│
      └────────────────────────────────────────────────────────────┘
      1                                                          120
                               nz = 275
```

See also
=========

[`Plot`](@ref), [`scatterplot`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
spy(A::AbstractMatrix; kwargs...) = spy(size(A)..., _findnz(A)...; kwargs...)

function spy(
    nrow::Int,
    ncol::Int,
    rows::AbstractArray{<:Integer},
    cols::AbstractArray{<:Integer},
    vals::AbstractArray;
    maxwidth::Int  = 0,
    maxheight::Int = 0,
    title = "Sparsity Pattern",
    out_stream::Union{Nothing,IO} = nothing,
    width::Int  = 0,
    height::Int = 0,
    margin::Int  = 3,
    padding::Int = 1,
    color::UserColorType = :auto,
    canvas::Type{T} = BrailleCanvas,
    kw...
) where {T <: Canvas}
    if color == :automatic
        Base.depwarn("`color = :automatic` is deprecated, use `color = :auto` instead", :spy)
        color = :auto
    end
    width, height, _, _ = get_canvas_dimensions_for_matrix(
        canvas, nrow, ncol, maxwidth, maxheight, width, height, margin, padding, out_stream;
        extra_rows = 9, extra_cols = 6
    )
    can = T(width, height, width  = Float64(ncol) + 1, height = Float64(nrow) + 1)
    plot = Plot(can; title = title, margin = margin, padding = padding, kw...)
    height = nrows(plot.graphics)
    width = ncols(plot.graphics)
    plot = if color != :auto
        points!(
            plot, convert(Vector{Float64}, cols),
            nrow + 1 .- convert(Vector{Float64}, rows), color
        )
    else
        pos_idx = vals .> 0
        neg_idx = (!).(pos_idx)
        pos_cols = cols[pos_idx]
        pos_rows = rows[pos_idx]
        neg_cols = cols[neg_idx]
        neg_rows = rows[neg_idx]
        points!(
            plot, convert(Vector{Float64}, pos_cols),
            nrow + 1 .- convert(Vector{Float64}, pos_rows), :red
        )
        points!(
            plot, convert(Vector{AbstractFloat}, neg_cols),
            nrow + 1 .- convert(Vector{Float64}, neg_rows), :blue
        )
        label!(plot, :r, 1, "> 0", :red)
        label!(plot, :r, 2, "< 0", :blue)
    end
    label!(plot, :l, 1, "1", :light_black)
    label!(plot, :l, height, string(nrow), :light_black)
    label!(plot, :bl, "1", :light_black)
    label!(plot, :br, string(ncol), :light_black)
    xlabel!(plot, string("nz = ", length(vals)))
    return plot
end

function _findnz(A::AbstractMatrix)
    I = findall(!iszero, A)
    getindex.(I, 1), getindex.(I, 2), A[I]
end

_findnz(A::AbstractSparseMatrix) = findnz(A)


text:
# Send message to API
function send_message(message, task)
    path = string(
        "/api/v1/organizations/",
        task["organization_id"],
        "/locations/",
        task["location_id"],
        "/roles/",
        task["role_id"],
        "/schedules/",
        task["schedule_id"],
        "/messages/",
    )
    results = Requests.post(
        build_uri(config, path),
        json={"message" => message}
    )

    if results.status != 201
        code = results.status
        reason = results.data
        Logging.err( "Unable to create message $message - status %code - response $reason")
        error()
    end
    Logging.info( "Sent message - $message")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
using Test
using JSON
using MianalyzerBackend
const KMeansController = MianalyzerBackend.KMeansController

include("../TestSupport.jl")

@testset "KMeanController" begin
  @testset ".call" begin
    request = Dict([
                    ("records", [[1.0, 1.0], [2.0, 2.0]]),
                    ("cluster_num", 2)
                   ])

    actual = JSON.parse(IOBuffer(KMeansController.call(request)))

    @test actual["clusters"] == [[[1.0, 1.0]], [[2.0, 2.0]]]
    @test actual["centers"] == [[1.0, 1.0], [2.0, 2.0]]
  end
end


text:
#HullReach
using HullReachNV
using Test
using LazySets
import HullReachNV: ReLU

nnet = read_nnet("nnet/toy_nnet.nnet", last_layer_activation = ReLU())

solver = HullReach(0.001, true)

center = [1.0, 1.0]
radius = [1.0, 1.0]

in_hyper = Hyperrectangle(center, radius)

lower = [0.0, 0.0]
upper = [1.6, 2.1]

out_hyper = Hyperrectangle(low=lower, high=upper)

problem = Problem(nnet, in_hyper, out_hyper)

print("HullReach - test")
timed_result =@timed solve(solver, problem)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1].status)
print("\n")


text:
@testset "Covariances" begin
  comps = rand(Composition{5}, 100)
  data  = reduce(hcat, components.(comps))
  table = (; zip([:N₂,:O₂,:CO₂,:S,:C], eachrow(data))...)

  Τ = variation(table)
  Σ = alrcov(table)
  Γ = clrcov(table)
  A = lrarray(table)

  # Testing matrix forms transformations (Aitchison 1986, pages 82 and 83)
  ## Test matrix transformation from Τ to Σ (4.25)
  @test Σ ≈ -0.5 * F * Τ * F'

  ## Test matrix transformation from Τ to Γ (4.26)
  @test Γ ≈ -0.5 * G * Τ * G

  ## Test matrix transformation from Σ to Τ (4.27)
  Γ₀ = F' * (inv(H) * Σ * inv(H)) * F
  @test Τ ≈ J * Diagonal(Γ₀) + Diagonal(Γ₀) * J - 2 * Γ₀

  ## Test matrix transformation from Σ to Γ (4.28)
  @test Γ ≈ F' * (inv(H) * Σ * inv(H)) * F

  ## Test matrix transformation from Γ to Τ (4.29)
  @test Τ ≈ J * Diagonal(Γ) + Diagonal(Γ) * J - 2 * Γ

  ## Test matrix transformation from Γ to Σ (4.30)
  @test Σ ≈ F * Γ * F'
end


text:
"""
    Zeros(::Type{T}, shape...) where T
return an all-zero-elements-array of type T which has shape `shape...`

# Example
    julia> Zeros(Array{Float64}, 2, 5)
    2×5 Array{Float64,2}:
     0.0  0.0  0.0  0.0  0.0
     0.0  0.0  0.0  0.0  0.0
 """
function Zeros(::Type{T}, shape...) where T
    return fill!(T(undef, shape...), 0.0)
end

"""
    Ones(::Type{T}, shape...) where T
return an all-one-elements-array of type T which has shape `shape...`

# Example
    julia> 7Ones(Array{Float64}, 2, 5)
    2×5 Array{Float64,2}:
     7.0  7.0  7.0  7.0  7.0
     7.0  7.0  7.0  7.0  7.0
 """
function Ones(::Type{T}, shape...) where T
    return fill!(T(undef, shape...), 1.0)
end

"""
    mutable struct Variable{T}
# Fields
+ `value::T`                : value in forward
+ `delta::Union{Nothing,T}` : gradients collected in backprop
+ `shape::Tuple`            : shape of `value`
+ `keepsgrad::Bool`         : whether keeps grad after backprop
+ `isleaf::Bool`            : whether leaf node
"""
mutable struct Variable{T}
    value::T
    delta::Union{Nothing,T}
    shape::Tuple
    backprop::Bool      # 是否反向传播，依赖的子节点都跟其一样
    keepsgrad::Bool     # 是否保留梯度，可能用于训练或者其他
    isleaf::Bool        # 可训练的，代表可学习的参数
    function Variable{T}(x, backprop::Bool=true,
                            keepsgrad::Bool=false,
                            isleaf::Bool=false) where T
        s = size(x)
        δ = nothing
        new{T}(x, δ, s, backprop, keepsgrad, isleaf)
    end
end

# Convenient abstract-type-specilized constructors for data on GPU/CPU/xPU etc....
function Variable(x; backprop::Bool=true,
                     keepsgrad::Bool=false,
                     type::Type=Array{Float32})
    isleaf = true    # any user defined Variable is a leaf Variable
    return Variable{type}(x, backprop, keepsgrad, isleaf)
end

# pretty printing
function Base.show(io::IO, var::Variable{T}) where T
    if  var.isleaf println(cyan("\n≡≡≡ Leaf Variable ≡≡≡")) end
    if !var.isleaf println(cyan("\n≡≡≡ None Leaf Variable ≡≡≡")) end

    print(blue("\nvalue is "))
    display(var.value)
    print(green("\ndelta is "))
    display(var.delta)
end

function clone(x::Variable; type::Type=Array{Float32})
    return Variable{type}(type(x.value), x.backprop, x.keepsgrad, x.isleaf)
end

function zeroDelta(var::Variable{T}) where T
    # 要切断某些反向传播路径的时候将其初始化为零
    if var.delta === nothing
        var.delta = Zeros(T, var.shape);
    end
end

function need2computeδ!(var::Variable{T}) where T
    # 需要反向传播的时候就需要初始化
    if !(var.isleaf && !var.keepsgrad)
        if var.delta===nothing
            var.delta = Zeros(T, var.shape);
        end
        return true
    else
        return false
    end
end

function ifNotKeepδThenFreeδ!(var::Variable{T}) where T
    if !var.keepsgrad
        var.delta = nothing
    end
end

Base.sizeof(x::Variable)         =  sizeof(x.value)
Base.size(x::Variable)           =    size(x.value)
Base.size(x::Variable, dim::Int) =    size(x.value, dim)
Base.ndims(x::Variable)          =   ndims(x.value)
Base.length(x::Variable)         =  length(x.value)
Base.strides(x::Variable)        = strides(x.value)
Base.eltype(x::Variable)         =  eltype(x.value)
Base.similar(x::Variable{T})  where T = Variable{T}( similar(x.value), x.backprop, x.keepsgrad, x.isleaf)
Base.copy(x::Variable{T})     where T = Variable{T}(    copy(x.value), x.backprop, x.keepsgrad, x.isleaf)
Base.deepcopy(x::Variable{T}) where T = Variable{T}(deepcopy(x.value), x.backprop, x.keepsgrad, x.isleaf)

Base.setindex!(x::Variable, v::Number,        k...) = (x.value[k...] .= v)
Base.setindex!(x::Variable, v::AbstractArray, k...) = (x.value[k...]  = v)

function Base.getindex(x::Variable{T}, k...) where T
    y = Variable{T}(x.value[k...], x.backprop, x.keepsgrad, x.isleaf)
    if x.backprop
        function getindexBackward()
            if need2computeδ!(x)
                x.delta[k...] .+= y.delta
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, getindexBackward)
    end
    return y
end

function (v::Variable)(i...)
    if v.delta ≠ nothing
        return v.delta[i...]
    else
        return nothing
    end
end

"""
    to(type::Type, var::Variable{T}, show::Bool=false) where T -> Variable{type}
Type conversions between data types in `Variable`. (like from Array{Float64} to CuArray{Float32} )
# Example
    x = Variable(ones(1,4), type=Array{Float64})
    y = to(CuArray{Float16}, x, true)
"""
function to(type::Type, var::Variable{T}, show::Bool=false) where T
    if type <: T || type >: T
        COLOR_TYPE = yellow(T)
        COLOR_DEST = yellow(type)
        @info "you don't need to convert data from $COLOR_TYPE to $COLOR_DEST.\n"
        return var
    else
        if !show
            return Variable{type}(var.value, var.backprop, var.keepsgrad, var.isleaf)
        else
            COLOR_TYPE = green(T)
            COLOR_DEST = green(type)
            @info "converting data from $COLOR_TYPE to $COLOR_DEST.\n"
            return Variable{type}(var.value, var.backprop, var.keepsgrad, var.isleaf)
        end
    end
end

"""
    to(type::Type, vars::Vector{Variable}, show::Bool=false)
Type conversions between data types in `Vector{Variable}`
# Example
    julia> v = Vector{Variable}(undef,2);
    julia> v[1] = Variable(4ones(2,4),type=Array{Float64});
    julia> v[2] = Variable(3ones(2,3),type=Array{Float64});

    julia> typeof(v[1])
    Variable{Array{Float64,N} where N}

    julia> v = to(Array{Float32},v);

    julia> typeof(v[1])
    Variable{Array{Float32,N} where N}
"""
function to(type::Type, vars::Vector{Variable}, show::Bool=false)
    for i = 1:length(vars)
        vars[i] = to(type, vars[i], show);
    end
    return vars
end

function to!(type::Type, vars::Vector{Variable}, show::Bool=false)
    vars = to(type, vars, show)
    return nothing
end

export Variable
export zeroDelta
export Zeros, Ones
export to, to!
export clone
export need2computeδ!
export ifNotKeepδThenFreeδ!
export elsizeof
export value, delta, ᵛ, ᵟ, δ

export XVariable, VarOrNil
const  XVariable = Tuple{Char, Variable}
const  VarOrNil  = Union{Variable, Nothing}

# pretty printing
function Base.show(io::IO, xvar::XVariable)
    c, var = xvar
    if  var.isleaf println(cyan("\n≡≡≡ Leaf Variable ($c) ≡≡≡")) end
    if !var.isleaf println(cyan("\n≡≡≡ None Leaf Variable ≡≡≡")) end

    print(blue("\nvalue is "))
    display(var.value)
    print(green("\ndelta is "))
    display(var.delta)
end

elsizeof(x::Variable) = sizeof(eltype(x))

# lazy showing way of Variable's main vars
@inline ᵛ(x::Variable) = x.value
@inline ᵟ(x::Variable) = x.delta
@inline δ(x::Variable) = x.delta
@inline value(x::Variable) = x.value
@inline delta(x::Variable) = x.delta

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
module SQS

using LightXML
using AWS.Crypto
using Requests
using AWS.AWSEnv
using AWS
using Compat

import AWS: AbstractAWSType, AWSError, AWSResponse, _ename, _parse_from_xml, _prep_params, aws_execute

const SQSError = AWSError
const SQSResponse = AWSResponse

parse_from_xml{T<:AbstractAWSType}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd::Union{LightXML.XMLElement,Void}) = _parse_from_xml(T, types, pd; member_pfx="")
prep_params{T<:AbstractAWSType}(obj::T) = _prep_params(obj; member_pfx="")

sqs_execute(env::AWSEnv, action::String, params::Vector{Tuple}=Vector{Tuple}()) = aws_execute(env, action, "sqs", SQSApiVer, params)
sqs_execute(env::AWSEnv, action::String, req::AbstractAWSType) = sqs_execute(env, action, prep_params(req))
sqs_execute(env::AWSEnv, action::String, req::Void) = sqs_execute(env, action)

include("sqs_types.jl")
include("sqs_operations.jl")
include("sqs_spec.jl")

export SQSError, SQSResponse

end # module SQS


text:
module PoissonTests

using Gridap
using Gridap.Algebra
using Gridap.FESpaces
using GridapDistributed
using PartitionedArrays
using Test

function main(parts)

  output = mkpath(joinpath(@__DIR__,"output"))

  domain = (0,4,0,4)
  cells = (4,4)
  model = CartesianDiscreteModel(parts,domain,cells)

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet",[1,2,3,5,7])
  add_tag_from_tags!(labels,"neumann",[4,6,8])

  Ω = Triangulation(model)
  Γn = Boundary(model,tags="neumann")
  n_Γn = get_normal_vector(Γn)

  k = 2
  u((x,y)) = (x+y)^k
  f(x) = -Δ(u,x)
  g = n_Γn⋅∇(u)

  reffe = ReferenceFE(lagrangian,Float64,k)
  V = TestFESpace(model,reffe,dirichlet_tags="dirichlet")
  U = TrialFESpace(u,V)

  dΩ = Measure(Ω,2*k)
  dΓn = Measure(Γn,2*k)

  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
  l(v) = ∫( v*f )dΩ + ∫( v*g )dΓn
  op = AffineFEOperator(a,l,U,V)

  uh = solve(op)
  eh = u - uh
  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9

  # Now with DG
  h = 2
  γ = 10

  V_dg = FESpace(model,reffe,conformity=:L2)

  Λ = Skeleton(model)
  Γd = Boundary(model,tags="dirichlet")

  dΛ = Measure(Λ,2*k)
  dΓd = Measure(Γd,2*k)

  n_Γd = get_normal_vector(Γd)
  n_Λ = get_normal_vector(Λ)

  a_dg(u,v) =
    ∫( ∇(v)⋅∇(u) )*dΩ +
    ∫( (γ/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u )*dΓd +
    ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) -
       jump(v*n_Λ)⋅mean(∇(u)) -
       mean(∇(v))⋅jump(u*n_Λ) )*dΛ

  l_dg(v) =
    ∫( v*f )*dΩ +
    ∫( v*g )dΓn +
    ∫( (γ/h)*v*u - (n_Γd⋅∇(v))*u )*dΓd

  op = AffineFEOperator(a_dg,l_dg,V_dg,V_dg)
  uh = solve(op)
  eh = u - uh
  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9

end

end # module


text:
abstract Enumerable{T}
Base.eltype{T}(source::Enumerable{T}) = T


text:
# based on https://dlmf.nist.gov/11.2.i
# for |z| << 1

const order_small = 10

_inner_small(n, ν, z) = (0.5z)^(2n) / (gamma(n + 1.5) * gamma(n + ν + 1.5))

_H_small(ν, z) = (0.5z)^(ν + 1) * sum((-1)^n * _inner_small(n, ν, z) for n in 0:order_small)

_L_small(ν, z) = (0.5z)^(ν + 1) * sum(_inner_small(n, ν, z) for n in 0:order_small)


text:
#=
Copyright (c) 2016 Saurav Sachidanand

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
=#

function encode_unsigned_with_type(type_bits::UInt8,
                                   num::Unsigned,
                                   bytes::Array{UInt8, 1})
    if num < SINGLE_BYTE_UINT_PLUS_ONE
        byte_len = 0
        addntl_info = num
    elseif num < UINT8_MAX_PLUS_ONE
        byte_len = sizeof(UInt8)
        addntl_info = ADDNTL_INFO_UINT8
    elseif num < UINT16_MAX_PLUS_ONE
        byte_len = sizeof(UInt16)
        addntl_info = ADDNTL_INFO_UINT16
    elseif num < UINT32_MAX_PLUS_ONE
        byte_len = sizeof(UInt32)
        addntl_info = ADDNTL_INFO_UINT32
    elseif num < UINT64_MAX_PLUS_ONE
        byte_len = sizeof(UInt64)
        addntl_info = ADDNTL_INFO_UINT64
    else
        error("128-bits ints can't be encoded in the CBOR format.")
    end

    push!(bytes, type_bits | addntl_info)

    i = length(bytes) + 1
    for _ in 1:byte_len
        insert!(bytes, i, num & LOWEST_ORDER_BYTE_MASK)
        num >>>= BITS_PER_BYTE
    end
end

# ------- straightforward encoding for a few Julia types

function encode(bool::Bool, bytes::Array{UInt8, 1})
    push!(bytes, CBOR_FALSE_BYTE + bool)
end

function encode(num::Unsigned, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_0, num, bytes)
end

function encode(num::Signed, bytes::Array{UInt8, 1})
    if num < 0
        return encode_unsigned_with_type(TYPE_1, Unsigned(-num - 1), bytes)
    else
        return encode_unsigned_with_type(TYPE_0, Unsigned(num), bytes)
    end
end

function encode(byte_string::AbstractVector{UInt8}, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_2, Unsigned(length(byte_string)), bytes)
    append!(bytes, byte_string)
end
if VERSION < v"0.5.0"
    function encode(string::Union{UTF8String, ASCIIString}, bytes::Array{UInt8, 1})
        encode_unsigned_with_type(TYPE_3, Unsigned(sizeof(string)), bytes)
        append!(bytes, string.data)
    end
else
    function encode(string::String, bytes::Array{UInt8, 1})
        encode_unsigned_with_type(TYPE_3, Unsigned(sizeof(string)), bytes)
        append!(bytes, Vector{UInt8}(string))
    end
end

function encode(list::Union{AbstractVector, Tuple}, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_4, Unsigned(length(list)), bytes)
    for e in list
        encode(e, bytes)
    end
end

function encode(map::AbstractDict, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_5, Unsigned(length(map)), bytes)
    for (key, value) in map
        encode(key, bytes)
        encode(value, bytes)
    end
end

function encode(big_int::BigInt, bytes::Array{UInt8, 1})
    if big_int < 0
        hex_str = hex(-big_int - 1)
        tag = NEG_BIG_INT_TAG
    else
        hex_str = hex(big_int)
        tag = POS_BIG_INT_TAG
    end

    encode_unsigned_with_type(TYPE_6, Unsigned(tag), bytes)

    if isodd(length(hex_str))
        hex_str = "0" * hex_str
    end

    encode(hex2bytes(hex_str), bytes)
end

function encode(float::Union{Float64, Float32, Float16}, bytes::Array{UInt8, 1})
    if typeof(float) == Float64 && float == Float32(float)
        float_bytes = hex2bytes(num2hex(Float32(float)) )
    else
        float_bytes = hex2bytes(num2hex(float))
    end
    float_bytes_len = length(float_bytes)

    if float_bytes_len == SIZE_OF_FLOAT64
        push!(bytes, TYPE_7 | ADDNTL_INFO_FLOAT64)
    elseif float_bytes_len == SIZE_OF_FLOAT32
        push!(bytes, TYPE_7 | ADDNTL_INFO_FLOAT32)
    else float_bytes_len == SIZE_OF_FLOAT16
        push!(bytes, TYPE_7 | ADDNTL_INFO_FLOAT16)
    end

    append!(bytes, float_bytes)
end

if VERSION < v"0.6.0"
    function encode(pair::Pair, bytes::Array{UInt8, 1})
        if typeof(pair.first) <: Integer && pair.first >= 0
            encode_with_tag(Unsigned(pair.first), pair.second, bytes)
        elseif typeof(pair.first) == Task
            encode_indef_length_collection(pair.first, pair.second, bytes)
        else
            encode_custom_type(pair, bytes)
        end
    end
else
    function encode(pair::Pair, bytes::Array{UInt8, 1})
        if typeof(pair.first) <: Integer && pair.first >= 0
            encode_with_tag(Unsigned(pair.first), pair.second, bytes)
        elseif typeof(pair.first) <: Channel
            encode_indef_length_collection(pair.first, pair.second, bytes)
        else
            encode_custom_type(pair, bytes)
        end
    end
end

function encode(data, bytes::Array{UInt8, 1})
    encode_custom_type(data, bytes)
end

# ------- encoding for indefinite length collections

function encode_indef_length_collection(producer::Channel, collection_type,
                                        bytes::Array{UInt8, 1})
    if collection_type <: AbstractVector{UInt8}
        typ = TYPE_2
    elseif collection_type <: String
        typ = TYPE_3
    elseif collection_type <: Union{AbstractVector, Tuple}
        typ = TYPE_4
    elseif collection_type <: AbstractDict
        typ = TYPE_5
    else
        error(@sprintf "Collection type %s is not supported for indefinite length encoding." collection_type)
    end

    push!(bytes, typ | ADDNTL_INFO_INDEF)

    count = 0
    for e in producer
        encode(e, bytes)
        count += 1
    end

    if typ == TYPE_5 && isodd(count)
        error(@sprintf "Collection type %s requires an even number of input data items in order to be consistent." collection_type)
    end

    push!(bytes, BREAK_INDEF)
end

# ------- encoding with tags

function encode_with_tag(tag::Unsigned, data, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_6, tag, bytes)
    encode(data, bytes)
end

# ------- encoding for user-defined types

if VERSION < v"0.5.0"
    function encode_custom_type(data, bytes::Array{UInt8, 1})
        type_map = Dict()

        type_map[UTF8String("type")] = UTF8String(string(typeof(data)) )

        for f in fieldnames(data)
            type_map[UTF8String(string(f))] = data.(f)
        end

        encode(type_map, bytes)
    end
else
    function encode_custom_type(data, bytes::Array{UInt8, 1})
        type_map = Dict()

        type_map[String("type")] = String(string(typeof(data)) )

        for f in fieldnames(data)
            type_map[String(string(f))] = data.(f)
        end

        encode(type_map, bytes)
    end
end

# ------- encoding for Simple types

struct Simple
    val::UInt8
end

Base.isequal(a::Simple, b::Simple) = Base.isequal(a.val, b.val)

function encode(simple::Simple, bytes::Array{UInt8, 1})
    encode_unsigned_with_type(TYPE_7, simple.val, bytes)
end

# ------- encoding for Null and Undefined

struct Null
end

struct Undefined
end

function encode(null::Null, bytes::Array{UInt8, 1})
    push!(bytes, CBOR_NULL_BYTE)
end

function encode(undef::Undefined, bytes::Array{UInt8, 1})
    push!(bytes, CBOR_UNDEF_BYTE)
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using BasicPOMCP
using POMCPOW
using ParticleFilters
using CPUTime
using VDPTag
using DataFrames
using ARDESPOT

N = 100

pomdp = VDPTagPOMDP()
@show mdp(pomdp).tag_radius
@show mdp(pomdp).tag_terminate

file_contents = readstring(@__FILE__())

solvers = Dict{String, Union{Solver,Policy}}(

    "pomcpow" => begin
        rng = MersenneTwister(13)
        ro = ToNextMLSolver(rng)
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(40.0),
                               final_criterion=MaxTries(),
                               max_depth=10,
                               max_time=0.1,
                               k_action=8.0,
                               alpha_action=1/20,
                               k_observation=4.0,
                               alpha_observation=1/20,
                               estimate_value=FORollout(ro),
                               check_repeat_act=true,
                               check_repeat_obs=true,
                               next_action=RootToNextMLFirst(rng),
                               default_action=ReportWhenUsed(1),
                               rng=rng
                              )
    end,

    "pomcp" => begin
        rng = MersenneTwister(13)
        ro = ToNextMLSolver(rng)
        POMCPSolver(max_depth=10,
                    max_time=0.1,
                    c=40.0,
                    tree_queries=typemax(Int),
                    default_action=ReportWhenUsed(1),
                    estimate_value=FORollout(ro),
                    rng=rng
                   )
    end,

    "despot" => begin
        rng = MersenneTwister(13)
        ro = ToNextMLSolver(rng)
        b = IndependentBounds(DefaultPolicyLB(ro), VDPUpper())
        DESPOTSolver(lambda=0.01,
                     K=500,
                     D=10,
                     max_trials=1_000_000,
                     T_max=0.1,
                     bounds=b,
                     random_source=MemorizingSource(500, 10, rng, min_reserve=8),
                     rng=rng)
    end
)

alldata = DataFrame()
for n_angles_float in logspace(0.5, 3, 6)
# for n_obs_angles_float in logspace(0.5, 3, 6)
# for n_angles_float in [10.0]
    n_angles = round(Int, n_angles_float)
    for (k, solver) in solvers
        println("$k ($n_angles)")
        dpomdp = ADiscreteVDPTagPOMDP(n_angles=n_angles)
        planner = solve(solver, dpomdp)
        sims = []
        for i in 1:N
            srand(planner, i+50_000)
            cplanner = translate_policy(planner, dpomdp, pomdp, dpomdp)

            filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                               LowVarianceResampler(10_000),
                                               0.05, MersenneTwister(i+90_000))

            sim = Sim(deepcopy(pomdp),
                      deepcopy(cplanner),
                      filter,
                      rng=MersenneTwister(i+70_000),
                      max_steps=100,
                      metadata=Dict(:solver=>k,
                                    :n_angles=>n_angles,
                                    :i=>i)
                     )

            push!(sims, sim)
        end

        data = run_parallel(sims)
        # data = run(sims)

        rs = data[:reward]
        println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
        alldata = vcat(alldata, data)
    end
end
# end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "vdp_a_discretization_$(Dates.format(now(), "E_d_u_HH_MM"))")
println("saving to $filename...")
writetable(filename*".csv", alldata)
write(filename*"_file.jl", file_contents)
println("done.")


text:
"""
Created on Thu Apr 23 2020

@author: Yoann Pradat

Cost function for nmf_α
"""

"""
    _α_cost(V, W, H, params)

Compute the cost function that is minimized by the NMF alogrithm running with parameters params. It adds up
the α_divergence and regularization terms.

See also: [`α_divergence`](@ref), [`CParams`](@ref)
"""
function _α_cost(V::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::CParams) where T <: Real
    regu_H = params.α_H * (params.l₁ratio_H * norm(H, 1) + (1. - params.l₁ratio_H) * norm(H, 2)) / length(H)
    regu_W = params.α_W * (params.l₁ratio_W * norm(W, 1) + (1. - params.l₁ratio_W) * norm(W, 2)) / length(W)
    α_div  = α_divergence(V, W*H, params.α)

    α_div + regu_H + regu_W
end


text:
# SPDX-License-Identifier: X11
# 2020-12-07
# Day 3, Part 2

const paths = [
  (1,1),
  (3,1),
  (5,1),
  (7,1),
  (1,2)
]

function main()
  input = readlines()
  # [cols][rows]

  (println ∘ prod ∘ map)(x -> counttrees(input, x),
                         paths::Vector{Tuple{Int,Int}})
end

function counttrees(input::Vector{String}, tup::Tuple{Int,Int})::Int
  right = tup[1]
  down = tup[2]

  cnt = 0
  x = 1 + right
  bound = length(input[1])

  for i ∈ 1 + down:down:length(input)
    input[i][x] == '#' && (cnt += 1)
    x += right
    x > bound && (x = x - bound)
  end

  return cnt
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Test
using ChemistryFeaturization.AbstractType:
    AbstractAtoms, AbstractFeaturization, AbstractFeatureDescriptor, AbstractCodec
using ChemistryFeaturization.Atoms: elements
using ChemistryFeaturization: elements

@testset "Modules and Abstract methods" begin
    struct FakeAtoms <: AbstractAtoms{Nothing} end
    struct FakeFD <: AbstractFeatureDescriptor end
    struct FakeFeaturization <: AbstractFeaturization end
    struct FakeCodec <: AbstractCodec end

    @testset "top-level module" begin
        # testing on `nothing` as example of ::Any
        @test_throws MethodError encodable_elements(nothing)
        @test_throws MethodError encode(nothing, nothing)
        @test_throws MethodError decode(nothing, nothing)
        @test_throws MethodError elements(nothing)
    end

    @testset "atoms module" begin
        fa = FakeAtoms()
        @test_throws MethodError elements(fa)
    end

    @testset "features module" begin
        fa = FakeAtoms()
        ffd = FakeFD()
        fc = FakeCodec()
        @test_throws MethodError encodable_elements(ffd)
        @test_throws MethodError get_value(ffd, fa)
        @test_throws ErrorException encode(ffd, fa)
        @test_throws ErrorException decode(ffd, nothing)
    end

    @testset "featurizations module" begin
        ff = FakeFeaturization()
        F2_atom = AtomGraph(Float32.([0 1; 1 0]), ["F", "F"])

        @test_throws MethodError encodable_elements(ff)
        @test_throws MethodError encode(ff, F2_atom)
        @test_throws MethodError decode(ff, nothing)
    end
end


text:
using Quantum, Test

@testset "Quantum" begin
    @testset "QuantumRegister" begin
        @testset "Construction" begin
            @test QuantumRegister(1) == [1, 0]
            @test QuantumRegister(2) == [1, 0, 0, 0]
            @test QuantumRegister(3) == [1, 0, 0, 0, 0, 0, 0, 0]
            @test QuantumRegister(2) == QuantumRegister(2)
            @test [1, 0, 0, 0] == QuantumRegister(2)
        end
        @testset "Length" begin
            @test length(QuantumRegister(1)) == 1
            @test length(QuantumRegister(2)) == 2
            @test length(QuantumRegister(3)) == 3
            @test length(QuantumRegister(4)) == 4
        end
    end
    @testset "Basic Gates" begin
        @testset "X" begin
            @test X(QuantumRegister(1), 1) == [0, 1]
            @test X(QuantumRegister(2), 1) == [0, 0, 1, 0]
            @test X(QuantumRegister(2), 2) == [0, 1, 0, 0]
            @test X(QuantumRegister(3), 2) == [0, 0, 1, 0, 0, 0, 0, 0]
        end
        @testset "H" begin
            @test H(QuantumRegister(1), 1) == [1 / √2, 1 / √2]
            @test H(X(QuantumRegister(1), 1), 1) == [1 / √2, -1 / √2]
        end
        @testset "CNOT" begin
            @test CNOT(QuantumRegister(2), 1, 2) == [1, 0, 0, 0]
            @test CNOT(X(QuantumRegister(2), 1), 1, 2) == [0, 0, 0, 1]
            @test CNOT(X(QuantumRegister(2), 2), 1, 2) == [0, 1, 0, 0]
            @test CNOT(X(QuantumRegister(2), 1), 2, 1) == [0, 0, 1, 0]
            @test CNOT(X(QuantumRegister(2), 2), 2, 1) == [0, 0, 0, 1]
            @test CNOT(X(QuantumRegister(3), 1), 1, 3) == [0, 0, 0, 0, 0, 1, 0, 0]
        end
        @testset "Swap" begin
            @test Swap(QuantumRegister(2), 1, 2) == [1, 0, 0, 0]
            @test Swap(QuantumRegister(2), 2, 1) == [1, 0, 0, 0]
            @test Swap(X(QuantumRegister(2), 1), 1, 2) == [0, 1, 0, 0]
            @test Swap(X(QuantumRegister(2), 2), 1, 2) == [0, 0, 1, 0]
            @test Swap(X(QuantumRegister(2), 1), 2, 1) == [0, 1, 0, 0]
        end
        @testset "Entanglement" begin
            @test CNOT(H(QuantumRegister(2), 1), 1, 2) == [1 / √2, 0, 0, 1 / √2]
            @test CNOT(H(X(QuantumRegister(2), 2), 1), 1, 2) == [0, 1 / √2, 1 / √2, 0]
        end
    end
    @testset "Single Qubit" begin
        @testset "Pauli Gates" begin
            @test Y(QuantumRegister(1), 1) == [0.0, im]
            @test Y(X(QuantumRegister(1), 1), 1) == [-im, 0.0]
            @test Z(QuantumRegister(1), 1) == [1, 0]
            @test Z(X(QuantumRegister(1), 1), 1) == [0, -1]
        end
        @testset "S(dag)" begin
            @test S(QuantumRegister(1), 1) == [1, 0]
            @test S(X(QuantumRegister(1), 1), 1) == [0, im]
            @test S(QuantumRegister(1), 1) == T(T(QuantumRegister(1), 1), 1)
            @test Sdag(QuantumRegister(1), 1) == [1, 0]
            @test Sdag(X(QuantumRegister(1), 1), 1) == [0, -im]
        end
        @testset "T(dag)" begin
            @test T(QuantumRegister(1), 1) == [1, 0]
            @test T(X(QuantumRegister(1), 1), 1) == [0, exp(im * pi / 4)]
            @test T(T(QuantumRegister(1), 1), 1) == S(QuantumRegister(1), 1)
            @test Tdag(QuantumRegister(1), 1) == [1, 0]
            @test Tdag(X(QuantumRegister(1), 1), 1) == [0, exp(-im * pi / 4)]
        end
    end
    @testset "Controlled Gates" begin
        @testset "Single Control" begin
            @test CZ(QuantumRegister(2), 1, 2) == [1, 0, 0, 0]
            @test CZ(X(QuantumRegister(2), 1), 1, 2) == [0, 0, 1, 0]
            @test CZ(X(QuantumRegister(2), 2), 1, 2) == [0, 1, 0, 0]
            @test CZ(X(QuantumRegister(2), 1), 2, 1) == [0, 0, 1, 0]
            @test CZ(X(QuantumRegister(2), 2), 2, 1) == [0, 1, 0, 0]
            @test CZ(X(QuantumRegister(3), 1), 1, 3) == [0, 0, 0, 0, 1, 0, 0, 0]
            @test CNOT(H(QuantumRegister(2), 2), 2, 1) == [1 / √2, 0, 0, 1 / √2]
            @test CNOT(H(X(QuantumRegister(2), 1), 2), 2, 1) == [0, 1 / √2, 1 / √2, 0]
        end
    end
end


text:
module OSMMakie

using LightOSM
using GraphMakie
using Graphs
using Makie
using Makie.Colors
using Makie.GeometryBasics

include("recipe.jl")
include("defaults/defaults.jl")
include("inspection.jl")

end # module


text:
_index_extrap(etp, inds) = etp(inds...)
function warp_and_resample(img, tfm; fillval=0.0)
    inds0 = axes(img)
    img = warp(img, tfm)
    itp = interpolate(img, BSpline(Interpolations.Linear()))
    etp = extrapolate(itp, fillval)
    return _index_extrap(etp, inds0) #Function barrier improves performance
end


text:
using GLAbstraction, GeometryTypes, ModernGL, Compat, FileIO, GLWindow, FixedSizeArrays, FixedPointNumbers, ColorTypes
using Base.Test

function is_ci()
	get(ENV, "TRAVIS", "") == "true" || get(ENV, "APPVEYOR", "") == "true" || get(ENV, "CI", "") == "true"
end

include("macro_test.jl")

if !is_ci() # only do test if not CI... this is for automated testing environments which fail for OpenGL stuff, but I'd like to test if at least including works

window = create_glcontext("test", resolution=(500,500))

include("accessors.jl")
include("uniforms.jl")
include("texture.jl")

# Test for creating a GLBuffer with a 1D Julia Array
# You need to supply the cardinality, as it can't be inferred
# indexbuffer is a shortcut for GLBuffer(GLUInt[0,1,2,2,3,0], 1, buffertype = GL_ELEMENT_ARRAY_BUFFER)
indexes = indexbuffer(GLuint[0,1,2])
# Test for creating a GLBuffer with a 1D Julia Array of Vectors
#v = Vec2f[Vec2f(0.0, 0.5), Vec2f(0.5, -0.5), Vec2f(-0.5,-0.5)]

v = Vec2f0[Vec2f0(0.0, 0.5), Vec2f0(0.5, -0.5), Vec2f0(-0.5,-0.5)]

verts = GLBuffer(v)
@test size(verts, 1) == 3
@test size(verts, 2) == 1
# lets define some uniforms
# uniforms are shader variables, which are supposed to stay the same for an entire draw call

const triangle = RenderObject(
	Dict(
		:vertex => verts,
		:name_doesnt_matter_for_indexes => indexes
	),
	LazyShader(load("test.vert"), load("test.frag"))
)

postrender!(triangle, render, triangle.vertexarray)

glClearColor(0,0,0,1)
while isopen(window)
  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
	render(triangle)
	swapbuffers(window)
	pollevents()
	sleep(0.01)
end

end


text:
function copy_self(obj::Hyperplane)
    return jcall(obj, "copySelf", Hyperplane, ())
end

function empty_hyperplane(obj::Hyperplane)
    return jcall(obj, "emptyHyperplane", SubHyperplane, ())
end

function get_offset(obj::Hyperplane, arg0::Point)
    return jcall(obj, "getOffset", jdouble, (Point,), arg0)
end

function get_tolerance(obj::Hyperplane)
    return jcall(obj, "getTolerance", jdouble, ())
end

function project(obj::Hyperplane, arg0::Point)
    return jcall(obj, "project", Point, (Point,), arg0)
end

function same_orientation_as(obj::Hyperplane, arg0::Hyperplane)
    return jcall(obj, "sameOrientationAs", jboolean, (Hyperplane,), arg0)
end

function whole_hyperplane(obj::Hyperplane)
    return jcall(obj, "wholeHyperplane", SubHyperplane, ())
end

function whole_space(obj::Hyperplane)
    return jcall(obj, "wholeSpace", Region, ())
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    x = y = -5:1:5
    z = [x^2 + y^2 for x in x, y in y]
    fig = Figure(resolution = (450, 400))
    ax = Axis(fig,  aspect = 1, xlabel = "x", ylabel = "y")
    hmap = heatmap!(x,y,z, colormap = :ice)
    cbar = Colorbar(fig, hmap, label = "z values", width = 15, ticksize=15,
                tickalign = 1, height = Relative(0.945))
    ax.xticks = -5:5
    ax.yticks = -5:5
    fig[1, 1] = ax
    fig[1, 2] = cbar
    colgap!(fig.layout, 7)
    fig
    save(joinpath(@__DIR__, "output", "heatmapxyz.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie