text:
#################
# DictionaryStyle
#################

export DictionaryStyle, UnknownDictionaryStyle, BasisStyle, FrameStyle
"""
    abstract type DictionaryStyle end

A dictionary style trait for platforms

See also, [`BasisStyle`](@ref), [`FrameStyle`](@ref),
[`UnknownDictionaryStyle`](@ref)
"""
abstract type DictionaryStyle end

"""
    struct UnknownDictionaryStyle <: DictionaryStyle end

The platform creates dictionaries of unknown format.
"""
struct UnknownDictionaryStyle <: DictionaryStyle end

"""
    struct BasisStyle <: DictionaryStyle end

The platform creates dictionaries that are bases.
"""
struct BasisStyle <: DictionaryStyle end

"""
    struct FrameStyle <: DictionaryStyle end

The platform creates dictionaries that are frames.
"""
struct FrameStyle <: DictionaryStyle end

DictionaryStyle(platform::Platform) = UnknownDictionaryStyle()
DictionaryStyle(platform::BasisPlatform) = BasisStyle()
DictionaryStyle(platform::FramePlatform) = FrameStyle()

#################
# SamplingStyle
#################

export SamplingStyle, DiscreteStyle, ProjectionStyle, InterpolationStyle,
    OversamplingStyle, GridStyle, GramStyle, DiscreteGramStyle,
    RectangularGramStyle, GenericSamplingStyle, ProductSamplingStyle
"""
    abstract type SamplingStyle end

A sampling style trait for platforms

See also [`DiscreteStyle`](@ref), [`ProjectionStyle`](@ref),
[`GenericSamplingStyle`](@ref), [`ProductSamplingStyle`](@ref)
"""
abstract type SamplingStyle end

"""
    abstract type DiscreteStyle <: SamplingStyle end

The sampling operator corresponds to evaluation in a grid. This grid is
specified by the listed styles here:

[`InterpolationStyle`](@ref), [`OversamplingStyle`](@ref),
[`GridStyle`](@ref)
"""
abstract type DiscreteStyle <: SamplingStyle end

"Discrete sampling without associated weights."
abstract type UnweightedDiscreteStyle <: DiscreteStyle end
"Discrete sampling with associated weights."
abstract type WeightedDiscreteStyle <: DiscreteStyle end

"The sampling grid is determined by invoking `interpolation_grid` on the dictionary or platform."
struct InterpolationStyle <: UnweightedDiscreteStyle end
"The sampling grid is determined by invoking `oversampling_grid` on the dictionary or platform."
struct OversamplingStyle <: UnweightedDiscreteStyle end
"""
The sampling grid is determined by invoking `platform_grid` on the platform,
or by providing a `grid=...` argument to the Fun constructor.
"""
struct GridStyle <: UnweightedDiscreteStyle end

"Discrete sampling followed by weighting the samples."
struct WeightedSampling <: WeightedDiscreteStyle
    style       # the discrete sampling style
end
unweighted(st::WeightedSampling) = st.style

"Discrete sampling followed by normalizing the samples."
struct NormalizedSampling <: WeightedDiscreteStyle
    style       # the discrete sampling style
end
unweighted(st::NormalizedSampling) = st.style

NormalizedSampling(style::NormalizedSampling) = style

"""
    abstract type ProjectionStyle <: SamplingStyle end

The sampling operator corresponds to projections (using inner products).
The kind of inner product is specified by
[`GramStyle`](@ref), [`DiscreteGramStyle`](@ref),
[`RectangularGramStyle`](@ref)
"""
abstract type ProjectionStyle <: SamplingStyle end
"""
    struct GramStyle <: ProjectionStyle end

Use the inner product implied by the measure obtained by invoking `measure` on the
dictionary or platform. This results in a Gram matrix.
"""
struct GramStyle <: ProjectionStyle end
"""
    struct DiscreteGramStyle <: ProjectionStyle end

Use the inner product implied by the measure obtained by invoking `discretemeasure` on the
dictionary or platform. This results in a discrete Gram matrix.
"""
struct DiscreteGramStyle <: ProjectionStyle end
"""
    struct RectangularGramStyle <: ProjectionStyle end

Not implemented (oversample with a projection on a larger dual dictionary)
"""
struct RectangularGramStyle <: ProjectionStyle end

"""
    struct GenericSamplingStyle <: SamplingStyle end

Write your own sampling_operator. The sampling operator corresponds the one implemented
by the user.
"""
struct GenericSamplingStyle <: SamplingStyle end

"""
    struct ProductSamplingStyle{STYLES} <: SamplingStyle

The sampling operator has product structure.
"""
struct ProductSamplingStyle{STYLES} <: SamplingStyle
    styles  :: STYLES
end
ProductSamplingStyle(styles::SamplingStyle...) = ProductSamplingStyle(styles)
components(style::ProductSamplingStyle) = style.styles
factors(style::ProductSamplingStyle) = components(style)

# TODO: remove this definition
export SamplingStyleSuper
const SamplingStyleSuper{TYPE} = Union{<:TYPE,ProductSamplingStyle{NTuple{N,<:TYPE}} where N} where {TYPE<:SamplingStyle}

#################
# SolverStyle
#################

export SolverStyle, DirectStyle, InverseStyle, DualStyle, IterativeStyle,
    TridiagonalProlateStyle, AZStyle, AZSmoothStyle, ProductSolverStyle
"""
    abstract type SolverStyle end

A solver style trait for platforms

See also [`DirectStyle`](@ref), [`InverseStyle`](@ref),
[`DualStyle`](@ref), [`IterativeStyle`](@ref),
[`TridiagonalProlateStyle`](@ref), [`AZStyle`](@ref),
[`AZSmoothStyle`](@ref), [`ProductSolverStyle`](@ref)
"""
abstract type SolverStyle end
"""
    struct DirectStyle <: SolverStyle end

Use a direct solver to solve the discretized approximation problem.
"""
struct DirectStyle <: SolverStyle end
"""
    struct InverseStyle <: SolverStyle end

Use the inverse operator to solve the discretized approximation problem.
"""
struct InverseStyle <: SolverStyle end
"""
    struct DualStyle <: SolverStyle end

Use projection on the dual dictionary to solve the discretized approximation problem.
"""
struct DualStyle <: SolverStyle end
"""
    struct IterativeStyle <: SolverStyle end

Use an iterative method to solve the discretized approximation problem.
"""
struct IterativeStyle <: SolverStyle end
"""
    struct TridiagonalProlateStyle <: SolverStyle end

Use a tridiagonal solver to solve the discretized approximation problem
(only for Fourier extension)
"""
struct TridiagonalProlateStyle <: SolverStyle end
"""
    struct AZStyle <: SolverStyle end

Use the AZ algorithm to solve the discretized approximation problem.
"""
struct AZStyle <: SolverStyle end
"""
    struct AZSmoothStyle <: SolverStyle end

Use the AZ algorithm to solve the discretized approximation problem,
but tweak the coefficients such that a smooth extension is obtained.
"""
struct AZSmoothStyle <: SolverStyle end

"""
    struct ProductSolverStyle <: SolverStyle

Use a tensor product solver to solve the discretized approximation problem.
"""
struct ProductSolverStyle <: SolverStyle
    styles
end
ProductSolverStyle(styles::SolverStyle...) = ProductSolverStyle(styles)
components(style::ProductSolverStyle) = style.styles
factors(style::ProductSolverStyle) = components(style)

#################
# Defaults
#################

# Defaults for platforms
SamplingStyle(p::Platform) = SamplingStyle(DictionaryStyle(p), p)

SamplingStyle(::BasisStyle, p::Platform) = InterpolationStyle()
SamplingStyle(::FrameStyle, p::Platform) = OversamplingStyle()
SamplingStyle(::UnknownDictionaryStyle, p::Platform) = OversamplingStyle()

SolverStyle(p::Platform, samplingstyle) =
    SolverStyle(DictionaryStyle(p), p, samplingstyle)

SolverStyle(::BasisStyle, p::Platform, ::InterpolationStyle) = InverseStyle()
SolverStyle(::BasisStyle, p::Platform, ::OversamplingStyle) = DirectStyle()
SolverStyle(::BasisStyle, p::Platform, ::GridStyle) = DirectStyle()
SolverStyle(::BasisStyle, p::Platform, ::SamplingStyle) = DirectStyle()
SolverStyle(::FrameStyle, p::Platform, ::SamplingStyle) = AZStyle()
SolverStyle(::UnknownDictionaryStyle, p::Platform, ::SamplingStyle) = DirectStyle()


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

"""
    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)

The classical inviscid Taylor-Green vortex.
"""
function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
  A  = 1.0 # magnitude of speed
  Ms = 0.1 # maximum Mach number

  rho = 1.0
  v1  =  A * sin(x[1]) * cos(x[2]) * cos(x[3])
  v2  = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
  v3  = 0.0
  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
  p   = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3]))

  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
end
initial_condition = initial_condition_taylor_green_vortex

volume_flux  = flux_ranocha
surface_flux = flux_lax_friedrichs
solver = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(),
                surface_integral= SurfaceIntegralWeakForm(surface_flux),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

mesh = DGMultiMesh(solver, cells_per_dimension=(8, 8, 8),
                   coordinates_min=(-pi, -pi, -pi), coordinates_max=(pi, pi, pi),
                   periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 25.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7,
            save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary

ext: jl


text:
#
# test LowRankMatrix
#
using ApproxFun, SingularIntegralEquations, Base.Test

L = lrzeros(Complex{Float64},10,20)
@test size(L) == (10,20)

@test full(L) == zeros(Complex{Float64},10,20)

@test_approx_eq full(lrones(10,10)) ones(10,10)

@test rank(lreye(10,10)) == 10

L = LowRankMatrix(1./collect(1:10),1./collect(1:10.0))

x = collect(linspace(-1,1,10))

@test norm(L*x-full(L)*x) < 200eps()

@test isa(L*L',LowRankMatrix)

A = rand(10,10)

@test isa(A*L,LowRankMatrix)
@test isa(L*A,LowRankMatrix)

@test isa(A+L,Matrix)
@test isa(L⊕A,LowRankMatrix)

@test rank(L+L) == 2rank(L⊕L)
# @test rank(L-L) == 2rank(L⊖L)
@test rank(2L+1+full(L)) ≥ rank(L⊕full(L))


text:
@testset "limit HMTL output" begin

@static VERSION ≥ v"1.4" && let

# no limit
doc = mock_run("""
```julia
using DataFrames
DataFrame(rand(10,3))
```
"""; doctype = "md2html")
@test isdefined(doc.chunks[1], :rich_output)
@test count("<tr>", doc.chunks[1].rich_output) == 12 # additonal 2 for name and type row

# limit
n = 100000
doc = mock_run("""
```julia
using DataFrames
DataFrame(rand($n,3))
```
"""; doctype = "md2html")
@test isdefined(doc.chunks[1], :rich_output)
@test count("<tr>", doc.chunks[1].rich_output) < n

end

end


text:
module Quoridor

import Base: (+), (-), isless, mean,convert, show, getindex, rand, start, next, done

import Graphs: AbstractGraph, Edge
using ImmutableArrays

const BOARD_MAX = 10
const BOARD_MIN = 1

@enum Player Player1=1 Player2=-1

function other(p::Player)
    return Player(-Int(p))
end

rand(::Type{Player}) = rand(Bool) ? Player1 : Player2

@enum Direction Horizontal=1 Vertical=-1

include("loc.jl")
include("wall.jl")
include("move.jl")
include("playerstate.jl")
include("game.jl")
include("graphimpl.jl")

end # module


text:
include("../TestFunctions.jl")
include("indexing.jl")
include("labelledmatrix.jl")
include("windows.jl")


text:
include("dependencies_for_runtests.jl")

using Statistics: mean
using Oceananigans.CubedSpheres
using Oceananigans.Models.HydrostaticFreeSurfaceModels
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField

@testset "Cubed spheres" begin

    @testset "Conformal cubed sphere grid" begin
        @info "  Testing conformal cubed sphere grid..."

        grid = ConformalCubedSphereGrid(face_size=(10, 10, 1), z=(-1, 0))
        @test try show(grid); println(); true; catch; false; end
    end

    for arch in archs

        @info "  Constructing a ConformalCubedSphereGrid from file [$(typeof(arch))]..."

        # Prototype grid and model for subsequent tests
        cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"
        grid = ConformalCubedSphereGrid(cs32_filepath, arch, Nz=1, z=(-1, 0))

        @info "  Constructing a HydrostaticFreeSurfaceModel on a ConformalCubedSphereGrid [$(typeof(arch))]..."

        free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1)
        model = HydrostaticFreeSurfaceModel(; grid, free_surface,
                                            momentum_advection = VectorInvariant(),
                                            coriolis = nothing,
                                            closure = nothing,
                                            tracers = :c,
                                            buoyancy = nothing)

        @testset "Constructing a grid from file [$(typeof(arch))]" begin
            @test grid isa ConformalCubedSphereGrid
        end

        @testset "CubedSphereData and CubedSphereFields [$(typeof(arch))]" begin
            @info "  Testing CubedSphereData and CubedSphereFields [$(typeof(arch))]..."
            c = model.tracers.c
            η = model.free_surface.η

            set!(c, 0)
            set!(η, 0)

            CUDA.allowscalar(true)
            @test all(all(face_c .== 0) for face_c in faces(c))
            @test all(all(face_η .== 0) for face_η in faces(η))
            CUDA.allowscalar(false)

            @test maximum(abs, c) == 0
            @test minimum(abs, c) == 0
            @test mean(c) == 0

            @test maximum(abs, η) == 0
            @test minimum(abs, η) == 0
            @test mean(η) == 0
        end

        @testset "Constructing a HydrostaticFreeSurfaceModel [$(typeof(arch))]" begin
            @test model isa HydrostaticFreeSurfaceModel
        end

        @testset "Time stepping a HydrostaticFreeSurfaceModel [$(typeof(arch))]" begin
            @info "  Time-stepping HydrostaticFreeSurfaceModel on a ConformalCubedSphereGrid [$(typeof(arch))]..."
            time_step!(model, 1)
            @test try time_step!(model, 1); true; catch; false; end
        end

        @testset "VerticalVorticityField on ConformalCubedSphereGrid [$(typeof(arch))]" begin
            @info "  Testing VerticalVorticityField on a ConformalCubedSphereGrid [$(typeof(arch))]..."
            ζ = VerticalVorticityField(model)

            @test ζ isa Field

            set!(model, u = (x, y, z) -> rand())

            @test try
                compute!(ζ)
                true
            catch err
                println(sprint(showerror, err))
                false
            end
            @test maximum(abs, ζ) > 0 # fingers crossed
        end
    end
end


text:
#!/usr/bin/env julia
# Simple REPL to test the parser and compare it against the base
# julia parser - ^R reloads the package, ^S switches between base
# and JuliaParser parsers. Entered lines are parsed and printed
# in tree forms.

using AbstractTrees
using JuliaParser
using TerminalUI
using VT100
using Colors

dosl = false
if length(ARGS) == 1 && ARGS[1] == "--sl"
    dosl = true
end

eval(Base,:(have_color = true))

import AbstractTrees: children, printnode

function printnode(io::IO, x::Expr, color = :default)
    print_with_color(:red, io, string(':', x.head))
end

function printnode(io::IO, x::LineNumberNode)
    print_with_color(:blue, io, string(x))
end

function printnode(io::IO, x::QuoteNode)
    print_with_color(:green, io, string(x))
end

function printnode(io::IO, x::GlobalRef)
    print_with_color(:yellow, io, string(x))
end

import Base: LineEdit, REPL

function reload_widget()
    eval(Main,:(module LocWidget
        include($(joinpath(dirname(@__FILE__),"..","src","interactiveutil.jl")))
    end))
end
reload_widget()

function RunShell()
    global hp
    prompt = "parse > "

    isBase = false

    # Setup debug panel
    panel = LineEdit.Prompt(prompt;
        prompt_prefix=()->isBase?"\e[38;5;107m":"\e[38;5;166m",
        prompt_suffix=Base.text_colors[:white],
        on_enter = dosl ? (s)->true : Base.REPL.return_callback)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:debug => panel))

    panel.on_done = (s,buf,ok)->begin
        line = String(take!(buf))
        if !ok || strip(line) == "q"
            LineEdit.transition(s, :abort)
        end
        try
            if isBase
                show(STDOUT, Tree(Base.parse(line)))
            elseif dosl
                eval(quote
                    ts = Main.JuliaParser.Lexer.TokenStream{Main.JuliaParser.Lexer.SourceLocToken}($line)
                    try
                        res = Main.JuliaParser.Parser.parse(ts)
                        if res != nothing
                            show(STDOUT,Tree(res.expr))
                            println($line)
                            w = Main.LocWidget.create_widget(res.loc,$(line))
                            TerminalUI.print_snapshot(TerminalUI.InlineDialog(w,
                                Base.Terminals.TTYTerminal("xterm", STDIN, STDOUT, STDERR)
                                ))
                        end
                    catch e
                        !isa(e, Main.JuliaParser.Diagnostics.Diagnostic) && rethrow(e)
                        Main.JuliaParser.Diagnostics.display_diagnostic(STDOUT, $line, e)
                    end
                end)
            else
                show(STDOUT,Tree(eval(:(Main.JuliaParser.Parser.parse($line)))))
            end
        catch err
            REPL.display_error(STDERR, err, Base.catch_backtrace())
            REPL.println(STDERR); REPL.println(STDERR)
        end
        LineEdit.reset_state(s)
        println()
    end

    hp = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:parse => panel))
    panel.hist = hp

    REPL.history_reset_state(hp)

    extra_keymap = Dict{Any,Any}(
      "^R" => (s,o...)->(print("Reloading... "); reload("JuliaParser");
        reload_widget();
        println("Done!"); LineEdit.refresh_line(s)),
      "^S" => (s,o...)->(isBase = !isBase; LineEdit.refresh_line(s))
    );
    search_prompt, skeymap = LineEdit.setup_prefix_keymap(hp, panel)

    b = Dict{Any,Any}[extra_keymap, skeymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap(b)

    t = Base.Terminals.TTYTerminal("xterm", STDIN, STDOUT, STDERR)
    Base.REPL.run_interface(t, LineEdit.ModalInterface([panel,search_prompt]))
end
RunShell()


text: # example config file


text:
using JuMP, Clp, Ipopt

ROOT=pwd()
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

"""
    get_JuMP_cartesian_model(problem_poly::Problem, mysolver)

Return JuMP cartesian model associated to `mysolver` defined by `problem_poly`, a polynomial problem in complex numbers

# Example
```jldoctest
V1 = Variable("VOLT_1",Complex)
V2 = Variable("VOLT_2",Complex)
p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
p_ctr1 = abs2(V1)
p_ctr2 = 3im * V1*conj(V2) + abs2(V1)

problem_poly=Problem()
add_variable!(problem_poly,V1)
add_variable!(problem_poly,V2)
add_constraint!(problem_poly, "ctr1", 0.95^2 << p_ctr1 << 1.05^2)
add_constraint!(problem_poly, "ctr2", p_ctr2==0)
set_objective!(problem_poly, p_obj)
print(problem_poly)
mysolver = ClpSolver()
m = get_JuMP_cartesian_model(problem_poly, mysolver)
print(m)
```
"""
function get_JuMP_cartesian_model(problem_poly::Problem, mysolver)
    pb_poly_real = pb_cplx2real(problem_poly)
    m = Model(solver = mysolver)
    variables_jump = Dict{String, JuMP.Variable}()
    for (varname, vartype) in pb_poly_real.variables
        if vartype<:Real
            variables_jump["$varname"] = @variable(m, varname, basename="$varname")
        end
    end
    ctr_jump = Dict{String,JuMP.ConstraintRef}()
    for (ctr, modeler_ctr) in pb_poly_real.constraints
        polynome = modeler_ctr.p
        lb = modeler_ctr.lb
        ub = modeler_ctr.ub

        for value in values(polynome)
            if imag(value)!=0
                error("Polynom coefficients have to be real numbers")
            end
        end
        ctr_jump[ctr] = @NLconstraint(m, lb <= sum(coeff*prod(variables_jump["$var"]^(exp[1]) for (var,exp) in monome) for (monome,coeff) in polynome) <= ub)
    end
    polynome_obj = pb_poly_real.objective
    @NLobjective(m,Min,sum(coeff*prod(variables_jump["$var"]^(exp[1]) for (var,exp) in monome) for (monome,coeff) in polynome_obj))
    return m
end

###TEST
# V1 = Variable("VOLT_1",Complex)
# V2 = Variable("VOLT_2",Complex)
# p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
# p_ctr1 = abs2(V1)
# p_ctr2 = 3im * V1*conj(V2) + abs2(V1)
#
# problem_poly=Problem()
# add_variable!(problem_poly,V1)
# add_variable!(problem_poly,V2)
# add_constraint!(problem_poly, "ctr1", 0.95^2 << p_ctr1 << 1.05^2)
# add_constraint!(problem_poly, "ctr2", p_ctr2==0)
# set_objective!(problem_poly, p_obj)
# print(problem_poly)
# mysolver = IpoptSolver()
# m = get_JuMP_cartesian_model(problem_poly, mysolver)
# print(m)

"""
    get_JuMP_polar_model(pb::Problem, mysolver)

Return JuMP polar model associated to `mysolver` defined by `pb`, a polynomial problem in complex numbers

# Example
```jldoctest
V1 = Variable("VOLT_1",Complex)
V2 = Variable("VOLT_2",Complex)
p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
p_ctr1 = abs2(V1)
p_ctr2 = Polynomial((3+4im) * V1*conj(V2))
pb=Problem()
add_variable!(pb,V1)
add_variable!(pb,V2)
add_constraint!(pb, "ctr1", 0 << p_ctr1 << 100)
# add_constraint!(pb, "ctr2", p_ctr2==3+4im)
set_objective!(pb, p_obj)
print(pb)
mysolver = IpoptSolver()
m = get_JuMP_polar_model(pb, mysolver)
print(m)
```
"""
function get_JuMP_polar_model(pb::Problem, mysolver)
    m = Model(solver = mysolver)
    jump_vars = Dict{String, JuMP.Variable}()
    for (varname, vartype) in pb.variables
        mod = "ρ_$varname"
        jump_vars["ρ_$varname"] = @variable(m, basename="ρ_$varname")
        jump_vars["θ_$varname"] = @variable(m, basename="θ_$varname")
    end
    ctr_jump = Dict{String,JuMP.ConstraintRef}()
    for (ctr, modeler_ctr) in pb.constraints
        p = modeler_ctr.p
        lb = modeler_ctr.lb
        ub = modeler_ctr.ub
        ps = Dict{Exponent, Tuple{Real, Real}}(exp=>(real(coeff), imag(coeff)) for (exp, coeff) in p)
        ## Real part of constraint
        @NLconstraint(m, real(lb) <= sum( coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        - coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= real(ub))

        ## Imag part of constraint
        @NLconstraint(m, imag(lb) <= sum( coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        + coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= imag(ub))

    polynome_obj = pb.objective
    ps = Dict{Exponent, Tuple{Real, Real}}(exp=>(real(coeff), imag(coeff)) for (exp, coeff) in polynome_obj)
    lb_real, ub_real = real(lb), real(ub)
    @NLobjective(m, :Min, lb_real <= sum( coeff[1] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * cos(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        - coeff[2] * prod(jump_vars["ρ_$var"]^(exp[1]+exp[2]) for (var,exp) in mon) * sin(sum( (exp[1]-exp[2])*jump_vars["θ_$var"] for (var,exp) in mon))
                                        for (mon, coeff) in ps) <= ub_real)
    end
    return m
end

####TEST
# V1 = Variable("VOLT_1",Complex)
# V2 = Variable("VOLT_2",Complex)
# p_obj = V1*conj(V2) + (2+im)*abs2(V1) + 1+2im
# p_ctr1 = abs2(V1)
# p_ctr2 = Polynomial((3+4im) * V1*conj(V2))
# pb=Problem()
# add_variable!(pb,V1)
# add_variable!(pb,V2)
# add_constraint!(pb, "ctr1", 0 << p_ctr1 << 100)
# # add_constraint!(pb, "ctr2", p_ctr2==3+4im)
# set_objective!(pb, p_obj)
# print(pb)
# mysolver = IpoptSolver()
# m = get_JuMP_polar_model(pb, mysolver)
# print(m)

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using YaoLang
using YaoLang.Compiler
using ExprTools
using IRTools
using IRTools: IR
using YaoArrayRegister
using FFTW
using Test

@device function qft(n::Int)
    1 => H
    for k in 2:n
        @ctrl k 1 => shift(2π / 2^k)
    end

    if n > 1
        2:n => qft(n - 1)
    end
end

# h q[0];
# cu1(pi/2) q[1],q[0];
# h q[1];
# cu1(pi/4) q[2],q[0];
# cu1(pi/2) q[2],q[1];
# h q[2];
# cu1(pi/8) q[3],q[0];
# cu1(pi/4) q[3],q[1];
# cu1(pi/2) q[3],q[2];
# h q[3];

@device function qft4()
    1 => H
    @ctrl 2 1 => shift(π / 2)
    @ctrl 3 1 => shift(π / 4)
    @ctrl 4 1 => shift(π / 8)

    2 => H
    @ctrl 3 2 => shift(π / 2)
    @ctrl 4 2 => shift(π / 4)

    3 => H
    @ctrl 4 3 => shift(π / 2)

    4 => H
end

@device function hadamard()
    1 => H
end

@testset "example/qft" begin
    r = rand_state(4)
    state_vec = statevec(r)
    a = invorder!(copy(r)) |> qft(4)
    kv = ifft(state_vec) * sqrt(length(state_vec))
    @test statevec(a) ≈ kv

    circ = qft(4)
    ir = @code_yao qft4()
    @test ir.qasm_compatible == true
    @test ir.pure_quantum == true
end

@testset "example/hadamard" begin
    ir = @code_yao hadamard()
    @test ir.pure_quantum == true
    @test ir.qasm_compatible == true
    r = rand_state(1)
    @test (copy(r) |> hadamard()) ≈ (copy(r) |> H())
end

@device function check_return(k::Int)
    c = @measure k
    return c
end

@testset "measure" begin
    r = ArrayReg(bit"000") + ArrayReg(bit"010") + ArrayReg(bit"110")
    @test bit"0" == check_return(1)(r)
end

@testset "printing" begin
    ir = @code_yao qft(3)
    println(ir)
end

@device function single_qubit_layer(configs::AbstractVector, gateset = (X, Y, Z))
    @inbounds for k in 1:size(configs, 1)
        idx = configs[k]
        if idx > 0
            k => gateset[idx]
        end
    end
end

@testset "default value (#19)" begin
    configs = rand(1:3, 4)
    gateset = (X, Y, Z)
    c = single_qubit_layer(configs)
    println(@code_yao single_qubit_layer(configs))
    r = rand_state(4)
    r1 = copy(r) |> c
    r2 = copy(r)
    for k in 1:4
        r2 |> gateset[configs[k]](Locations(k))
    end
    @test r1 ≈ r2
end

struct Foo
    a::Int
    b::Int
end

@device function (b::Foo)(theta)
    @ctrl b.a b.b => shift(theta)
    theta
end

@testset "callable" begin
    m = Foo(1, 2)
    r = rand_state(4)
    r1 = copy(r) |> m(1.2)
    r2 = copy(r) |> shift(1.2)(1, 2)
    @test r1 ≈ r2
end

@testset "purify" begin
    pure_qft = @quantum 3 qft(3)
    @test is_pure_quantum(@code_yao pure_qft())
end


text:
using ConformalMaps
using Test
using Requires

vertices = [1.0  0.0;
            0.0  1.0;
            -1.0  0.0;
            0.0 -1.0]
f = ConformalMap(vertices,0.0)
@test f(0.0) == 0

vertices = [1, im, -1, -im]
f = ConformalMap(vertices,0.0im)
g = inv(f)
z = 0.1 + 0.1im;
@test g(f(z)) ≈ z

vertices = Complex{BigFloat}[1, im, -1, -im]
f = ConformalMap(vertices,0.0im)
g = inv(f)
z = 0.1 + 0.1im;
@test g(f(z)) ≈ z

Requires.@require AsyPlots="77e5a97a-5ef9-58df-9d21-21957d92d960" begin
    slitdomain = [0.0, 0.495, 0.5 + 0.25*im, 0.505, 1.0, 1.0 + 1.0*im, im]
    f = ConformalMap(slitdomain,0.5+0.5im;resolution=40)
    g = inv(f)
    visualize(g)
end


text:
module MRIBloch

	using LinearAlgebra

	include("rotation_matrices.jl")

	# TODO: MAKE GENERATED?

	function simulate(mode::Val{:pulse_frame}, m0::NTuple{3, <: Real}, ω1::AbstractVector{<: Real}, Δω0::AbstractVector{<: Real}, dt::Real)
		# Memory for magnetisation
		m = Array{Float64}(undef, 3, 1+length(Δω0))
		m[:, 1] .= m0
		# Precompute flip angles
		α = @. dt * sqrt(Δω0^2 + ω1^2)
		# Run
		R = Matrix{Float64}(undef, 3, 3)
		for t = 1:length(α)
			# Compute rotation axis and matrix
			axis = (ω1[t], -Δω0[t])
			axis = axis ./ norm(axis)
			rotation_matrix(Val(:xz), axis, α[t]; out=R)
			# Apply
			@views m[:, t+1] = R * m[:, t]
		end
		return m
	end

	#function simulate(mode::Val{:B0_frame}, m0::NTuple{3, <: Real}, ω1::AbstractVector{<: Real}, ϕ::AbstractVector{<: Real}, dt::Real)
	#	# Memory for magnetisation
	#	m = Array{Float64}(undef, 3, 1+length(ϕ))
	#	m[:, 1] .= m0
	#	# Run
	#	R = Matrix{Float64}(undef, 3, 3)
	#	for t = 1:length(α)
	#		# Compute rotation axis and matrix
	#		α = dt * ω1[t]
	#		rotation_matrix(Val(:xy), ϕ[t], α; out=R)
	#		# Apply
	#		@views m[:, t+1] = R * m[:, t]
	#	end
	#	return m
	#end

	function simulate(mode::Val{:B0_frame}, m0::NTuple{3, <: Real}, ω1::AbstractVector{<: Complex}, dt::Real)
		# Memory for magnetisation
		m = Array{Float64}(undef, 3, 1+length(ω1))
		m[:, 1] .= m0
		# Run
		R = Matrix{Float64}(undef, 3, 3)
		for t = 1:length(ω1)
			# Compute rotation axis and matrix
			abs_ω1 = abs(ω1[t])
			α = dt * abs_ω1
			axis = ω1[t] / abs_ω1
			rotation_matrix(Val(:xy), reim(axis), α; out=R)
			# Apply
			@views m[:, t+1] = R * m[:, t]
		end
		return m
	end
end


text:
function mandelbrot(a)
    z = 0
    for i=1:50
        z = z^2 + a
    end
    return z
end

for y=1.0:-0.05:-1.0
    for x=-2.0:0.0315:0.5
        abs(mandelbrot(complex(x, y))) < 2 ? print("*") : print(" ")
    end
    println()
end

# Taken from: https://rosettacode.org/wiki/Mandelbrot_set#Julia


text:
using CUDA, Dates, Glob, SeisNoise, Serialization
include("MPI-CC-impl.jl")

# files and such
BASEDIR = "/media/four/"
FFTDIR = joinpath(BASEDIR,"LASSO/FFT/")
CORRDIR = joinpath(BASEDIR,"LASSO/CORR/")
if !isdir(CORRDIR)
    mkpath(CORRDIR)
end
files = glob("*",FFTDIR)
stas = [split(basename(s),'.')[2] |> (y -> parse(Int,y)) for s in files]
ind = sortperm(stas)
files = files[ind]

# create directory structure for saving CORRs
for sta in stas
    mkpath(joinpath(CORRDIR,"2A.$sta..DPZ"))
end

# parameters
maxlag = 200.
Nper = 18

# warm up functions
splits = splitFFT(files,Nper)
splitCC(splits[1:3],maxlag,CORRDIR)
testfiles = glob("*/*",CORRDIR)
rm.(testfiles)

t1 = now()
## xcorr on the GPU minimizing I/O
splitCC(splits,maxlag,CORRDIR)
t2 = now()
total = Dates.canonicalize(Dates.CompoundPeriod(t2-t1))
println("Total computation took $total")


text:
"""si basa sul file voxellizationProva ma usa il mio package invece di LinearAlgebraicAlgebra
e verrà usato per l'analisi delle prestazioni una volta implementata la parallellizazzione"""

using CuboidGrids, SparseArrays, DataStructures
using Distributed
Lar = CuboidGrids

using ViewerGL
GL = ViewerGL

# Dati da larGrid

function CV2FV( v::Array{Int64} )
	faces = [
		[v[1], v[2], v[3], v[4]], [v[5], v[6], v[7], v[8]],
		[v[1], v[2], v[5], v[6]],	[v[3], v[4], v[7], v[8]],
		[v[1], v[3], v[5], v[7]], [v[2], v[4], v[6], v[8]]]
end

function CV2EV( v::Array{Int64} )
	edges = [
		[v[1], v[2]], [v[3], v[4]], [v[5], v[6]], [v[7], v[8]], [v[1], v[3]], [v[2], v[4]],
		[v[5], v[7]], [v[6], v[8]], [v[1], v[5]], [v[2], v[6]], [v[3], v[7]], [v[4], v[8]]]
end

function K( CV )	#costruisce e ritorna una matrice sparsa a partire dalla collezione passata(matrice dei vertici, degli edges, delle facce, etc)
	I = vcat( [ [k for h in CV[k]] for k=1:length(CV) ]...) #righe
	J = vcat(CV...)											#colonne
	Vals = [1 for k=1:length(I)]							#valori della matrice
	return sparse(I,J,Vals)									#costruzione della matrice sparsa
end

V,CV = Lar.cuboidGrid([10,20,10])
#pointcloud,V,CV = random3cells([40,20,10],4_000)

VV = [[v] for v=1:size(V,2)]
FV = convert(Array{Array{Int64,1},1}, collect(Set(cat(pmap(CV2FV,CV))))) #verts faces
EV = convert(Array{Array{Int64,1},1}, collect(Set(cat(pmap(CV2EV,CV))))) #verts edges

#Cotruzine delle varie matrici tramite la funzione K
M_0 = K(VV)
M_1 = K(EV)
M_2 = K(FV)
M_3 = K(CV)

∂_1 = M_0 * M_1'
∂_2 = (M_1 * M_2') .÷ 2 #	.÷ sum(M_1,dims=2)
s = sum(M_2,dims=2)
∂_3 = (M_2 * M_3')
∂_3 = ∂_3 ./	s
∂_3 = ∂_3 .÷ 1	#	.÷ sum(M_2,dims=2)

S2 = sum(∂_3,dims=2)
inner = [k for k=1:length(S2) if S2[k]==2]
outer = setdiff(collect(1:length(FV)), inner)

# point cloud
#GL.VIEW([GL.GLPoints(convert(Lar.Points,pointcloud'), GL.COLORS[1])])

# voxels
GL.VIEW([GL.GLLar2gl(V,CV)])

#=
# edges
GL.VIEW([ GL.GLGrid(V,EV,GL.Point4d(1,1,1,1))
          GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1)) ]);

# inner faces
GL.VIEW([ GL.GLGrid(V,FV[inner],GL.Point4d(1,1,1,1))
          GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1)) ]);

# outer faces
GL.VIEW([ GL.GLGrid(V,FV[outer],GL.Point4d(1,1,1,1))
          GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1)) ]);

=#

#-------------------------------------------------------------------------------------------

V,CV = Lar.cuboidGrid([5,20,10])
#pointcloud,V,CV = random3cells([40,20,10],4_000)

VV = [[v] for v=1:size(V,2)]
FV = convert(Array{Array{Int64,1},1}, collect(Set(cat(pmap(CV2FV,CV)))))
EV = convert(Array{Array{Int64,1},1}, collect(Set(cat(pmap(CV2EV,CV)))))

M_0 = K(VV)
M_1 = K(EV)
M_2 = K(FV)
M_3 = K(CV)

∂_1 = M_0 * M_1'
∂_2 = (M_1 * M_2') .÷ 2 #	.÷ sum(M_1,dims=2)
s = sum(M_2,dims=2)
∂_3 = (M_2 * M_3')
∂_3 = ∂_3 ./	s
∂_3 = ∂_3 .÷ 1	#	.÷ sum(M_2,dims=2)

S2 = sum(∂_3,dims=2)
inner = [k for k=1:length(S2) if S2[k]==2]
outer = setdiff(collect(1:length(FV)), inner)

# point cloud
#GL.VIEW([GL.GLPoints(convert(Lar.Points,pointcloud'), GL.COLORS[1])])

# voxels
GL.VIEW([GL.GLLar2gl(V,CV)])


text:
# -----------------------------------------------------------
#
# Implementing a nodal basis for DG sparse grids
#
# -----------------------------------------------------------

#=

We begin by defining analogues of our DG functions
to satisfy the following conditions:

    1) There are a set of collocation points respecting the hierarchical structure
    2) The nodal basis functions at a given level are those that vanish
        on all nodal points at the above levels and vanish on all but one
        of the cells on the level below
    3) Third constraint?

To do this, we _restrict_ to the continuous part of the DG basis, and focus
on piecewise polynomials that vanish on their boundaries. We then pick the
collocation points to be exactly at these boundaries.

For this, it is useful to define lagrange polynomials, and moreover
use k=2, 3, 5.
The k=2 case has our basis is exactly the hat function basis constructed
previously.

For k=2, the coarsest level has nodal points at {0, 1},
    the next level (the first hat) has nodal point at 0.5
    the next level has two cells with points at 0.25, 0.75 etc.

For k=3 the coarsest level has nodal points at {0, 0.5, 1},
    the next level has two nodal functions with points at 0.25, 0.75,
    the next level has two cells and two nodes per cell at all odd multiples of 0.125 etc.

For k=5 the coarsest level has nodal points at {0, 0.25, 0.5, 0.75, 1}
    the next level has four nodal functions with points at all multiples of 0.125
    the next level has two cells and four nodes per cell at all odd multiples of 0.0625 etc.

Notice that k=2,3,5 work because at each cell of the multiresolution basis the space of
    picewise continuous polynomials is k-1 = 2^i for i=0,1,2 in this case.
    A construction for other k would be harder, if not impossible.

=#

# Lagrange polynomials interpolating at the coarsest level
function lag_nodal(k::Int, mode::Int, x::Real)
    (k < 2 || k > 5) && throw(ArgumentError("k is not between 2 and 5"))
    (mode > k || mode < 1) && throw(ArgumentError("mode is not in [1, k]"))

    (x < 0 || x > 1) && return zero(x)
    if k == 2
        mode == 1 && return x
        mode == 2 && return 1 - x
    elseif k == 3
        mode == 1 && return 2 * (x - 1//2) * (x - 1)
        mode == 2 && return -4 * x * (x - 1)
        mode == 3 && return 2 * x * (x - 1//2)
    elseif k == 4
        throw(ArgumentError("k = 4 not implemented"))
    elseif k == 5
        mode == 1 && return 32//3 * (x - 1//4) * (x - 1//2) * (x - 3//4) * (x - 1)
        mode == 2 && return -128//3 * x * (x - 1//2) * (x - 3//4) * (x - 1)
        mode == 3 && return 64 * x * (x - 1//4) * (x - 3//4) * (x - 1)
        mode == 4 && return -128//3 * x * (x - 1//4) * (x - 1//2) * (x - 1)
        mode == 5 && return 32//3 * x * (x - 1//4) * (x - 1//2) * (x - 3//4)
    end
end

# Multiresolution functions for the levels below
function h_nodal(k::Int, mode::Int, x::Real)
    (k < 2 || k > 5) && throw(ArgumentError("k is not between 2 and 5"))
    (mode > k || mode < 1) && throw(ArgumentError("mode is not in [1, k]"))

    if k == 2
        mode == 1 && return max(zero(x), 1-abs(x))
        if mode == 2
            (abs(x) > 1 || abs(x) == 0) && return zero(x)
            x < 0 && return (x + 1)
            x > 0 && return zero(x)# (x - 1)
        end
    elseif k == 3
        mode == 1 && return max(zero(x), -4 * x * (x + 1))
        mode == 2 && return max(zero(x), -4 * x * (x - 1))
        if mode == 3
            (abs(x) > 1 || abs(x) == 0) && return zero(x)
            x < 0 && return 2 * (x + 1) * (x + 1//2)
            x > 0 && return zero(x) # -2 * (x - 1) * (x - 1//2)
        end
    elseif k == 4
        throw(MethodError("k = 4 not implemented"))
    elseif k == 5
        mode == 1 && return (x < -1 || x > 0) ? zero(x) : -128//3 * x * (x + 1) * (x + 1//2) * (x + 1//4)
        mode == 2 && return (x < -1 || x > 0) ? zero(x) : -128//3 * x * (x + 1) * (x + 1//2) * (x + 3//4)
        mode == 3 && return (x > 1  || x < 0) ? zero(x) : -128//3 * x * (x - 1) * (x - 1//2) * (x - 3//4)
        mode == 4 && return (x > 1  || x < 0) ? zero(x) : -128//3 * x * (x - 1) * (x - 1//2) * (x - 1//4)
        if mode == 5
            (abs(x) > 1 || abs(x) == 0) && return zero(x)
            x < 0 && return 32//3 * (x + 1) * (x + 3//4) * (x + 1//2) * (x + 1//4)
            x > 0 && return zero(x) #32//3 * (x - 1) * (x - 3//4) * (x - 1//2) * (x - 1//4)
        end
    end
end

# This yields our nodal basis:
function v_nodal(k::Int, level::Int, cell::Int, mode::Int, x::Real)
    if level==0
        return lag_nodal(k, mode, x)
    else
        return h_nodal(k, mode, (1<<level)*x - (2*cell-1))
    end
end

function v_nodal(k::Int, level::Int, cell::Int, mode::Int)
    return x->v_nodal(k, level, cell, mode, x)
end

# -----------------------------------------------------------
# Made for evaluating the dicontinuous elements in the nodal
# basis at their points of discontinuity
# -----------------------------------------------------------

# Evaluates the modal basis at a given point, averaging over jumps
function eval_v_periodic(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = max(x - 1e-16, 0)
    xr = min(x + 1e-16, 1)
    (xl == 0) && (xl = 1 - 1e-16)
    (xr == 1) && (xr = 1e-16)
    return 1//2 * (v(k, level, cell, mode, xl) + v(k, level, cell, mode, xr))
end

# Find the size of the discontinuity of the modal basis at a given point
function gap_v(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = max(x - 1e-16, 0)
    xr = min(x + 1e-16, 1)
    return 1//2 * (v(k, level, cell, mode, xr) - v(k, level, cell, mode, xl))
end

# Evaluates the modal basis at a given point, averaging over jumps
function eval_v_nodal_periodic(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = max(x - 1e-30, 0)
    xr = min(x + 1e-30, 1)
    (xl == 0) && (xl = 1 - 1e-30)
    (xr == 1) && (xr = 1e-30)
    return 1//2 * (v_nodal(k, level, cell, mode, xl) + v_nodal(k, level, cell, mode, xr))
end

# Find the size of the discontinuity of the modal basis at a given point
function avg_v_nodal(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = x - 1e-16
    xr = x + 1e-16
    (xl < 0) && (xl = xr)
    (xr > 1) && (xr = xl)
    return 1//2 * (v_nodal(k, level, cell, mode, xr) + v_nodal(k, level, cell, mode, xl))
end

# Find the size of the discontinuity of the modal basis at a given point
function gap_v_nodal(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = x - 1e-16
    xr = x + 1e-16
    (xl < 0) && (xl = xr)
    (xr > 1) && (xr = xl)
    return 1//2 * (v_nodal(k, level, cell, mode, xr) - v_nodal(k, level, cell, mode, xl))
end

function eval_v_nodal_right(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xr = x + 1e-16
    return v_nodal(k, level, cell, mode, xr)
end

function eval_v_nodal_left(k::Int, level::Int, cell::Int, mode::Int, x::T) where {T<:Real}
    xl = x - 1e-16
    return v_nodal(k, level, cell, mode, xl)
end

# -----------------------------------------------------------
# Next we define the sparse matrices transforming between
# the heirarchical, nodal, and point bases
# -----------------------------------------------------------

# Huge TODO:
#     Replace all of these for loops and multiple function calls
#     by cleverly using the iterable interface in julia

# Evaluation of a given nodal function along all relevant gridpoints
function eval_points_1D(k::Int, max_level::Int, level::Int, cell::Int, mode::Int)
    i = 1
    I = Int[]
    V = Float64[]
    for l in 0:max_level
        l == 0 ? (node_range = 0:(1//(k-1)):(1 - 1//(k-1))) :
                 (node_range = 1//(2*(k-1)):1//(k-1):(1-1//(2*(k-1))))

        for c in 1:1<<max(0, l-1)
            for node in node_range
                x = (c - 1 + node)//(1<<max(0, l-1))
                val = eval_v_nodal_right(k, level, cell, mode, x)
                (val != zero(x)) && (push!(I, i); push!(V, val))
                i += 1
            end

            if l == 0
                x = 1
                val = eval_v_nodal_left(k, level, cell, mode, x)
            else
                x = (c - 1//2)/(1<<max(0, l-1))
                val = eval_v_nodal_left(k, level, cell, mode, x)
            end
            (val != zero(x)) && (push!(I, i); push!(V, val))
            i += 1
        end
    end
    return sparsevec(I, V)
end

# Sparse matrix converting nodal basis element to its value
# on the collocation points
function nodal2points_1D(k::Int, max_level::Int; atol=1e-15)
    I = Int[]
    J = Int[]
    V = Float64[]
    j = 1

    for level in 0:max_level
        mode_range = 1:k

        for cell in 1:1<<max(0, level-1)
            for mode in mode_range
                coeffs = eval_points_1D(k, max_level, level, cell, mode)
                for (i, val) in zip(findnz(coeffs)...)
                    push!(I, i)
                    push!(J, j)
                    push!(V, val)
                end
                j += 1
            end
        end
    end
    return threshold(sparse(I, J, V), atol)
end

# Sparse matrix inverting the above construction
# Converting a set of values on the collocation points to the DG space
# in terms of the nodal basis
function points2nodal_1D(k::Int, max_level::Int; atol=1e-15)
    n2p = nodal2points_1D(k, max_level)
    return threshold(inv(Matrix(n2p)), atol)
end

# Converts a given nodal basis element to a sum in the standard
# hieararchical DG scheme
# TODO: Remove hquadrature and do this directly
function nodal2modal_1D(k::Int, level::Int, cell::Int, mode::Int; rtol=1e-10, atol=1e-15, max_evals=1500)
    I = Int[]
    V = Float64[]
    i = 1

    for hier_level in 0:level
        hier_cells = 1<<max(0, hier_level-1)
        for hier_cell in 1:hier_cells
            for hier_mode in 1:k
                x_min = (hier_cell-1)//hier_cells
                x_max = (hier_cell)//hier_cells
                val = hquadrature(x->v_nodal(k, level, cell, mode, x)*
                                    v(k, hier_level, hier_cell, hier_mode, x),
                                    x_min, x_max,
                                    reltol=rtol, abstol=atol, maxevals=max_evals)[1]
                (abs(val) > eps(atol)) && (push!(I, i); push!(V, val))
                i += 1
            end
        end
    end
    return dropzeros!(sparsevec(I, V))
end

function nodal2modal_1D(k::Int, max_level::Int; rtol=1e-10, atol=1e-15, max_evals=1500)
    I = Int[]
    J = Int[]
    V = Float64[]
    j = 1

    for nodal_level in 0:max_level
        num_nodes = k
        nodal_cells = 1<<max(0, nodal_level-1)
        for nodal_cell in 1:nodal_cells
            for nodal_mode in 1:num_nodes
                coeffs = nodal2modal_1D(k, nodal_level, nodal_cell, nodal_mode;
                                        rtol=rtol,
                                        atol=atol,
                                        max_evals=max_evals)
                for (i, val) in zip(findnz(coeffs)...)
                    push!(I, i)
                    push!(J, j)
                    push!(V, val)
                end
                j += 1
            end
        end
    end
    return threshold(sparse(I, J, V), atol)
end

# Accounting for ambiguity in the definition of the hierarchical
# basis across regions
function eval_v(k, level, cell, mode, x)
    xl = max(x - 1e-16, 0)
    xr = min(x + 1e-16, 1)
    (xl == 0) && (xl = 1 - 1e-16)
    (xr == 1) && (xr = 1e-16)
    return 0.5 * (v(k, level, cell, mode, xl) + v(k, level, cell, mode, xr))
end

function nodal2pos_1D(k, max_level, nodal_level, nodal_cell, nodal_mode;
                   rtol=1e-10, atol=1e-15, max_evals=1500)
    I = Int[]
    V = Float64[]
    i = 1

    nodal_min = (nodal_cell-1)/(1<<max(0, nodal_level-1))
    nodal_med = (nodal_cell-1//2)/(1<<max(0, nodal_level-1))
    nodal_max = (nodal_cell)/(1<<max(0, nodal_level-1))
    for c in 1:1<<max_level
        pos_min = (c-1)/(1<<max_level)
        pos_med = (c-1//2)/(1<<max_level)
        pos_max = c/(1<<max_level)

        sm = 1e-20

        if (pos_min > nodal_max) || (pos_max < nodal_min)
            i += k
            continue
        end
        for m in 1:k
            val1  = hquadrature(x->v_nodal(k, nodal_level, nodal_cell, nodal_mode, x)*basis(max_level, c, m, x),
                              pos_min + sm, pos_med - sm,
                              rtol=rtol, atol=atol, maxevals=max_evals)[1]

            val2 = hquadrature(x->v_nodal(k, nodal_level, nodal_cell, nodal_mode, x)*basis(max_level, c, m, x),
                              pos_med + sm, pos_max - sm,
                              rtol=rtol, atol=atol, maxevals=max_evals)[1]

            val = val1 + val2
            (abs(val) > eps(atol)) && (push!(I, i); push!(V, val))
            i += 1
        end
    end
    return sparsevec(I, V)
end

function nodal2pos_1D(k, max_level; rtol=1e-10, atol=1e-15, max_evals=1500)
    I = Int[]
    J = Int[]
    V = Float64[]
    j = 1

    for nodal_level in 0:max_level
        num_nodes = k
        nodal_cells = 1<<max(0, nodal_level-1)
        for nodal_cell in 1:nodal_cells
            for nodal_mode in 1:num_nodes
                coeffs = nodal2pos_1D(k, max_level, nodal_level, nodal_cell, nodal_mode;
                                        rtol=rtol,
                                        atol=atol,
                                        max_evals=max_evals)
                for (i, val) in zip(findnz(coeffs)...)
                    push!(I, i)
                    push!(J, j)
                    push!(V, val)
                end
                j += 1
            end
        end
    end
    return threshold(sparse(I, J, V), atol)
end

function pos2nodal_1D(k, max_level; rtol=1e-10, atol=1e-12, max_evals=1500)
    n2pos = nodal2pos_1D(k, max_level; rtol=rtol, atol=atol, max_evals=max_evals)
    return threshold(inv(Matrix(n2pos)), atol)
end

# Final function, combining all of the above:
function transform_1D(k::Int, n::Int, from::String, to::String; atol=1e-12)
    if from == "nodal" && to == "pos"
        return nodal2pos_1D(k, n; atol=atol)
    elseif from == "pos" && to == "nodal"
        return pos2nodal_1D(k, n; atol=atol)
    elseif from == "nodal" && to == "points"
        return nodal2points_1D(k, n; atol=atol)
    elseif from == "points" && to == "nodal"
        return points2nodal_1D(k, n; atol=atol)
    elseif from == "nodal" && to == "modal"
        return threshold(pos2hier(k, n) * nodal2pos_1D(k, n), atol)
    elseif from == "modal" && to == "nodal"
        return threshold(pos2nodal_1D(k, n) * hier2pos(k, n), atol)
    elseif from == "modal" && to == "pos"
        return hier2pos(k, n; atol=atol)
    elseif from == "pos" && to == "modal"
        return pos2hier(k, n; atol=atol)
    elseif from == "points" && to == "pos"
        return threshold(nodal2pos_1D(k, n) * points2nodal_1D(k, n), atol)
    elseif from == "pos" && to == "points"
        return threshold(nodal2points_1D(k, n) * pos2nodal_1D(k, n), atol)
    elseif from == "modal" && to == "points"
        m2n = transform_1D(k, n, "modal", "nodal"; atol=atol)
        n2p = transform_1D(k, n, "nodal", "points"; atol=atol)
        return threshold(n2p * m2n, atol)
    elseif from == "points" && to == "modal"
        p2n = transform_1D(k, n, "points", "nodal"; atol=atol)
        n2m = transform_1D(k, n, "nodal", "modal"; atol=atol)
        return threshold(n2m * p2n, atol)
    else throw(ArgumentError("A basis was undefined"))
    end
end


text:
export PrecomputedInitialisation, RandomInitialisation, KMeansInitialisation, IncrementalInitialisation

mutable struct PrecomputedInitialisation <: InitialisationType
  Z::Array{Int}
  PrecomputedInitialisation(Z::Array{Int}) = new(Z)
end

mutable struct RandomInitialisation <: InitialisationType
  k::Int
  RandomInitialisation(;k = 2) = new(k)
end

mutable struct KMeansInitialisation <: InitialisationType
  k::Int
  maxiterations::Int
  KMeansInitialisation(;k = 2, maxiterations = 1000) = new(k, maxiterations)
end

mutable struct IncrementalInitialisation <: InitialisationType end


text:
using Dash
function todo_app(content_callback)
    app = dash()

    content = html_div() do
        html_div("Dash To-Do list"),
        dcc_input(id="new-item"),
        html_button("Add", id="add"),
        html_button("Clear Done", id="clear-done"),
        html_div(id="list-container"),
        html_hr(),
        html_div(id="totals")
    end

    if content_callback
        app.layout = html_div() do
            html_div(id="content"),
            dcc_location(id="url")
        end
        callback!(app, Output("content", "children"), Input("url", "pathname")) do v
            return content
        end
    else
        app.layout = content
    end

    style_todo = Dict("display" => "inline", "margin" =>"10px")
    style_done = Dict("textDecoration" => "line-through", "color" => "#888")
    merge!(style_done, style_todo)

    callback!(app, Output("list-container", "children"), Output("new-item", "value"),
            Input("add", "n_clicks"),
            Input("new-item", "n_submit"),
            Input("clear-done", "n_clicks"),
            State("new-item", "value"),
            State((item = ALL,), "children"),
            State((item = ALL, action = "done"), "value"),
    ) do add, add2, clear, new_item, items, items_done
        triggered = [t.prop_id for t in callback_context().triggered]
        adding = any((id)->id in ("add.n_clicks", "new-item.n_submit"), triggered)
        clearing = any((id)->id == "clear-done.n_clicks", triggered)

        new_spec = Vector{Any}()

        push!.(Ref(new_spec),
            Iterators.filter(zip(items, items_done)) do (text, done)
                return !(clearing && !isempty(done))
            end
        )

        adding && push!(new_spec, (new_item, []))
        new_list = [
            html_div(style = (clear = "both",)) do
                dcc_checklist(
                    id=(item = i, action = "done"),
                    options=[(label = "", value = "done")],
                    value=done,
                    style=(display = "inline",),
                ),
                html_div(
                    text, id=(item = i,), style= !isempty(done) ? style_done : style_todo
                ),
                html_div(id=(item = i, preceding = true), style=style_todo)
            end
            for (i, (text, done)) in enumerate(new_spec)
        ]
        return [new_list, adding ? "" : new_item]
    end

    callback!(app,
        Output((item = MATCH,), "style"),
        Input((item = MATCH, action = "done"), "value")
    ) do done
        println(done)
        return !isempty(done) ? style_done : style_todo
    end

    callback!(app,
        Output((item = MATCH, preceding = true), "children"),
        Input((item = ALLSMALLER, action = "done"), "value"),
        Input((item = MATCH, action = "done"), "value"),

    ) do done_before, this_done
        !isempty(this_done) && return ""
        all_before = length(done_before)
        done_before = count((d)->!isempty(d), done_before)
        out = "$(done_before) of $(all_before) preceding items are done"
        (all_before == done_before) && (out *= " DO THIS NEXT!")
        return out
    end

    callback!(app,
        Output("totals", "children"), Input((item = ALL, action = "done"), "value")
    ) do done
        count_all = length(done)
        count_done = count((d)->!isempty(d), done)
        result = "$(count_done) of $(count_all) items completed"
        (count_all > 0) && (result *= " - $(floor(Int, 100 * count_done / count_all))%")

        return result
    end

    return app

end


text:
#this file holds functions related to data processing

#reads in a CSV type
function readSingleTypeCSV(path::String; colType=Union{Missing, String},
    streamer::Function = I_CSV_STREAMER)::DataFrame

  #count the number of cols
  iStream = streamer(path)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)
  nCols::Int = size(df,2)

  #set the type
  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  df = CSV.read(iStream, strings=:raw, categorical=false, types=types)
  close(iStream)

  return df
end

#same as the read version only returns nothing and just checks for errors
function validateSingleTypeCSV(path::String; colType=Union{Missing, String},
  streamer::Function = I_CSV_STREAMER, missingStrings=MISSING_STRINGS)::Nothing

  #count the number of cols

  iStream = streamer(path)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)
  nCols::Int = size(df,2)

  #set the type

  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  CSV.validate(iStream, strings=:raw, categorical=false, types=types,
    missingstrings=missingStrings, strict=false)
  close(iStream)

  return nothing

end

#this fucntion is almsot a script which analyzes data
function preProcessData(; refreshData::Bool=true, validateDownloads::Bool = true,
  nccsPath::String = NCCS_PATH,  refreshMetadata::Bool = true, reserializeData::Bool = true,
  refreshSummary = true)::Vector{MetaDictionary}

  downloadNCCSCSVs(NCCS_PATH, refreshData=refreshData, #=delay=0.001,=# validate=validateDownloads)
  metadata::Vector{MetaDictionary} = getMetadata(refreshMetadata = refreshMetadata)
  loadNCCSFiles(metadata, reserializeData=reserializeData, serializeOnly = true, mode=:parallel)
  consolidateSummary(metadata, refreshSummary=refreshSummary)

  #=println("fileName: $(metadata[1].fileName)")
  willFail(metadata[1])=#

  return metadata
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
__precompile__()

module PlotlyWebIO

export WebIOPlot, render

using Reexport
@reexport using PlotlyBase
using WebIO
using JSON
using JSExpr
using JSExpr: @var, @new

# need to import some functions because methods are meta-generated
import PlotlyBase:
    restyle!, relayout!, update!, addtraces!, deletetraces!, movetraces!,
    redraw!, extendtraces!, prependtraces!, purge!, to_image, download_image,
    restyle, relayout, update, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces, prep_kwargs, sizes, savefig, _tovec

include("plot_type.jl")
include("plotly_api.jl")

end # module


text:
function geomgoogle(json_response::Dict{String, Any})
    gem = json_response["results"][1]["geometry"]["location"]
    return((gem))
end

function geomOSM(json_response::Dict{String, Any})
    gem = json_response["results"][1]["locations"][1]["latLng"]
    return((gem))
end

"""
Print all the longitudes from geocode output.

    Usage:
    `julia> longitude(geocode_output)`

Note: This function can only be used if there were more than one addresses in `geocode()`.
"""

function longitude(geocode_output)
    lng = getindex.(geocode_output, "lng")
    return(lng)
end

"""
Print all the latitudes from geocode output.

    Usage:
    `julia> latitude(geocode_output)`

Note: This function can only be used if there were more than one addresses in `geocode()`.
"""

function latitude(geocode_output)
    lat = getindex.(geocode_output, "lat")
    return(lat)
end


text:
export deep_ope_d4rl_evaluate

using ReinforcementLearningBase
using ReinforcementLearningEnvironments
using PyCall
using UnicodePlots
using Random
using ProgressMeter

"""
    deep_ope_d4rl_evaluate(env_name, agent, epoch; <keyword arguments>)

Return the `UnicodePlot` for the `env_name`, `agent`, `epoch` that is given. Provide `gym_env_name` for specifying the environment explicitly.
`γ` is the discount factor which defaults to 1. Seed of the env can be provided in `env_seed`.
"""
function deep_ope_d4rl_evaluate(
    env_name::String,
    agent::String,
    epoch::Int;
    gym_env_name::Union{String, Nothing}=nothing,
    rng::AbstractRNG=MersenneTwister(123),
    num_evaluations::Int=10,
    γ::Float64=1.0,
    noisy::Bool=false,
    env_seed::Union{Int, Nothing}=nothing
)
    policy_folder = "$(env_name)_$(agent)_$(epoch)"

    if gym_env_name === nothing
        for policy in D4RL_POLICIES
            policy_file = split(policy["policy_path"], "/")[end]
            if chop(policy_file, head=0, tail=4) == policy_folder
                gym_env_name = policy["task.task_names"][1]
                break
            end
        end

        if gym_env_name === nothing error("invalid parameters") end
    end

    env = GymEnv(gym_env_name; seed=env_seed)

    model = d4rl_policy(env_name, agent, epoch)
    scores = Vector{Float64}(undef, num_evaluations)

    @showprogress for eval in 1:num_evaluations
        score = 0
        reset!(env)
        while !is_terminated(env)
            s = state(env)
            a = model(s;rng=rng, noisy=noisy)[1]
            s, a , env(a)
            r = reward(env)
            t = is_terminated(env)
            score += r*γ*(1-t)
        end
        scores[eval] = score
    end
    plt = lineplot(1:length(scores), scores, title = "$(gym_env_name) scores", name = "scores", xlabel = "episode", canvas = DotCanvas, ylabel = "score", border=:ascii)
    plt
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ b7c63d26-0827-491f-908a-02c3c8bf5d0c
begin
	using BenchmarkTools, LinearAlgebra, OffsetArrays, SpecialFunctions, StaticArrays
	using Plots, Profile, ProfileSVG, PlutoUI
	using LoopVectorization
end

# ╔═╡ 7f1163a2-5a4b-4892-a5ed-cac70c41d3b7
md"# Standalone Quantum Chemistry Molecular Integrals"

# ╔═╡ e109b7eb-725d-499c-b136-9552fff9da38
PlutoUI.TableOfContents()

# ╔═╡ 927fecfa-bb26-11eb-3126-6393fad683fd
md"
The following file is the time-intensive part of the [MolecularIntegrals.jl]()
package. This is called thousands of times in a typical calculation, and every
bit of speed is important here.

I'd be grateful for any help people could give me on speeding things up. Thus far,
I have not had any luck speeding this up further using macros like `@simd`
or `@turbo`. Often one of these macros speeds up the timing of the `vrr!` routine, but slows down the overall code for [reasons I don't understand](https://discourse.julialang.org/t/turbo-speeds-routine-slows-down-everything-else/62163).

I've read through and checked everything in the [Julia Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/) page and @ChrisRackauckas' [7 Julia Gotchas and How to Handle Them](https://www.stochasticlifestyle.com/7-julia-gotchas-handle/) blog post.

I believe there are still performance gains to be had, because the head to head timing of my code against those C/C++ libraries shows that the Julia code is still 5-10 times slower.

The best algorithms in common use make use of recurrence relations to generate integrals for higher angular momentum basis functions in terms of the lower angular momentum integrals, which are ultimate worked on in terms of incomplete error functions. These methods involve a certain amount of irregular memory access.
"

# ╔═╡ 9282d572-29e6-4c48-8981-4280cf0de817
md"## vrr routine"

# ╔═╡ 94bcfcc1-5216-4a2d-a136-0c8d3dae675a
md"""
Here's the version that has @turbo turned on. It's faster for the tests here, but
for some reason slows down the full 2e integral code.
"""

# ╔═╡ d4aa62ea-7bd3-4094-ac77-1d50a97d3b6d
md"""
Here's a version that has @inbounds turned on. It has similar results to @turbo: ie faster for the single-routine behchmark, but slower on the full 2e integral code.
"""

# ╔═╡ 7fe4e843-7543-45ea-8765-ed04dfab157e
md"## Timing results
The standalone `vrr` routine is coming in at ~28 μs on my development box (old mac mini); the one with `@turbo` clocks at ~19 μs, the one with @inbounds is ~25 μs.
"

# ╔═╡ 1fb56ae0-654d-4288-9da7-1ac4869115a4
md"""## Profiling
Profiling is not working well from Pluto. For comparison, here's a recent screenshot of the full MolecularIntegrals.jl profile from vscode.
![MolecularIntegrals.jl profile](https://raw.githubusercontent.com/rpmuller/MolecularIntegrals.jl/master/profile-2021-05-22.png)
`vrr!` is toward the bottom. `Fgamma` takes up a big chunk, and much of the rest are standard Julia functions like `*` and `Base.math.exp` and `getindex`.
"""

# ╔═╡ 63626a39-616b-4264-932d-114e3e7aca30
ProfileSVG.view(maxdepth=65)

# ╔═╡ 92bf1847-1b11-4a2d-9231-ef629ccbf0f6
md"## Utility functions
"

# ╔═╡ d01acace-1209-425f-826a-5efadb27b947
@inline factorial2(n::Int64) = prod(n:-2:1); # double factorial !!

# ╔═╡ 9174abad-51b6-4dbe-b326-2ee25b56b8af
@inline dist2(dxyz) = dot(dxyz,dxyz);

# ╔═╡ 0b75d079-bdb6-48c2-bb7f-d3d257bc17f6
@inline dist2(xyz1,xyz2) = dist2(xyz1-xyz2);

# ╔═╡ 7111c63d-bbc7-4343-8941-0b87228a9694
"gammainc - return the lower incomplete gamma function"
@inline gammainc(a,x) = gamma(a)*gamma_inc(a,x)[1];

# ╔═╡ b4dde3a6-7ab4-4f4f-8270-ff32d8d3e12b
"boys_array_gamma - Compute the Boys Fm(T) function.
This is currently one of the rate determining steps of the program, but
it's a *separate* problem than the one I want to explore in this notebook."
function boys_array_gamma(mmax,T,SMALL=1e-18)
    T = max(T,SMALL) # needs underflow protection because of inverse
    boys_array = zeros(Float64,mmax)
    ooT = 1/T
    denom = sqrt(ooT)
    boys_array[1] = 0.5*denom*gammainc(0.5,T)
    for m in 2:mmax
        denom *= ooT
        # Could speed this up more by expressing gamma(m) in terms of gamma(m±1)
        boys_array[m] = 0.5*denom*gammainc(m-0.5,T)
    end
    return boys_array
end;

# ╔═╡ ce86e563-3dee-4c0c-9cc7-a3efbfc3e326
md"## Global data
I define global data, which I understand is discouraged,
but which are all defined `const`. Timing increases
when I inline these in `vrr`."

# ╔═╡ 828390fb-fe8c-4a8d-8882-642ea78e3555
const shell_indices = Dict(
	# map from a shell l-value to the Cartesian version of m-values that
	# are the powers of the Cartesian Gaussian basis functions.
	# Uses `StaticArrays` for some speed gains."
    0 => [MVector(0,0,0)], # 1
    1 => [MVector(1,0,0),MVector(0,1,0),MVector(0,0,1)], # 3
    2 => [MVector(2,0,0),MVector(1,1,0),MVector(1,0,1),MVector(0,2,0),
			MVector(0,1,1),MVector(0,0,2)],
    3 => [MVector(3,0,0),MVector(2,1,0),MVector(2,0,1),
            MVector(1,2,0),MVector(1,0,2),MVector(1,1,1),
            MVector(0,3,0),MVector(0,2,1),MVector(0,1,2),MVector(0,0,3)], # 10
    4 => [MVector(4,0,0),MVector(3,1,0),MVector(3,0,1),MVector(2,2,0),
			MVector(2,1,1),MVector(2,0,2),
            MVector(1,3,0),MVector(1,2,1),MVector(1,1,2),
			MVector(1,0,3),
            MVector(0,4,0),MVector(0,3,1),MVector(0,2,2),
			MVector(0,1,3),MVector(0,0,4)] # 15
);

# ╔═╡ 3b321740-b76a-4b81-896d-a328398e6477
"make_m2ao - Map between a sequential list of `(mx,my,mz)` values and ao indices"
function make_m2ao(lmax=4)
    m2ao = Dict{Vector{Int64}, Int64}()
    iao = 0
    for l in 0:lmax
        for ms in shell_indices[l]
            iao += 1
            m2ao[ms] = iao
        end
    end
    return m2ao
end;

# ╔═╡ f364fe16-17d3-42d0-9ad4-f562168a4a35
function make_ao2m(lmax=4)
    ao2m = MVector{3, Int64}[]
    for l in 0:lmax
        for ms in shell_indices[l]
            push!(ao2m,ms)
        end
    end
    return ao2m
end;

# ╔═╡ f6b8a5d1-18fd-41d3-bd2d-9eb85001af2c
"make_nao - Number of AOs for system with l shells"
make_nao(l) = sum(length(shell_indices[i]) for i in 0:l);

# ╔═╡ c6d84970-9d30-4214-a456-eb03db812aa9
const ao2m = make_ao2m();

# ╔═╡ ec367e77-8eb0-4e66-89a1-249841e1ecfe
function make_shift_direction()
    n = length(ao2m)
    shift_direction = zeros(Int,n)
    for a in 1:n
        shift_direction[a] = argmax(ao2m[a])
    end
    return shift_direction
end;

# ╔═╡ b9142292-e0e8-4bb7-98df-16ab5c23a8c6
const shift_direction = make_shift_direction();

# ╔═╡ 82f11fe4-e95d-4f3d-be43-e65b9c0b3058
function make_shell_number()
    n = length(ao2m)
    shell_number = zeros(Int,n)
    for a in 1:n
        shell_number[a] = sum(ao2m[a])
    end
    return shell_number
end;

# ╔═╡ 58071e7c-13d3-406d-9e93-4ca5e0906b6e
const shell_number = make_shell_number();

# ╔═╡ 474085ee-1a0f-458a-af9a-c746e988dff9
const nao = OffsetArray([make_nao(l) for l in 0:4],0:4);

# ╔═╡ 6a50233d-3948-4a0f-835e-171587e615fc
begin
	# Define variables for the timing functions
	amax = cmax = 4  # Max angular momentum of the Gaussian function
	aexp = bexp = cexp = dexp = 1.0 # exponent of the Gaussian function
	A = [0.0,0.0,0.0] # coordinates of the Gaussian function centers.
	B = [0.0,0.0,1.0]
	C = [0.0,1.0,0.0]
	D = [1.0,0.0,0.0]

	# Allocate main space
    vrrs = zeros(Float64,amax+cmax+1,nao[amax],nao[cmax])
end;

# ╔═╡ 21bf15f1-661e-4859-86d5-12ace608a7b9
const m2ao = make_m2ao();

# ╔═╡ c376bf6b-a153-4efa-9f4a-0e299377a170
function make_shift_index()
    n = length(ao2m)
    shift_index = zeros(Int,n,3)
    for a in 1:n
        m = ao2m[a]
        for i in 1:3
            if m[i] == 0
                shift_index[a,i] = 0
            else
                mm = copy(m)
                mm[i] -= 1
                am = m2ao[mm]
                shift_index[a,i] = am
            end
        end
    end
    return shift_index
end;

# ╔═╡ 21fafb54-d7c7-4dfc-83ad-7142a618ab91
const shift_index = make_shift_index();

# ╔═╡ 52972c5d-2d1c-487f-826e-a9d1196219f7
"vrr!(vrrs, amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns an `n`x`m` array, where `n` is the number
of aos in the `a+b` shell, and `m` is the number of aos in the
`c+d` shell.

Equations refer to
[A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations](https://doi.org/10.1063/1.455553). Martin Head-Gordon and John A. Pople. JCP, 89 (9), 5777, 1988."
function vrr!(vrrs,amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1

	# Prefactors arising from multiplying Gaussians together
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze # "ooz = one over zeta, etc."
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*oortze*exp(
		-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12 using an array of calls to
	# the incomplete gamma function. This is ~40% of the time here,
	# but can be replaced with table lookup or interpolation. The
	# simple version is included here.
    boys_array = boys_array_gamma(mmax,T)
    for m in 1:mmax
        vrrs[m,1,1] = KabKcd_rtze*boys_array[m]
    end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        if aminus > 0
            a_i = 0.5*ooz*ao2m[a][i]
            for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1] +
					a_i*(vrrs[m,aminus,1]-eta*ooze*vrrs[m+1,aminus,1])
            end
        else
            for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1]
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        cminus = shift_index[c,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            aminus = shift_index[a,i]
            if cminus > 0
                c_i = 0.5*ooe*ao2m[c][i]
                if aminus > 0
                    a_i = 0.5*ooze*ao2m[a][i]
                    for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            a_i*vrrs[m+1,aminus,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                else
                    for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                end
            elseif aminus > 0
                a_i = 0.5*ooze*ao2m[a][i]
                for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                        a_i*vrrs[m+1,aminus,c]
                end
            else
                for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c]
                end
            end
        end
    end
    return nothing
end;

# ╔═╡ 7cb6a914-ea03-4a30-956e-67be78c4547c
@btime vrr!($vrrs, $amax,$cmax, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ 53b65f59-8fcd-4124-b1d1-3e5f8fd85417
@btime vrr!($vrrs, 1,1, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ fe06d209-8826-4410-b83e-ccf6e28136bf
@profile (for _ in 1:500; vrr!(vrrs, amax,cmax, aexp,bexp,cexp,dexp, A,B,C,D);end)

# ╔═╡ 7ebda10b-54e4-4def-bb72-831dc9637183
"vrr_turbo!(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns an `n`x`m` array, where `n` is the number
of aos in the `a+b` shell, and `m` is the number of aos in the
`c+d` shell.

Equations refer to
[A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations](https://doi.org/10.1063/1.455553). Martin Head-Gordon and John A. Pople. JCP, 89 (9), 5777, 1988."
function vrr_turbo!(vrrs,amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1

	# Prefactors arising from multiplying Gaussians together
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze # "ooz = one over zeta, etc."
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*oortze*exp(
		-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12 using an array of calls to
	# the incomplete gamma function. This is ~40% of the time here,
	# but can be replaced with table lookup or interpolation. The
	# simple version is included here.
    boys_array = boys_array_gamma(mmax,T)
    @turbo for m in 1:mmax
        vrrs[m,1,1] = KabKcd_rtze*boys_array[m]
    end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        if aminus > 0
            a_i = 0.5*ooz*ao2m[a][i]
            @turbo for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1] +
					a_i*(vrrs[m,aminus,1]-eta*ooze*vrrs[m+1,aminus,1])
            end
        else
            @turbo for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1]
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        cminus = shift_index[c,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            aminus = shift_index[a,i]
            if cminus > 0
                c_i = 0.5*ooe*ao2m[c][i]
                if aminus > 0
                    a_i = 0.5*ooze*ao2m[a][i]
                    @turbo for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            a_i*vrrs[m+1,aminus,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                else
                    @turbo for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                end
            elseif aminus > 0
                a_i = 0.5*ooze*ao2m[a][i]
                @turbo for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                        a_i*vrrs[m+1,aminus,c]
                end
            else
                @turbo for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c]
                end
            end
        end
    end
    return nothing
end;

# ╔═╡ d07a14a5-6a5f-4889-bdd3-89e8b134ffdf
@btime vrr_turbo!($vrrs, $amax,$cmax, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ 79740deb-f4d1-4c55-8a9c-f709ebcbfbe2
@btime vrr_turbo!($vrrs, 1,1, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ 581b59ec-4d43-4601-bd9a-990a50261603
"vrr_inbounds!(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns an `n`x`m` array, where `n` is the number
of aos in the `a+b` shell, and `m` is the number of aos in the
`c+d` shell.

Equations refer to
[A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations](https://doi.org/10.1063/1.455553). Martin Head-Gordon and John A. Pople. JCP, 89 (9), 5777, 1988."
function vrr_inbounds!(vrrs,amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1

	# Prefactors arising from multiplying Gaussians together
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze # "ooz = one over zeta, etc."
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*oortze*exp(
		-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12 using an array of calls to
	# the incomplete gamma function. This is ~40% of the time here,
	# but can be replaced with table lookup or interpolation. The
	# simple version is included here.
    boys_array = boys_array_gamma(mmax,T)
    for m in 1:mmax
        @inbounds vrrs[m,1,1] = KabKcd_rtze*boys_array[m]
    end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        if aminus > 0
            a_i = 0.5*ooz*ao2m[a][i]
            for m in 1:lim
                @inbounds vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1] +
					a_i*(vrrs[m,aminus,1]-eta*ooze*vrrs[m+1,aminus,1])
            end
        else
            for m in 1:lim
                @inbounds vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1]
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        cminus = shift_index[c,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            aminus = shift_index[a,i]
            if cminus > 0
                c_i = 0.5*ooe*ao2m[c][i]
                if aminus > 0
                    a_i = 0.5*ooze*ao2m[a][i]
                    for m in 1:lim
                        @inbounds vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            a_i*vrrs[m+1,aminus,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                else
                    for m in 1:lim
                        @inbounds vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                end
            elseif aminus > 0
                a_i = 0.5*ooze*ao2m[a][i]
                for m in 1:lim
                    @inbounds vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                        a_i*vrrs[m+1,aminus,c]
                end
            else
                for m in 1:lim
                    @inbounds vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c]
                end
            end
        end
    end
    return nothing
end;

# ╔═╡ ea3639bd-89c5-4b8f-b2d0-066c677824af
@btime vrr_inbounds!($vrrs, $amax,$cmax, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ efb4e0c9-f620-4a38-bf60-660c646a4663
@btime vrr_inbounds!($vrrs, 1,1, $aexp,$bexp,$cexp,$dexp, $A,$B,$C,$D)

# ╔═╡ Cell order:
# ╟─7f1163a2-5a4b-4892-a5ed-cac70c41d3b7
# ╠═b7c63d26-0827-491f-908a-02c3c8bf5d0c
# ╠═e109b7eb-725d-499c-b136-9552fff9da38
# ╟─927fecfa-bb26-11eb-3126-6393fad683fd
# ╟─9282d572-29e6-4c48-8981-4280cf0de817
# ╠═52972c5d-2d1c-487f-826e-a9d1196219f7
# ╟─94bcfcc1-5216-4a2d-a136-0c8d3dae675a
# ╠═7ebda10b-54e4-4def-bb72-831dc9637183
# ╠═d4aa62ea-7bd3-4094-ac77-1d50a97d3b6d
# ╠═581b59ec-4d43-4601-bd9a-990a50261603
# ╠═7fe4e843-7543-45ea-8765-ed04dfab157e
# ╠═6a50233d-3948-4a0f-835e-171587e615fc
# ╠═7cb6a914-ea03-4a30-956e-67be78c4547c
# ╠═d07a14a5-6a5f-4889-bdd3-89e8b134ffdf
# ╠═ea3639bd-89c5-4b8f-b2d0-066c677824af
# ╠═53b65f59-8fcd-4124-b1d1-3e5f8fd85417
# ╠═79740deb-f4d1-4c55-8a9c-f709ebcbfbe2
# ╠═efb4e0c9-f620-4a38-bf60-660c646a4663
# ╟─1fb56ae0-654d-4288-9da7-1ac4869115a4
# ╠═fe06d209-8826-4410-b83e-ccf6e28136bf
# ╠═63626a39-616b-4264-932d-114e3e7aca30
# ╟─92bf1847-1b11-4a2d-9231-ef629ccbf0f6
# ╠═d01acace-1209-425f-826a-5efadb27b947
# ╠═9174abad-51b6-4dbe-b326-2ee25b56b8af
# ╠═0b75d079-bdb6-48c2-bb7f-d3d257bc17f6
# ╠═b4dde3a6-7ab4-4f4f-8270-ff32d8d3e12b
# ╠═7111c63d-bbc7-4343-8941-0b87228a9694
# ╠═3b321740-b76a-4b81-896d-a328398e6477
# ╠═f364fe16-17d3-42d0-9ad4-f562168a4a35
# ╠═f6b8a5d1-18fd-41d3-bd2d-9eb85001af2c
# ╠═c376bf6b-a153-4efa-9f4a-0e299377a170
# ╠═ec367e77-8eb0-4e66-89a1-249841e1ecfe
# ╠═82f11fe4-e95d-4f3d-be43-e65b9c0b3058
# ╟─ce86e563-3dee-4c0c-9cc7-a3efbfc3e326
# ╠═828390fb-fe8c-4a8d-8882-642ea78e3555
# ╠═58071e7c-13d3-406d-9e93-4ca5e0906b6e
# ╠═b9142292-e0e8-4bb7-98df-16ab5c23a8c6
# ╠═21fafb54-d7c7-4dfc-83ad-7142a618ab91
# ╠═c6d84970-9d30-4214-a456-eb03db812aa9
# ╠═474085ee-1a0f-458a-af9a-c746e988dff9
# ╠═21bf15f1-661e-4859-86d5-12ace608a7b9


text:
# f(x)=mean(abs(x))

#TODO sum(abs)

function FmeanAbs(x...)
    tmp=0.0
    for i in 1:length(x)
        tmp+=mean(abs(x[i]))
    end
    return ([tmp/length(x)],nothing)
end

function FmeanAbs_inplace(value::Array,auxvalue,x...) # inplace
    tmp=0.0
    for i in 1:length(x)
        tmp+=mean(abs(x[i]))
    end
    fill!(value,tmp/length(x))
end

function DmeanAbs(derivativeIDX,f_c,faux_c,grad_c,grad_n,x...)
    axpy!(grad_c[1]/(length(x)*length(x[derivativeIDX])),sign(x[derivativeIDX]),grad_n)

end

if PROC=="GPU"

    function meanAbs(A::CudaArray)
        return CudaArray(CUBLAS.asum(A)/length(A)*ones(1,1))
    end
    export meanAbs

    function meanAbs!(A::CudaArray,Out::CudaArray)
        copy!(Out,meanAbs(A))
    end
    export meanAbs!

    function FmeanAbs(x::CudaArray...)
        tmp=CudaArray(zeros(1,1))
        for i in 1:length(x)
            axpy!(1.0/length(x),meanAbs(x[i]),tmp)
        end
        return (tmp,nothing)
    end

    function FmeanAbs_inplace(value::CudaArray,auxvalue,x::CudaArray...) # inplace
        fill!(value,0.0)
        for i in 1:length(x)
            axpy!(1.0/length(x),meanAbs(x[i]),value)
        end
    end

    function DmeanAbs(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray...)
        tmp=CudaArray(Float64,size(grad_n))
        vsign!(x[derivativeIDX],tmp)
        alphaaxpy!(1.0/(length(x)*length(x[derivativeIDX])),grad_c,tmp,grad_n)
        free(tmp)
    end

end

Derivative[FmeanAbs]=DmeanAbs # Define dictionary lookup
Inplace[FmeanAbs]=FmeanAbs_inplace

meanAbs(n::ADnode)=ADnode(FmeanAbs,n)

function meanAbs(n::ArrayADnode)
    return ADnode(FmeanAbs,n)
end

#mean(A::ADtrans)=ADnode(Fmean, ftranspose(node[A.parent]))
meanAbs(A::ADtrans)=ADnode(FmeanAbs, node[A.parent]) # mean(A')=mean(A)

export meanAbs

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
struct QueryableMap <: Queryable
    source
    f_func
    f_expr
    getiterator
end

function QueryOperators.map(source::Queryable, f::Function, f_expr::Expr)
    return QueryableMap(source, f, f_expr, source.getiterator)
end


text:
# Autogenerated wrapper script for AsmJit_jll for armv7l-linux-musleabihf
export libasmjit

JLLWrappers.@generate_wrapper_header("AsmJit")
JLLWrappers.@declare_library_product(libasmjit, "libasmjit.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libasmjit,
        "lib/libasmjit.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Learn
using Base.Test

failed = []
@testset "Learn" begin
    for pkg in ["LearnBase",
                "LossFunctions",
                "Transformations",
                "PenaltyFunctions",
                "ObjectiveFunctions",
                "StochasticOptimization",
                "MLMetrics",
                "ValueHistories",
                "MLDataUtils",
                ]
        try
            Pkg.test(pkg)
            @test true
        catch err
            warn(pkg, " had test errors")
            push!(failed, pkg)
            @test false
        end
    end

    if !isempty(failed)
        warn("The following packages failed the tests:\n    $(join(failed, "\n    "))")
    end
end

# @test passed


text:
using Debugger: Debugger, @enter, execute_command, RunDebugger
import JuliaInterpreter: JuliaInterpreter, JuliaStackFrame, @lookup, Compiled, pc_expr, @make_stack, finish!

using Test
using .Meta: isexpr
import REPL

#@testset "Main tests" begin
    include("utils.jl")
    include("misc.jl")
    include("stepping.jl")
    include("ui.jl")
#end