text:
import Test: @test, @testset, @test_throws
import Luna: Polarisation, Output, Modes
import LinearAlgebra: norm

@testset "Jones" begin
    h = Polarisation.H()
    @test norm(Polarisation.Stokes(h) - [1; 1; 0; 0]) < 1e-15

    QWP = Polarisation.rotate(Polarisation.WP(π/2), 0)
    @test norm(Polarisation.Stokes(QWP*h) - [1; 1; 0; 0]) < 1e-15
    QWP = Polarisation.rotate(Polarisation.WP(π/2), π/4)
    @test norm(Polarisation.Stokes(QWP*h) - [1; 0; 0; 1]) < 1e-15
    QWP = Polarisation.rotate(Polarisation.WP(π/2), -π/4)
    @test norm(Polarisation.Stokes(QWP*h) - [1; 0; 0; -1]) < 1e-15

    HWP = Polarisation.rotate(Polarisation.WP(π), 0)
    @test norm(Polarisation.Stokes(HWP*h) - [1; 1; 0; 0]) < 1e-15
    HWP = Polarisation.rotate(Polarisation.WP(π), π/4)
    @test norm(Polarisation.Stokes(HWP*h) - [1; -1; 0; 0]) < 1e-15
    HWP = Polarisation.rotate(Polarisation.WP(π), π/8)
    @test norm(Polarisation.Stokes(HWP*h) - [1; 0; 1; 0]) < 1e-15
    HWP = Polarisation.rotate(Polarisation.WP(π), -π/8)
    @test norm(Polarisation.Stokes(HWP*h) - [1; 0; -1; 0]) < 1e-15

    @test isapprox(Polarisation.ellipticity(Polarisation.Stokes(h)), 0.0)
    QWP = Polarisation.rotate(Polarisation.WP(π/2), π/4)
    @test isapprox(Polarisation.ellipticity(Polarisation.Stokes(QWP*h)), 1.0)
    QWP = Polarisation.rotate(Polarisation.WP(π/2), -π/4)
    @test isapprox(Polarisation.ellipticity(Polarisation.Stokes(QWP*h)), 1.0)

    QWP = Polarisation.rotate(Polarisation.WP(π/2), π/4)
    C = QWP*h
    L = Polarisation.LP()*C
    @test norm(Polarisation.Stokes(L, normalise=true) - [1; 1; 0; 0]) < 1e-15
    @test isapprox(Polarisation.Stokes(L)[1], 0.5)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "string indexing" begin
    let str = "this is a test\xed\x80"
        @test iterate(str, 15) == (reinterpret(Char, 0xed800000), 17)
        @test_throws BoundsError getindex(str, 0:3)
        @test_throws BoundsError getindex(str, 17:18)
        @test_throws BoundsError getindex(str, 2:17)
        @test_throws BoundsError getindex(str, 16:17)
        @test string(Char(0x110000)) == "\xf4\x90\x80\x80"
    end
end

@testset "string reverse" begin
    @test reverse("") == ""
    @test reverse("a") == "a"
    @test reverse("abc") == "cba"
    @test reverse("xyz\uff\u800\uffff\U10ffff") == "\U10ffff\uffff\u800\uffzyx"
    for (s, r) in [
        "xyz\xc1"         => "\xc1zyx",
        "xyz\xd0"         => "\xd0zyx",
        "xyz\xe0"         => "\xe0zyx",
        "xyz\xed\x80"     => "\xed\x80zyx",
        "xyz\xf0"         => "\xf0zyx",
        "xyz\xf0\x80"     => "\xf0\x80zyx",
        "xyz\xf0\x80\x80" => "\xf0\x80\x80zyx",
    ]
        @test reverse(s) == r
    end
end

@testset "string convert" begin
    @test String(b"this is a test\xed\x80\x80") ==
                  "this is a test\xed\x80\x80"  ==
                  "this is a test\ud000"
    # Specifically check UTF-8 string whose lead byte is same as a surrogate
    @test String(b"\xed\x9f\xbf") == "\xed\x9f\xbf" == "\ud7ff"
end

ext: jl

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
"""
KernelFunctions. [Github](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl)
[Documentation](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/)
"""
module KernelFunctions

export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!
export transform
export duplicate, set! # Helpers

export Kernel
export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel, WienerKernel
export CosineKernel
export SqExponentialKernel, RBFKernel, GaussianKernel, SEKernel
export LaplacianKernel, ExponentialKernel, GammaExponentialKernel
export ExponentiatedKernel
export MaternKernel, Matern32Kernel, Matern52Kernel
export LinearKernel, PolynomialKernel
export RationalQuadraticKernel, GammaRationalQuadraticKernel
export MahalanobisKernel, GaborKernel, PiecewisePolynomialKernel
export PeriodicKernel, NeuralNetworkKernel
export KernelSum, KernelProduct
export TransformedKernel, ScaledKernel
export TensorProduct

export Transform, SelectTransform, ChainTransform, ScaleTransform, LinearTransform,
    ARDTransform, IdentityTransform, FunctionTransform

export NystromFact, nystrom

export spectral_mixture_kernel, spectral_mixture_product_kernel

using Compat
using Requires
using Distances, LinearAlgebra
using SpecialFunctions: loggamma, besselk, polygamma
using ZygoteRules: @adjoint, pullback
using StatsFuns: logtwo
using InteractiveUtils: subtypes
using StatsBase

"""
Abstract type defining a slice-wise transformation on an input matrix
"""
abstract type Transform end

abstract type Kernel end
abstract type BaseKernel <: Kernel end
abstract type SimpleKernel <: BaseKernel end

include("utils.jl")
include("distances/pairwise.jl")
include("distances/dotproduct.jl")
include("distances/delta.jl")
include("distances/sinus.jl")
include("transform/transform.jl")

for f in readdir(joinpath(@__DIR__, "basekernels"))
    endswith(f, ".jl") && include(joinpath("basekernels", f))
end

include("kernels/transformedkernel.jl")
include("kernels/scaledkernel.jl")
include("matrix/kernelmatrix.jl")
include("kernels/kernelsum.jl")
include("kernels/kernelproduct.jl")
include("kernels/tensorproduct.jl")
include("approximations/nystrom.jl")
include("generic.jl")

include("zygote_adjoints.jl")

function __init__()
    @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" include("matrix/kernelkroneckermat.jl")
    @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" include("matrix/kernelpdmat.jl")
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("trainable.jl")
end

end


text:
using JuMP, Cbc, CSV, PrettyTables

Raw_Plants_Data = CSV.read("Raw-Plants.csv", header=true)
Plants_DC_Data = CSV.read("Plants-DC.csv", header=true)
DC_Wares_Data = CSV.read("DC-Wares.csv", header=true)
Wares_Stores_Data = CSV.read("Wares-Stores.csv", header=true)

Raw_Plants_Costs = Matrix{Float64}(Raw_Plants_Data[1:4,2:5])
Plants_DC_Costs = Matrix{Float64}(Plants_DC_Data[1:4,2:5])
DC_Wares_Costs = Matrix{Float64}(DC_Wares_Data[1:4,2:5])
Wares_Stores_Costs = Matrix{Float64}(Wares_Stores_Data[1:4,2:5])

Raw_Fixed = Array{Float64}(Raw_Plants_Data[1:4,6])
Plants_Fixed = Array{Float64}(Plants_DC_Data[1:4,6])
DC_Fixed = Array{Float64}(DC_Wares_Data[1:4,6])
Wares_Fixed = Array{Float64}(Wares_Stores_Data[1:4,6])
Stores_Fixed = 9500

Raw_Supply = Array{Float64}(Raw_Plants_Data[1:4,7])
Store_Demand = transpose(Array{Float64}(Wares_Stores_Data[5,2:5]))

numF = 4

caseModel = Model(Cbc.Optimizer)

@variables(caseModel, begin
    RawToPlants[1:numF,1:numF] >= 0, Int
    PlantsToDC[1:numF,1:numF] >= 0, Int
    DCToWares[1:numF,1:numF] >= 0, Int
    WaresToStores[1:numF,1:numF] >= 0, Int
end)

@variables(caseModel, begin
    RawOpen[1:numF], Bin
    PlantsOpen[1:numF], Bin
    DCOpen[1:numF], Bin
    WaresOpen[1:numF], Bin
    StoresOpen[1:numF], Bin
end)

fixedCosts = @expression(caseModel, sum(RawOpen.*Raw_Fixed) + sum(PlantsOpen.*Plants_Fixed) + sum(DCOpen.*DC_Fixed) + sum(WaresOpen.*Wares_Fixed) + sum(StoresOpen.*Stores_Fixed))

variableCosts = @expression(caseModel, sum(RawToPlants.*Raw_Plants_Costs) + sum(PlantsToDC.*Plants_DC_Costs) + sum(DCToWares.*DC_Wares_Costs) + sum(WaresToStores.*Wares_Stores_Costs))

@objective(caseModel, Min, fixedCosts + variableCosts)

for n = 1:numF
    @constraints(caseModel, begin
        sum(WaresToStores[1:numF,n]) == Store_Demand[n]
        sum(RawToPlants[n,1:numF]) <= Raw_Supply[n]*.9*RawOpen[n]
        PlantsToDC[:,n] in SOS1()
    end)

    @constraints(caseModel, begin
        sum(RawToPlants[:,n]) == sum(PlantsToDC[n,:])
        sum(PlantsToDC[:,n]) == sum(DCToWares[n,:])
        sum(DCToWares[:,n]) == sum(WaresToStores[n,:])
    end)

    @constraints(caseModel, begin
        sum(RawToPlants[:,n]) <= PlantsOpen[n]*typemax(Int)
        sum(PlantsToDC[:,n]) <= DCOpen[n]*typemax(Int)
        sum(DCToWares[:,n]) <= WaresOpen[n]*typemax(Int)
        sum(WaresToStores[:,n]) <= StoresOpen[n]*typemax(Int)
    end)
end

JuMP.optimize!(caseModel)

RTP=JuMP.value.(RawToPlants)
PTD=JuMP.value.(PlantsToDC)
DTW=JuMP.value.(DCToWares)
WTS=JuMP.value.(WaresToStores)

println("Epa! no se imprime el primer print")
println("Amounts")
println(JuMP.value.(RawToPlants))
println(JuMP.value.(PlantsToDC))
println(JuMP.value.(DCToWares))
println(JuMP.value.(WaresToStores))
println("Opens")
println(JuMP.value.(RawOpen))
println(JuMP.value.(PlantsOpen))
println(JuMP.value.(DCOpen))
println(JuMP.value.(WaresOpen))
println(JuMP.value.(StoresOpen))

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
# https://tutorials.yaoquantum.org/dev/
using YaoExtensions, IBMQJulia, Yao, Test

@testset "QFT" begin
    for i in [2, 4, 5, 7, 8]
        q = qft_circuit(i)
        qobj = IBMQJulia.yaotoqobj([q], "foo_device")
        exp_1 = qobj.data["experiments"]
        ins = exp_1[1]["instructions"]
        @test qobj isa IBMQJulia.Qobj
        @test qobj.data isa Dict{String,Any}
        @test exp_1 isa Array{Dict{String,Any},1}
        @test ins isa Array{Any, 1}
        @test length(ins) == i*(i+1)/2
        for j in ins
            @test j isa Dict{String, Any}
        end
    end
    c = qft_circuit(4)
    inst = generate_inst(c)
    c2 = inst |> inst2qbir
    @test operator_fidelity(c, c2) ≈ 1
end

@testset "Quantum Circuit Born Machine" begin
    layer(nbit::Int, x::Symbol) = layer(nbit, Val(x))
    layer(nbit::Int, ::Val{:first}) = chain(nbit, put(i=>chain(Rx(0), Rz(0))) for i = 1:nbit)
    layer(nbit::Int, ::Val{:last}) = chain(nbit, put(i=>chain(Rz(0), Rx(0))) for i = 1:nbit)
    layer(nbit::Int, ::Val{:mid}) = chain(nbit, put(i=>chain(Rz(0), Rx(0), Rz(0))) for i = 1:nbit)
    entangler(pairs) = chain(control(ctrl, target=>X) for (ctrl, target) in pairs)

    function build_circuit(n, nlayers, pairs)
        circuit = chain(n)
        push!(circuit, layer(n, :first))
        for i in 2:nlayers
            push!(circuit, cache(entangler(pairs)))
            push!(circuit, layer(n, :mid))
        end
        push!(circuit, cache(entangler(pairs)))
        push!(circuit, layer(n, :last))
        return circuit
    end

    qc = build_circuit(4, 1, [1=>2, 2=>3, 3=>4])
    qobj = IBMQJulia.yaotoqobj([qc], "foo_device")
    exp_1 = qobj.data["experiments"]
    ins = exp_1[1]["instructions"]
    @test qobj isa IBMQJulia.Qobj
    @test qobj.data isa Dict{String,Any}
    @test exp_1 isa Array{Dict{String,Any},1}
    @test ins isa Array{Any,1}
    for j in ins
        @test j isa Dict{String,Any}
    end
    c = build_circuit(4, 1, [1=>2, 2=>3, 3=>4])
    inst = generate_inst(c)
    c2 = inst |> inst2qbir
    @test operator_fidelity(c, c2) ≈ 1
end

@testset "Variational Quantum Eigen Solver" begin
    for n in [3, 6, 7, 8, 9]
        circuit = dispatch!(variational_circuit(n, n+1),:random)
        qobj = IBMQJulia.yaotoqobj([circuit], "foo_device")
        exp_1 = qobj.data["experiments"]
        ins = exp_1[1]["instructions"]
        @test qobj isa IBMQJulia.Qobj
        @test qobj.data isa Dict{String,Any}
        @test exp_1 isa Array{Dict{String,Any},1}
        @test ins isa Array{Any,1}
        for j in ins
            @test j isa Dict{String,Any}
        end
    end
    c = dispatch!(variational_circuit(4, 5), :random)
    inst = generate_inst(c)
    c2 = inst |> inst2qbir
    @test operator_fidelity(c, c2) ≈ 1
end

@testset "Misc" begin
    qc = chain(3, put(1=>X), put(2=>Y) ,put(3=>Z),
                put(2=>T), swap(1,2), put(3=>Ry(0.7)),
                control(2, 1=>Y), control(3, 2=>Z))
    qobj = IBMQJulia.yaotoqobj([qc], "foo_device")
    exp_1 = qobj.data["experiments"]
    ins = exp_1[1]["instructions"]
    @test qobj isa IBMQJulia.Qobj
    @test qobj.data isa Dict{String,Any}
    @test exp_1 isa Array{Dict{String,Any},1}
    @test ins isa Array{Any,1}
    for j in ins
        @test j isa Dict{String,Any}
    end

    c = chain(3, put(1=>X), put(2=>Y) ,put(3=>Z),
                put(2=>T), swap(1,2), put(3=>Ry(0.7)),
                control(2, 1=>Y), control(3, 2=>Z))
    inst = generate_inst(c)
    c2 = inst |> inst2qbir
    @test operator_fidelity(c, c2) ≈ 1
end


text:
config = MOI.Test.TestConfig()
# The test does not check the solution so we just set zeros.
optimize!(mock) = MOIU.mock_optimize!(mock, zeros(MOI.get(mock, MOI.NumberOfVariables())))
for mock in mocks(optimize!)
    Tests.sosdemo9_test(mock, config)
end


text:
module ProofOfConcept

using LinearAlgebra

import Base: *, ==

export Triangle
export Point
export get_barycentric_coordinates
export *
export ==
export is_inside
export refinement_criteria
export TerrainMap
export look_for_points_inside_triangle
export distribute_and_mark_for_refinement

# Write your package code here.
struct Point
    coords::Vector{Real}
    height::Real
end

function Point(coords)
    return Point(coords, NaN)
end

function get_coords(p::Point)
    return p.coords
end

function ==(p::Point, v)
return all(p.coords ≈ v)
end

struct Triangle
    coords::Matrix{Real}
    height::Vector{Real}
end

function Triangle(m)
    t = Triangle(m, similar(m, length(m)))
end

function Triangle(p1, p2, p3)
    return Triangle([p1 p2 p3])
end

function Triangle(p1::Point, p2, p3)
    return Triangle([get_coords(p1) p2 p3])
end

function get_v1(t)
    return t.coords[:,1]
end

function get_v2(t)
    return t.coords[:,2]
end

function get_v3(t)
    return t.coords[:,3]
end

function get_area(t)
    t_v1 = get_v2(t) - get_v1(t)
    t_v2 = get_v3(t) - get_v1(t)

    return det([t_v1 t_v2])/2
end

function get_barycentric_coordinates(p, t)
    t_area = get_area(t)
    t1 = Triangle(p, get_v2(t), get_v3(t))
    t1_area = get_area(t1)
    t2 = Triangle(p, get_v3(t), get_v1(t))
    t2_area = get_area(t2)
    t3 = Triangle(p, get_v1(t), get_v2(t))
    t3_area = get_area(t3)

    return [t1_area/t_area, t2_area/t_area, t3_area/t_area]
end

function ⪅(a, b)
return (a < b) || (a ≈ b)
end

function is_inside(bc)
    return all(0 .⪅ bc .⪅ 1)
end

function is_inside(p::Point, t::Triangle)
    bc = get_barycentric_coordinates(p, t)
    return is_inside(bc)
end

function *(bc, t::Triangle)
return t.coords * bc
end

function refinement_criteria(bc, p::Point, t::Triangle, tol::Real)
    t_height = dot(bc, t.height)
    return abs(p.height - t_height) > tol
end

function refinement_criteria(p::Point, t::Triangle, tol::Real)
    bc = get_barycentric_coordinates(p, t)
    return refinement_criteria(bc, p, t, tol)
end

# The map starts at the point (lon_min, lat_min)
# and moves first on longitude and then on latitude
struct TerrainMap
    lat_min::Real
    lon_min::Real
    Δ_lat::Real
    Δ_lon::Real
    height::Matrix{Real}
end

function get_point(m::TerrainMap, idx::CartesianIndex)
    lon = m.lon_min + (idx[1] - 1) * m.Δ_lon
    lat = m.lat_min + (idx[2] - 1) * m.Δ_lat
    height = m.height[idx]
    return Point([lon, lat], height)
end

function get_point(idx::CartesianIndex, m::TerrainMap)
    return get_point(m, idx)
end

function is_inside(m::TerrainMap, idx::CartesianIndex, t::Triangle)
    return is_inside(get_point(m, idx), t)
end

function refinement_criteria(m::TerrainMap, idx::CartesianIndex, t::Triangle, tol::Real)
    return refinement_criteria(get_point(m, idx), t, tol)
end

function look_for_points_inside_triangle(m::TerrainMap, t::Triangle)
    return idx -> is_inside(m, idx, t)
end

function get_barycentric_coordinates(m::TerrainMap, idx::CartesianIndex, t::Triangle)
    return get_barycentric_coordinates(get_point(m, idx), t)
end

function distribute_and_mark_for_refinement(m::TerrainMap, parent_indices, son_t::Triangle, tol::Real)
    pts = map(x -> get_point(m, x), parent_indices)
    bc = map(x -> get_barycentric_coordinates(x, son_t), pts)
    local_indices = findall(is_inside, bc)
    son_indices = view(parent_indices, local_indices)
    refine = any(i -> refinement_criteria(bc[i], pts[i], son_t, tol), local_indices)
    return (refine, son_indices)
end

end


text:
@testset "solve_advection.jl" begin
    domain = PeriodicIntervalDomain(0.0, 2π)

    T = 1.0

    # Exact solutions
    u(ξ, t) = sin(ξ - t) + 3 / 5 * cos(5(ξ - t)) + 1 / 25 * sin(20(ξ - 1 - t))
    u_int(ξ, t) = -cos(ξ - t) + 3 / 25 * sin(5(ξ - t)) - 1 / 25 / 20 * cos(20(ξ - 1 - t))

    params = (; subspacedim = 50, reltol = 1e-7, abstol = 1e-9)

    # Number of mesh points
    NN = floor.(Int, 10 .^ LinRange(2.5, 4.0, 4))

    err_ref = 0
    err_bar_ref = 0
    err_allbar_ref = 0
    n_ref = 0
    for (i, N) ∈ enumerate(NN)
        M = N

        # Discretization
        x = discretize(domain, M)
        ξ = discretize(domain, N)
        Δx = (domain.right - domain.left) / N

        # Filter
        h(x) = Δx / 2

        # Exact filtered solution
        ū(x, t) = 1 / 2h(x) * (u_int(x + h(x), t) - u_int(x - h(x), t))

        # Equations
        equation = AdvectionEquation(domain, IdentityFilter())
        equation_filtered = AdvectionEquation(domain, TopHatFilter(h))

        # Solve discretized problem
        sol = solve(
            equation,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "discretizefirst",
            params...,
        )

        # Solve filtered-then-discretized problem
        sol_bar = solve(
            equation_filtered,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "filterfirst",
            params...,
        )

        # Solve discretized-then-filtered problem
        sol_allbar = solve(
            equation_filtered,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "discretizefirst",
            params...,
        )

        ## Relative error
        u_exact = u.(ξ, T)
        ū_exact = ū.(x, T)
        err = norm(sol(T) - u_exact) / norm(u_exact)
        err_bar = norm(sol_bar(T) - ū_exact) / norm(ū_exact)
        err_allbar = norm(sol_allbar(T) - ū_exact) / norm(ū_exact)

        if i == 1
            err_ref = err
            err_bar_ref = err_bar
            err_allbar_ref = err_allbar
            n_ref = N
        end

        # Test for quadratic convergence
        @test err < 2 * err_ref * (n_ref / N)^2
        @test err_bar < 2 * err_bar_ref * (n_ref / N)^2
        @test err_allbar < 2 * err_allbar_ref * (n_ref / N)^2
    end
end


text:
"""
`Quaternion` represents ``q_w + q_x \\bm{i} + q_y \\bm{j} + q_z \\bm{k}``.
The salar part and vector part can be accessed by `q.scalar` and `q.vector`, respectively.

# Examples
```jldoctest
julia> Quaternion(1,2,3,4)
1 + 2𝙞 + 3𝙟 + 4𝙠

julia> Quaternion(1)
1 + 0𝙞 + 0𝙟 + 0𝙠

julia> Quaternion(Vec(1,2,3))
0 + 1𝙞 + 2𝙟 + 3𝙠
```

See also [`quaternion`](@ref).

!!! note

    `Quaternion` is experimental and could change or disappear in future versions of Tensorial.
"""
struct Quaternion{T} <: Number
    data::NTuple{4, T}
    function Quaternion{T}(data::NTuple{4, Real}) where {T}
        new{T}(convert_ntuple(T, data))
    end
end

@inline Quaternion(data::NTuple{4, Any}) = Quaternion{promote_ntuple_eltype(data)}(data)
@inline (::Type{T})(data::Vararg{Any}) where {T <: Quaternion} = T(data)

# from scalar and vector
@inline Quaternion{T}(r::Real, v::Vec{3}) where {T} = @inbounds Quaternion{T}(r, v[1], v[2], v[3])
@inline Quaternion{T}(r::Real, v::Vec{2}) where {T} = @inbounds Quaternion{T}(r, v[1], v[2], zero(eltype(v)))
@inline Quaternion{T}(r::Real, v::Vec{1}) where {T} = @inbounds Quaternion{T}(r, v[1], zero(eltype(v)), zero(eltype(v)))
@inline Quaternion(r::Real, v::Vec) = Quaternion{promote_type(typeof(r), eltype(v))}(r, v)

# from vector
@inline Quaternion{T}(v::Vec{4}) where {T} = Quaternion{T}(Tuple(v))
for dim in 1:3
    @eval @inline Quaternion{T}(v::Vec{$dim}) where {T} = Quaternion{T}(zero(eltype(v)), v)
end
@inline Quaternion(v::Vec) = Quaternion{eltype(v)}(v)

# from scalar
@inline Quaternion{T}(r::Real) where {T} = (z = zero(r); Quaternion{T}(r, z, z, z))
@inline Quaternion(r::Real) = Quaternion{typeof(r)}(r)

@inline Vec(q::Quaternion) = Vec(Tuple(q))

Base.Tuple(q::Quaternion) = getfield(q, :data)

@inline function Base.getproperty(q::Quaternion, name::Symbol)
    name == :scalar && return @inbounds q[1]
    name == :vector && return @inbounds Vec(q[2], q[3], q[4])
    getfield(q, name)
end

Base.propertynames(q::Quaternion) = (:scalar, :vector, :data)

# conversion
Base.convert(::Type{Quaternion{T}}, x::Quaternion{T}) where {T} = x
Base.convert(::Type{Quaternion{T}}, x::Quaternion{U}) where {T, U} = Quaternion(map(T, Tuple(x)))
Base.convert(::Type{Quaternion{T}}, x::Real) where {T} = convert(Quaternion{T}, Quaternion(x))

# promotion
Base.promote_rule(::Type{Quaternion{T}}, ::Type{T}) where {T <: Real} = Quaternion{T}
Base.promote_rule(::Type{Quaternion{T}}, ::Type{U}) where {T <: Real, U <: Real} = Quaternion{promote_type(T, U)}
Base.promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{T}}) where {T <: Real} = Quaternion{T}
Base.promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{U}}) where {T <: Real, U <: Real} = Quaternion{promote_type(T, U)}

# used for `isapprox`
Base.real(q::Quaternion) = q.scalar
Base.isfinite(q::Quaternion) = prod(map(isfinite, Tuple(q)))

"""
    quaternion(θ, n::Vec; [normalize = true, degree = false])

Construct `Quaternion` from angle `θ` and axis `n` as

```math
q = \\cos\\frac{\\theta}{2} + \\bm{n} \\sin\\frac{\\theta}{2}
```

The constructed quaternion is normalized such as `norm(q) ≈ 1` by default.

# Examples
```jldoctest
julia> q = quaternion(π/4, Vec(0,0,1))
0.9238795325112867 + 0.0𝙞 + 0.0𝙟 + 0.3826834323650898𝙠

julia> x = rand(Vec{3})
3-element Vec{3, Float64}:
 0.5908446386657102
 0.7667970365022592
 0.5662374165061859

julia> (q * x / q).vector ≈ rotmatz(π/4) ⋅ x
true
```
"""
function quaternion(::Type{T}, θ::Real, x::Vec{3}; normalize::Bool = true, degree::Bool = false) where {T}
    if degree
        θ = deg2rad(θ)
    end
    ϕ = θ / 2
    if normalize
        n = LinearAlgebra.normalize(x) * sin(ϕ)
    else
        n = x * sin(ϕ)
    end
    @inbounds Quaternion{T}(cos(ϕ), n)
end
quaternion(T::Type, θ::Real, x::Vec{2}; normalize::Bool = true, degree::Bool = false) =
    @inbounds quaternion(T, θ, Vec(x[1], x[2], 0); normalize, degree)
quaternion(θ::Real, x::Vec; normalize::Bool = true, degree::Bool = false) =
    quaternion(promote_type(typeof(θ), eltype(x)), θ, x; normalize, degree)

Base.length(::Quaternion) = 4
Base.size(::Quaternion) = (4,)

@inline function Base.getindex(q::Quaternion, i::Int)
    @boundscheck 1 ≤ i ≤ 4 || throw(BoundsError(q, i))
    @inbounds Tuple(q)[i]
end

# quaternion vs quaternion
@inline Base.:-(q::Quaternion) = Quaternion(-Vec(q))
@inline Base.:+(q::Quaternion, p::Quaternion) = Quaternion(Vec(q) + Vec(p))
@inline Base.:-(q::Quaternion, p::Quaternion) = Quaternion(Vec(q) - Vec(p))
@inline Base.:/(q::Quaternion, p::Quaternion) = q * inv(p)
@inline function Base.:*(q::Quaternion, p::Quaternion)
    q₁, q₂, q₃, q₄ = Tuple(q)
    A = @Mat [ q₁ -q₂ -q₃ -q₄
               q₂  q₁ -q₄  q₃
               q₃  q₄  q₁ -q₂
               q₄ -q₃  q₂  q₁ ]
    Quaternion(A ⋅ Vec(p))
end

# quaternion vs number
@inline Base.:*(a::Number, q::Quaternion) = Quaternion(a * Vec(q))
@inline Base.:*(q::Quaternion, a::Number) = Quaternion(Vec(q) * a)
@inline Base.:/(q::Quaternion, a::Number) = Quaternion(Vec(q) / a)

# quaternion vs vector
@inline Base.:*(q::Quaternion, v::Vec) = q * Quaternion(v)
@inline Base.:*(v::Vec, q::Quaternion) = Quaternion(v) * q
@inline Base.:/(v::Vec, q::Quaternion) = v * inv(q)

function angleaxis(q::Quaternion)
    a = norm(q.vector)
    θ = 2atan(a, q.scalar)
    n = q.vector / a
    θ, n
end

"""
    rotate(x::Vec, q::Quaternion)

Rotate `x` by quaternion `q`.

# Examples
```jldoctest
julia> v = Vec(1.0, 0.0, 0.0)
3-element Vec{3, Float64}:
 1.0
 0.0
 0.0

julia> rotate(v, quaternion(π/4, Vec(0,0,1)))
3-element Vec{3, Float64}:
 0.7071067811865475
 0.7071067811865476
 0.0
```
"""
@inline rotate(v::Vec{3}, q::Quaternion) = (q * v / q).vector
@inline rotate(v::Vec{2}, q::Quaternion) = (v = (q * v / q).vector; @inbounds Vec(v[1], v[2]))

@inline Base.conj(q::Quaternion) = Quaternion(q.scalar, -q.vector)
@inline Base.abs2(q::Quaternion) = (v = Vec(q); dot(v, v))
@inline Base.abs(q::Quaternion) = sqrt(abs2(q))
@inline norm(q::Quaternion) = abs(q)
@inline inv(q::Quaternion) = conj(q) / abs2(q)

"""
    exp(::Quaternion)

Compute the exponential of quaternion as

```math
\\exp(q) = e^{q_w} \\left( \\cos\\| \\bm{v} \\| + \\frac{\\bm{v}}{\\| \\bm{v} \\|} \\sin\\| \\bm{v} \\| \\right)
```
"""
function Base.exp(q::Quaternion)
    v = q.vector
    v_norm = norm(v)
    if v_norm > 0
        n = v / v_norm
    else
        n = zero(v)
    end
    exp(q.scalar) * quaternion(2*v_norm, n; normalize = false)
end

"""
    log(::Quaternion)

Compute the logarithm of quaternion as

```math
\\ln(q) = \\ln\\| q \\| + \\frac{\\bm{v}}{\\| \\bm{v} \\|} \\arccos\\frac{q_w}{\\| q \\|}
```
"""
function Base.log(q::Quaternion)
    q_norm = norm(q)
    Quaternion(log(q_norm), normalize(q.vector) * acos(q.scalar/q_norm))
end

@inline normalize(q::Quaternion) = q / norm(q)

function rotmat_normalized(q::Quaternion)
    q₁, q₂, q₃, q₄ = Tuple(q)
    q₁² = q₁ * q₁
    q₂² = q₂ * q₂
    q₃² = q₃ * q₃
    q₄² = q₄ * q₄
    q₁q₂ = q₁ * q₂
    q₂q₃ = q₂ * q₃
    q₃q₄ = q₃ * q₄
    q₁q₃ = q₁ * q₃
    q₁q₄ = q₁ * q₄
    q₂q₄ = q₂ * q₄
    @Mat [q₁²+q₂²-q₃²-q₄² 2(q₂q₃-q₁q₄)    2(q₂q₄+q₁q₃)
          2(q₂q₃+q₁q₄)    q₁²-q₂²+q₃²-q₄² 2(q₃q₄-q₁q₂)
          2(q₂q₄-q₁q₃)    2(q₃q₄+q₁q₂)    q₁²-q₂²-q₃²+q₄²]
end

"""
    rotmat(::Quaternion)

Construct rotation matrix from quaternion.

# Examples
```jldoctest
julia> q = quaternion(π/4, Vec(0,0,1))
0.9238795325112867 + 0.0𝙞 + 0.0𝙟 + 0.3826834323650898𝙠

julia> rotmat(q)
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.707107  -0.707107  0.0
 0.707107   0.707107  0.0
 0.0        0.0       1.0
```
"""
@inline rotmat(q::Quaternion) = rotmat_normalized(normalize(q))

function Base.show(io::IO, q::Quaternion)
    pm(x) = x < 0 ? " - $(-x)" : " + $x"
    print(io, q[1], pm(q[2]), "𝙞", pm(q[3]), "𝙟", pm(q[4]), "𝙠")
end


text:
using Test, MyAwesomePackage

out = plusTwo(3)
greeting = MyAwesomePackage.greet()

@test out == 5

@test greeting == nothing


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, Langevin, LinearAlgebra, Plots
using ProgressMeter: @showprogress

function Filter(u::Array{Float64,2},lambdaX::Float64,lambdaXi::Float64,filterType::String,PhiL1=0.0,filterfunction=0.0)
    if filterType == "L2"
        for i = 1:size(u,1)
            for j = 1:size(u,2)
                u[i,j] = u[i,j]/(1.0+lambdaX*i^2*(i-1)^2+lambdaXi*j^2*(j-1)^2);
            end
        end
    elseif filterType == "Lasso"
        #PhiL1 = ones(size(u))
        N1 = size(u,1); # number moments in x
        N2 = size(u,2); # number moments in xi
        lambda1 = abs(u[N1,1])/(N1*(N1-1)*PhiL1[N1,1]);
        lambda2 = abs(u[1,N2])/(N2*(N2-1)*PhiL1[1,N2]);
        for i = 1:N1
            for j = 1:N2
                scL1 = 1.0-(lambda1*i*(i-1)+lambda2*j*(j-1))*PhiL1[i,j]/abs(u[i,j]);
                if scL1 < 0 || abs(u[i,j]) < 1.0e-7
                    scL1 = 0;
                end
                u[i,j] = scL1*u[i,j];
            end
        end
    elseif filterType == "HouLi"
        N1 = size(u,1); # number moments in x
        N2 = size(u,2); # number moments in xi
        for i = 1:N1
            for j = 1:N2
                u[i,j] = filterfunction[i,j]*u[i,j];
            end
        end
    end
    return u;
end

begin
    x0 = 0
    x1 = 1
    ncell = 100
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 3 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.05
    dt = cfl * dx
    t = 0.0
end

ps = FRPSpace1D(x0, x1, ncell, deg)

begin
    uqMethod = "galerkin"
    nr = 9
    nRec = 18
    opType = "uniform"
    parameter1 = 0.95
    parameter2 = 1.05
end

uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)

function dudt!(du, u, p, t)
    du .= 0.0
    J, ll, lr, lpdm, dgl, dgr, γ, uq = p

    nm = uq.nm
    nq = uq.nq

    ncell = size(u, 1)
    nsp = size(u, 2)

    u_ran = zeros(ncell, nsp, 3, nq)
    for i = 1:ncell, j = 1:nsp, k = 1:3
        u_ran[i, j, k, :] .= chaos_ran(u[i, j, k, :], uq)
    end

    f = zeros(ncell, nsp, 3, nm+1)
    for i = 1:ncell, j = 1:nsp
        _f = zeros(3, nq)
        for k = 1:nq
            _f[:, k] .= euler_flux(u_ran[i, j, :, k], γ)[1] ./ J[i]
        end

        for k = 1:3
            f[i, j, k, :] .= ran_chaos(_f[k, :], uq)
        end
    end

    f_face = zeros(ncell, 3, nm+1, 2)
    for i = 1:ncell, j = 1:3, k = 1:nm+1
        f_face[i, j, k, 1] = dot(f[i, :, j, k], lr)
        f_face[i, j, k, 2] = dot(f[i, :, j, k], ll)
    end

    u_face = zeros(ncell, 3, nq, 2)
    for i = 1:ncell, j = 1:3, k = 1:nq
        u_face[i, j, k, 1] = dot(u_ran[i, :, j, k], lr)
        u_face[i, j, k, 2] = dot(u_ran[i, :, j, k], ll)
    end

    fq_interaction = zeros(ncell + 1, 3, nq)
    for i = 2:ncell, j = 1:nq
        fw = @view fq_interaction[i, :, j]
        flux_hll!(fw, u_face[i-1, :, j, 1], u_face[i, :, j, 2], γ, 1.0)
    end

    f_interaction = zeros(ncell + 1, 3, nm+1)
    for i = 2:ncell, j = 1:3
        f_interaction[i, j, :] .= ran_chaos(fq_interaction[i, j, :], uq)
    end

    rhs1 = zero(u)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        rhs1[i, ppp1, k, l] = dot(f[i, :, k, l], lpdm[ppp1, :])
    end

    idx = 2:ncell-1 # ending points are Dirichlet
    for i in idx, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        du[i, ppp1, k, l] =
            -(
                rhs1[i, ppp1, k, l] +
                (f_interaction[i, k, l] / J[i] - f_face[i, k, l, 2]) * dgl[ppp1] +
                (f_interaction[i+1, k, l] / J[i] - f_face[i, k, l, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 0.15)
p = (ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, γ, uq)

function HouLiFilterFilterFunction( eta,gammaHL )
    if eta <= 2.0 / 3.0
        return 0.0;
    else
        return eta^gammaHL;
    end
end

begin
    V = vandermonde_matrix(ps.deg,ps.xpl)
    VInv = inv(V)

    Nx = ncell
    nLocal = deg+1
    nStates = 3
    nMoments = uq.nm+1 # number of moments

    lambdaX = 0.01#1e-2
    lambdaXi = 0.01#1e-5
    filterType = "HouLi"#"Lasso"

    # compute L1 norms of basis
    NxHat = 100
    xHat = collect(range(-1,stop = 1,length = NxHat))
    dxHat = xHat[2]-xHat[1];
    PhiL1 = zeros(nLocal,nMoments)
    for i = 1:nLocal
        for j = 1:nMoments
            PhiJ = sum(@. abs(uq.op.quad.weights * uq.phiRan[:, j])) / (uq.t2Product[j-1, j-1] + 1.e-7)
            PhiI = dxHat*sum(abs.(JacobiP(xHat, 0, 0, i-1)))
            PhiL1[i,j] = PhiI*PhiJ
        end
    end

    filterfunction = ones(nLocal,nMoments)
    gammaHL = 1.0#36;
    epsilonXi   = 1.0 / lambdaXi;
    epsilonX   = 1.0 / lambdaX;
    for i = 1:nLocal
        for j = 1:nMoments
            filterfunction[i,j] = exp( -HouLiFilterFilterFunction( j-1 / ( nMoments + 1 ), gammaHL ) / epsilonXi - HouLiFilterFilterFunction( i-1 / ( nLocal + 1 ), gammaHL ) / epsilonX )^(dt*(lambdaXi+lambdaX))
            #filterfunction[i,j] = exp( -HouLiFilterFilterFunction( j-1 / ( nMoments + 1 ), gammaHL ) / epsilonXi) #* exp(-HouLiFilterFilterFunction( i-1 / ( nLocal + 1 ), gammaHL ) / epsilonX )#^36#(dt*(lambdaXi+lambdaX))
        end
    end
end

HouLiFilterFilterFunction( (10-1) / ( nMoments ), gammaHL )

u = zeros(ncell, nsp, 3, uq.nm+1)
# stochastic density
#=for i = 1:ncell, j = 1:nsp
    prim = zeros(3, uq.nm+1)
    if ps.x[i] <= 0.5
        prim[1, :] .= uq.pce
        prim[2, 1] = 0.0
        prim[3, 1] = 0.5
    else
        prim[:, 1] .= [0.3, 0.0, 0.625]
    end

    u[i, j, :, :] .= uq_prim_conserve(prim, γ, uq)
end=#

# stochastic location
for i = 1:ncell, j = 1:nsp
    prim = zeros(3, uq.nq)

    for k = 1:uq.nq
        if ps.x[i] <= 0.5 + 0.05 * uq.op.quad.nodes[k]
            prim[:, k] .= [1.0, 0.0, 0.5]
        else
            prim[:, k] .= [0.4, 0.0, 0.625]
        end
    end

    prim_chaos = zeros(3, uq.nm+1)
    for k = 1:3
        prim_chaos[k, :] .= ran_chaos(prim[k, :], uq)
    end

    u[i, j, :, :] .= uq_prim_conserve(prim_chaos, γ, uq)
end

u1 = deepcopy(u)
#=for j = 1:size(u,1)
    for s = 1:size(u,3)
        uModal = VInv*u[j,:,s,:]
        u1[j,:,s,:] .= V*Filter(uModal,lambdaX,lambdaXi,filterType,PhiL1)
    end
end=#

prob = ODEProblem(dudt!, u1, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:50#nt
    step!(itg)

    # filter
    for j = 1:size(itg.u,1)
        for s = 1:size(itg.u,3)
            uModal = VInv*itg.u[j,:,s,:]
            itg.u[j,:,s,:] .= V*Filter(uModal,lambdaX,lambdaXi,filterType,PhiL1,filterfunction)
        end
    end
end

sol = zeros(ncell, nsp, 3, 2)
for i in axes(sol, 1), j in axes(sol, 2)
    p1 = zeros(3, uq.nm+1)
    p1 = uq_conserve_prim(itg.u[i, j, :, :], γ, uq)
    p1[end, :] .= lambda_tchaos(p1[end, :], 1.0, uq)

    for k = 1:3
        sol[i, j, k, 1] = mean(p1[k, :], uq.op)
        sol[i, j, k, 2] = std(p1[k, :], uq.op)
    end
end

pic1 = plot(ps.x, sol[:, 2, 3, 1], label="mean", xlabel="x", ylabel="ρ")
pic2 = plot(ps.x, sol[:, 2, 1, 2], label="std")
plot(pic1, pic2)

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
# Expected Output According to Reference python : an array of 4841.08061686
lat_ref = 0.7098
lon_ref = 1.2390

lat_test = ones(10) * 0.069
lon_test = ones(10) * 0.069

dlat = sin.((lat_test .- lat_ref) / 2) .^ 2
dlon = sin.((lon_test .- lon_ref) / 2) .^ 2

a = cos.(lat_ref) * cos.(lat_test) .* dlon + dlat

c = asin.(sqrt.(a)) * 2 * 3959.0 # 3959.0 is miles conversion I think
print(c)

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using MultiscaleGraphSignalTransforms, JLD, Plots, LightGraphs, MultivariateStats
using LaTeXStrings
push!(LOAD_PATH, @__DIR__)
using pSGWT

G = loadgraph("../datasets/simple_tree_graph.lgz")
X = load("../datasets/simple_tree_xy.jld", "xy")
N = nv(G)
L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L)
standardize_eigenvectors!(𝚽)
W = 1.0 * adjacency_matrix(G)
Q = incidence_matrix(G; oriented = true)

ib1 = 36:56
ib2 = 21:35
ib3 = 71:100
ib4 = 57:70
ijc = [3,5,12,16]
ir = setdiff(1:N, ib1, ib2, ib3, ib4, ijc)

##
P = 𝚽.^2
# P = exp.(𝚽) ./ sum(exp.(𝚽), dims = 1)

dist_sROT, Ws, Xs, 𝚯 = eigsROT_Distance(P, W, X; α = 1.0)

##
function find_active_eigenvectors(P, interest_locs; threshold = 0.5)
    N = size(P, 1)
    energy = zeros(N)
    for k=1:N
        energy[k] = norm(P[interest_locs, k], 1) / norm(P[:, k], 1)
    end
    ind = findall(energy .> threshold)
    return ind
end

# index of eigenvectors active at branch k (k = 1,2,3,4)
ieb1 = find_active_eigenvectors(𝚽.^2, ib1)
ieb2 = find_active_eigenvectors(𝚽.^2, ib2)
ieb3 = find_active_eigenvectors(𝚽.^2, ib3)
ieb4 = find_active_eigenvectors(𝚽.^2, ib4)
iejc = find_active_eigenvectors(𝚽.^2, ijc; threshold = 0.1)

function simpletree_mds_plot(E, ieb1, ieb2, ieb3, ieb4, iejc)
    scatter_gplot(E'; c = :grey, ms = 2)
    scatter_gplot!(E[:, ieb1]'; c = :pink, ms = 2)
    scatter_gplot!(E[:, ieb2]'; c = :orange, ms = 2)
    scatter_gplot!(E[:, ieb3]'; c = :green, ms = 2)
    scatter_gplot!(E[:, ieb4]'; c = :yellow, ms = 2)
    scatter_gplot!(E[:, iejc]'; c = :red, ms = 2)
    scatter_gplot!(E[:, 1:1]'; c = :magenta, ms = 4)
    plt = plot!(xaxis = "X₁", yaxis = "X₂", zaxis = "X₃", legend = false,
                cbar = false, grid = true)
    return plt
end

##
f = zeros(N)
ib1 = 36:56; ib2 = 21:35; ib3 = 71:100; ib4 = 57:70; ijc = [3,5,12,16]
f[ib1] = sin.(0.3 * (0:(length(ib1) - 1)))
f[ib2] = cos.(0.4 * (0:(length(ib2) - 1)))
f[ib3] = sin.(0.5 * (0:(length(ib3) - 1)))
f[ib4] = cos.(0.6 * (0:(length(ib4) - 1)))
f[ijc] .= 1;


text:
#
# Cubic fit
#

struct Cubic
  a :: Float64
  b :: Float64
  c :: Float64
  d :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

function fitcubic(X :: AbstractVector, Y :: AbstractVector, options :: Options)
  X, Y = checkdata(X,Y)
  @. model(x,p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
  fit = find_best_fit(model, X, Y, 4, options)
  R = pearson(X,Y,model,fit)
  x, y, ypred = finexy(X,options.fine,model,fit)
  return Cubic(fit.param...,R,x,y,ypred,fit.resid)
end
fitcubic(X :: AbstractVector, Y :: AbstractVector) = fitcubic(X,Y,Options())

function Base.show( io :: IO, fit :: Cubic )
  println("")
  println(" ------------------- Cubic Fit ----------------- ")
  println("")
  println(" Equation: y = ax^3 + bx^2 + cx + d ")
  println("")
  println(" With: a = ", fit.a)
  println("       b = ", fit.b)
  println("       c = ", fit.c)
  println("       d = ", fit.d)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitcubic


text:
export BatchExplorer

using Flux

"""
    BatchExplorer(explorer::AbstractExplorer)
"""
struct BatchExplorer{E} <: AbstractExplorer
    explorer::E
end

BatchExplorer(explorers::AbstractExplorer...) = BatchExplorer(explorers)

"""
    (x::BatchExplorer)(values::AbstractMatrix)

Apply inner explorer to each column of `values`.
"""
(x::BatchExplorer)(values::AbstractMatrix) = [x.explorer(v) for v in eachcol(values)]

(x::BatchExplorer{<:Tuple})(values::AbstractMatrix) =
    [explorer(v) for (explorer, v) in zip(x.explorer, eachcol(values))]

Flux.testmode!(x::BatchExplorer, mode = true) = testmode!(x.explorer, mode)

function Flux.testmode!(x::BatchExplorer{<:Tuple}, mode = true)
    for p in x.explorer
        testmode!(p, mode)
    end
end


text:
@testset "Field broadcasting" begin
    @info "  Testing broadcasting with fields..."

    for arch in archs

        #####
        ##### Basic functionality tests
        #####

        grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
        a, b, c = [CenterField(grid) for i = 1:3]

        Nx, Ny, Nz = size(a)

        a .= 1
        @test all(a .== 1)

        b .= 2

        c .= a .+ b
        @test all(c .== 3)

        c .= a .+ b .+ 1
        @test all(c .== 4)

        # Halo regions
        fill_halo_regions!(c, arch) # Does not happen by default in broadcasting now
        @test c[1, 1, 0] == 4
        @test c[1, 1, Nz+1] == 4

        #####
        ##### Broadcasting with interpolation
        #####

        three_point_grid = RectilinearGrid(arch, size=(1, 1, 3), extent=(1, 1, 1))

        a2 = CenterField(three_point_grid)

        b2_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), top=OpenBoundaryCondition(0), bottom=OpenBoundaryCondition(0))
        b2 = ZFaceField(three_point_grid, boundary_conditions=b2_bcs)

        b2 .= 1
        fill_halo_regions!(b2, arch) # sets b2[1, 1, 1] = b[1, 1, 4] = 0

        @test b2[1, 1, 1] == 0
        @test b2[1, 1, 2] == 1
        @test b2[1, 1, 3] == 1
        @test b2[1, 1, 4] == 0

        a2 .= b2
        @test a2[1, 1, 1] == 0.5
        @test a2[1, 1, 2] == 1.0
        @test a2[1, 1, 3] == 0.5

        a2 .= b2 .+ 1
        @test a2[1, 1, 1] == 1.5
        @test a2[1, 1, 2] == 2.0
        @test a2[1, 1, 3] == 1.5

        #####
        ##### Broadcasting with ReducedField
        #####

        r, p, q = [Field{Center, Center, Nothing}(grid) for i = 1:3]

        r .= 2
        @test all(r .== 2)

        p .= 3

        q .= r .* p
        @test all(q .== 6)

        q .= r .* p .+ 1
        @test all(q .== 7)

        #####
        ##### Broadcasting with arrays
        #####

        two_two_two_grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(1, 1, 1))

        c = CenterField(two_two_two_grid)
        random_column = arch_array(arch, reshape(rand(2), 1, 1, 2))

        c .= random_column # broadcast to every horizontal column in c

        c_cpu = Array(interior(c))
        random_column_cpu = Array(random_column)

        @test all(c_cpu[1, 1, :] .== random_column_cpu[:])
        @test all(c_cpu[2, 1, :] .== random_column_cpu[:])
        @test all(c_cpu[1, 2, :] .== random_column_cpu[:])
        @test all(c_cpu[2, 2, :] .== random_column_cpu[:])
    end
end


text:
module FluxGAN

include("models.jl")
using .Models

export GANHyperParams
export GAN

include("image_utilities.jl")
using .ImageUtilities

export save_image_grid
export image_grid_figure
export image_grid_tensor
export image_grid
export color_image
export color_image_tensor

include("train.jl")
using .Train

export train!

end

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
export norm_relation

function _add_sunits_from_norm_relation!(c, UZK, N)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  for i = 1:length(N)
    k, mk = subfield(N, i)
    @vprint :NormRelation 1 "Computing maximal order ..."
    zk = maximal_order(k)
    @vprint :NormRelation 1 "Computing lll basis ... "
    zk = lll(zk)
    @vprint :NormRelation 1 "Computing class group of $k... "
    class_group(zk)
    @vprint :NormRelation 1 "done"
    lpk = NfOrdIdl[ P[1] for p in cp for P = prime_decomposition(zk, p)]
    @vprint :NormRelation 1 "Now computing the S-unit group for lp of length $(length(lpk))"
    @assert length(lpk) > 0
    Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)

    D = Dict{nf_elem, nf_elem}()
    function N_mk(x, D, i)
      if haskey(D, x)
        return D[x]
      else
        y = N(x, i)
        D[x] = y
        return y
      end
    end

    #D = Dict{nf_elem, nf_elem}()
    cc = Vector{Tuple{Int, fmpz}}[]
    for j in 1:length(N.coefficients_gen[i])
      @vprint :NormRelation 1 "Inducing the action ... "
      z = induce_action(N, i, j, lpk, c.FB, cc)

      @vprint :NormRelation "Feeding in the S-units of the small field ... "

      for l=1:ngens(Szk)
        u = mS(Szk[l])  #do compact rep here???
        valofnewelement = mul(mS.valuations[l], z)
        el_to_add = FacElem(Dict{nf_elem, AnticNumberField}((N(x, i, j), v) for (x,v) = u.fac))
        Hecke.class_group_add_relation(c, el_to_add, valofnewelement)#, always = false)
      end
    end

    # Skipping the units
    U, mU = unit_group_fac_elem(zk)
    for j=2:ngens(U) # I cannot add a torsion unit. He will hang forever.
      u = mU(U[j])
      Hecke.add_unit!(UZK, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)))
    end
    UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end
end

function _compute_sunit_and_unit_group!(c, U, N, saturate = true)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  skipped_units = FacElem{nf_elem, AnticNumberField}[]
  autos = automorphisms(field(N), copy = false)

  for i = 1:length(N)
    k, mk = subfield(N, i)
    zk = lll(maximal_order(k))
    print("Computing class group of $k... ")
    @time class_group(zk, redo = true, use_aut = true)
    println("done")
    lpk = NfOrdIdl[ P[1] for p in cp for P = prime_decomposition(zk, p)]
    println("Now computing the S-unit group for lp of length $(length(lpk))")
    @assert length(lpk) > 0
    @time Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)
    @show length(N.coefficients_gen[i])

    #D = Dict{nf_elem, nf_elem}()
    cc = Vector{Tuple{Int, fmpz}}[]
    induced = induce_action_from_subfield(N, i, lpk, c.FB, cc)
    for l in 1:ngens(Szk)
      u = mS(Szk[l])
      for j in 1:length(induced)
        aut = autos[j]
        valofnewelement = mul(mS.valuations[l], induced[j])
        el_to_add = FacElem(Dict{nf_elem, AnticNumberField}((aut(embedding(N, i)(x)), v) for (x,v) = u.fac))
        Hecke.class_group_add_relation(c, el_to_add, valofnewelement)#, always = false)
      end
    end

    #for j in 1:length(N.coefficients_gen[i])
    #  println("Inducing the action ... ")
    #  @time z = induce_action(N, i, j, lpk, c.FB, cc)

    #  print("Feeding in the S-units of the small field ... ")

    #  for l=1:ngens(Szk)
    #    u = mS(Szk[l])  #do compact rep here???
    #    valofnewelement = mul(mS.valuations[l], z)
    #    Hecke.class_group_add_relation(c, FacElem(Dict((N(x, i, j), v) for (x,v) = u.fac)), valofnewelement)
    #  end
    #end

    println("done")

    u, mu = unit_group_fac_elem(zk)
    for n=2:ngens(u) # I cannot add a torsion unit. He will hang forever.
      uelem = mu(u[n])
      for j in 1:length(induced)
        aut = autos[j]
        lifted_unit = FacElem(Dict((aut(embedding(N, i)(x)), v) for (x,v) = uelem.fac))
        bb = Hecke.add_unit!(U, lifted_unit)
        @show bb
        if !bb
          push!(skipped_units, lifted_unit)
        end
      end
      U.units = Hecke.reduce(U.units, U.tors_prec)
    end

    #for j in 1:length(N.coefficients_gen[i])
    #  for n=2:ngens(u) # I cannot add a torsion unit. He will hang forever.
    #    uelem = mu(u[n])
    #    lifted_unit = FacElem(Dict((N(x, i, j), v) for (x,v) = uelem.fac))
    #    bb = Hecke._add_unit(U, lifted_unit)
    #    @show bb
    #    if !bb
    #      push!(skipped_units, lifted_unit)
    #    end
    #  end
    #  U.units = Hecke.reduce(U.units, U.tors_prec)
    #end
  end

  @show length(skipped_units)

  for lifted_unit in skipped_units
    @show Hecke.regulator(U.units, 64)
    Hecke._add_dependent_unit(U, lifted_unit)
  end

  U.tentative_regulator = Hecke.regulator(U.units, 64)

  println("Now saturating ... ")

  if saturate
    for (p, e) in factor(index(N))
      b = Hecke.saturate!(c, U, Int(p))
      while b
        b = Hecke.saturate!(c, U, Int(p))
      end
    end
  end
  return c, U
end

################################################################################
#
#  Redo the S-unit computation for Brauer norm relation
#
################################################################################

function _embed(N::NormRelation, i::Int, a::nf_elem)
  b = get(N.embed_cache_triv[i], a, nothing)
  if b === nothing
    _, mk = subfield(N, i)
    c = mk(a)
    N.embed_cache_triv[i][a] = c
    return c
  else
    return b
  end
end

# pure

function _add_sunits_from_brauer_relation!(c, UZK, N; invariant = false, compact::Int = 0)
  # I am assuming that c.FB.ideals is invariant under the action of the Galois group
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  add_unit_later = FacElem{nf_elem, AnticNumberField}[]
  for i = 1:length(N)
    if isdefined(N, :nonredundant) && !(i in N.nonredundant)
      continue
    end
    k, mk = subfield(N, i)
    @vprint :NormRelation 1 "Looking at the subfield $i / $(length(N)) with defining equation $(k.pol) \n"
    @vprint :NormRelation 1 "Computing lll basis of maximal order ...\n"
    zk = maximal_order(k)
    zk = lll(zk)
    @vprint :NormRelation 1 "Computing class group ... \n"
    class_group(zk, use_aut = true)
    lpk = NfOrdIdl[]
    for p in cp
      for P in prime_decomposition(zk, p)
        push!(lpk, P[1])
      end
    end
	  #lpk = unique!(NfOrdIdl[ intersect_prime(mk, P) for P in c.FB.ideals])
    @assert length(lpk) > 0
    @vprint :NormRelation 1 "Computing sunit group for set of size $(length(lpk)) ... \n"
    Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)

    @vprint :NormRelation 1 "Coercing the sunits into the big field ...\n"
    z = induce_action_just_from_subfield(N, i, lpk, c.FB, invariant)

    for l=1:ngens(Szk)
      @vprint :NormRelation 3 "Sunits $l / $(ngens(Szk))\n"
      u = mS(Szk[l])  #do compact rep here???
      if compact != 0
        @vprint :NormRelation 3 "  Compact presentation ...\n"
        u = Hecke.compact_presentation(u, compact)
      end
      valofnewelement = mul(mS.valuations[l], z)
      @hassert :NormRelation 1 begin zz = mk(evaluate(u)); true; end
      @hassert :NormRelation 1 sparse_row(FlintZZ, [ (i, valuation(zz, p)) for (i, p) in enumerate(c.FB.ideals) if valuation(zz, p) != 0]) == valofnewelement
      @vtime :NormRelation 4 img_u = FacElem(Dict{nf_elem, fmpz}((_embed(N, i, x), v) for (x,v) = u.fac))
      @vtime :NormRelation 4 Hecke.class_group_add_relation(c, img_u, valofnewelement)#, always = false)
    end

    @vprint :NormRelation 1 "Coercing the units into the big field ... \n"
    U, mU = unit_group_fac_elem(zk)
    for j=2:ngens(U) # I cannot add a torsion unit. He will hang forever.
      @vprint :NormRelation 3 "Unit $j / $(ngens(U))\n"
      u = mU(U[j])
      if compact != 0
        @vprint :NormRelation 3 "  Compact presentation ...\n"
        u = Hecke.compact_presentation(u, compact)
      end
      @vtime :NormRelation 4 img_u = FacElem(Dict{nf_elem, fmpz}((_embed(N, i, x), v) for (x,v) = u.fac))
      has_full_rk = Hecke.has_full_rank(UZK)
      @vtime :NormRelation 4 ff = Hecke.add_unit!(UZK, img_u)
      if !has_full_rk && !ff
        push!(add_unit_later, img_u)
      end
    end
    @vprint :NormRelation 4 "Reducing the units\n"
    @vtime :NormRelation 4 UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end

  if length(add_unit_later) > 0
    for uu in add_unit_later
      Hecke.add_unit!(UZK, uu)
    end
    UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end

  return nothing
end

# TODO: Encode that the i-th field is normal over Q
function induce_action_just_from_subfield(N::NormRelation, i, s, FB, invariant = false)
  S = FB.ideals
  ZK = order(S[1])

  z = zero_matrix(SMat, FlintZZ, 0, length(S))

  mk = embedding(N, i)
  zk = order(s[1])

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    v = Tuple{Int, fmpz}[]
    P = s[l]
    genofsl = elem_in_nf(_embed(N, i, P.gen_two.elem_in_nf))
    pmin = minimum(P, copy = false)
    # Compute the number of ideals
    # Assume that L/K and L/Q are normal
    rele = divexact(ramification_index((FB.fb[pmin].lp)[1][2]), ramification_index(P))
    relf = divexact(degree((FB.fb[pmin].lp)[1][2]), degree(P))
    @assert mod(reldeg, rele * relf) == 0
    numb_ideal = divexact(reldeg, rele * relf)
    found = 0
    for k in 1:length(S)
      Q = S[k]
      if minimum(Q, copy = false) == pmin
        if genofsl in Q
          found += 1
          @assert mod(ramification_index(Q), ramification_index(s[l])) == 0
          push!(v, (k, divexact(ramification_index(Q), ramification_index(s[l]))))
        end
      end
      if invariant && N.isnormal[i]
				if found == numb_ideal
					break
				end
			end
    end
    sort!(v, by = x -> x[1])
    push!(z, sparse_row(FlintZZ, v))
  end
  return z
end

function norm_relation(K::AnticNumberField, coprime::Int = 0; small_degree = true)
  if !small_degree
    @assert coprime == 0
    N = _norm_relation_setup_generic(K, pure = true, small_degree = false)
    return true, N
  end
  local N
  try
    N = _get_nf_norm_relation(K)::Vector{NormRelation{Int}}
    if coprime == 0
      return true, N[1]::NormRelation{Int}
    else
      for i in 1:length(N)
        if coprime == N[i].isoptimal
          return true, N[i]::NormRelation{Int}
        end
      end
      fl, M = has_coprime_norm_relation(K, fmpz(coprime))
      if fl
        push!(N, M)
        return true, M
      end
      return false, NormRelation{Int}()
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    if coprime == 0
      if !has_useful_brauer_relation(automorphism_group(K)[1])
        return false, NormRelation{Int}()
      end
      M = _norm_relation_setup_generic(K, pure = true, small_degree = true)
      _set_nf_norm_relation(K, NormRelation{Int}[M])
      return true, M::NormRelation{Int}
    else
      fl, M = has_coprime_norm_relation(K, fmpz(coprime))
      if fl
        _set_nf_norm_relation(K, NormRelation{Int}[M])
        return true, M::NormRelation{Int}
      end
      return false, NormRelation{Int}()
    end
  end
end

function _sunit_group_fac_elem_quo_via_brauer(K::AnticNumberField, S, n::Int, invariant::Bool = false)
  @vprint :NormRelation 1 "Setting up the norm relation context ... \n"
  fl, N = norm_relation(K, n)
  if !fl
    fl, N = norm_relation(K, 0)
  end
  @assert fl
  @vprint :NormRelation 1 "Using norm relation $N\n"
  return __sunit_group_fac_elem_quo_via_brauer(N, S, n, invariant)
end

function sunit_group_fac_elem_via_brauer(N::NormRelation, S, invariant::Bool = false)
  return _sunit_group_fac_elem_quo_via_brauer(N, S, 0, invariant)
end

function _class_group(O::NfOrd; saturate::Bool = true, redo::Bool = false)
  if !redo
    c = Hecke._get_ClassGrpCtx_of_order(O, false)
    if c !== nothing
      c::Hecke.ClassGrpCtx{SMat{fmpz}}
      @assert c.finished
      return class_group(c, O)
    end
  end
  K = nf(O)
  G, mG = automorphism_group(K)
  if degree(K) > 25 && has_useful_brauer_relation(G)
    lll(O)
    fl, N = norm_relation(K, small_degree = false)
    @assert fl
    return class_group_via_brauer(O, N, recursive = true, compact = false, saturate = saturate)
  end
  return class_group(O, use_aut = true)
end

function get_sunits_from_subfield_data(OK, N; recursive::Bool = false, compact::Bool = true)
  K = nf(OK)
  fbbound = Hecke.factor_base_bound_grh(OK)
  S = prime_ideals_up_to(OK, fbbound)
  @vprint :NormRelation 1 "Factor base bound: $fbbound\n"
  c, UZK = _setup_for_norm_relation_fun(K, S)

  docompact = false
  onlyp = 0

  if Hecke.isprime_power(index(N)) && compact
    docompact = true
    _, onlyp = ispower(index(N))
  end

  if recursive
    for i = 1:length(N.subfields)
      if !(i in N.nonredundant)
        continue
      end
      L = N.subfields[i][1]
      if N.isnormal[i] && degree(L) > 25
        @vprint :NormRelation 1 "Computing class group of $(L.pol)\n"
        _class_group(lll(maximal_order(L)), saturate = false)
      end
    end
  end

  @vprint :NormRelation 1 "Doing something in the subfields\n"
  if !docompact
    _add_sunits_from_brauer_relation!(c, UZK, N)
  else
    @vprint :NormRelation 1 "Using the compact presentation\n"
    _add_sunits_from_brauer_relation!(c, UZK, N, compact = onlyp)
  end
  return c, UZK
end

function class_group_via_brauer(O::NfOrd, N::NormRelation; saturate::Bool = true, recursive::Bool = false, do_lll::Bool = true, compact::Bool = true, stable = 3.5)
  K = N.K
  if do_lll
    OK = lll(maximal_order(nf(O)))
  else
    OK = O
  end

  c, UZK = get_sunits_from_subfield_data(OK, N, recursive = recursive, compact = compact)
  #auts = automorphisms(K)
  #c.aut_grp = Hecke.class_group_add_auto(c, auts)

  if saturate || isone(index(N))
    for (p, e) in factor(index(N))
      @vprint :NormRelation 1 "Saturating at $p \n"
      b = Hecke.saturate!(c, UZK, Int(p), stable)
      while b
        idx = Hecke._validate_class_unit_group(c, UZK)[1]
        @vprint :NormRelation 1 "Index bound from analysis $idx\n"
        b = Hecke.saturate!(c, UZK, Int(p), stable)
      end
    end
    idx = Hecke._validate_class_unit_group(c, UZK)[1]
    @vprint :NormRelation 1 "Index is $idx (should be 1)!\n"
    if idx != 1
      @vprint :NormRelation 1 "Index is $idx (should be 1)!\n"
    end
    @assert idx == 1

    @vprint :NormRelation 1 "\n"
    c, _ = simplify(c)

    c.finished = true
    UZK.finished = true
    Hecke._set_ClassGrpCtx_of_order(OK, c)
    Hecke._set_UnitGrpCtx_of_order(OK, UZK)
    return class_group(c, O)
  else
    # We finish off the computation by a relation search in OK
    Hecke._set_ClassGrpCtx_of_order(OK, c)
    Hecke._set_UnitGrpCtx_of_order(OK, UZK)
    c, UZK, d = Hecke._class_unit_group(OK, saturate_at_2 = false)
    @assert isone(d)

    return Hecke.class_group(c, O)
  end
end

function __sunit_group_fac_elem_quo_via_brauer(N::NormRelation, S, n::Int, invariant::Bool = false)
  O = order(S[1])

  K = N.K

  if invariant
    c, UZK = Hecke._setup_for_norm_relation_fun(K, S)
    _add_sunits_from_brauer_relation!(c, UZK, N, invariant = true)
  else
    cp = sort!(collect(Set(minimum(x) for x = S)))
    Sclosed = NfOrdIdl[]
    for p in cp
      lp = prime_decomposition(O, p)
      for (P, _) in lp
        push!(Sclosed, P)
      end
    end
    if length(Sclosed) == length(S)
      invariant = true
      Sclosed = S
    end
    @vprint :NormRelation 1 "I am not Galois invariant. Working with S of size $(length(Sclosed))\n"
    c, UZK = _setup_for_norm_relation_fun(K, Sclosed)
    _add_sunits_from_brauer_relation!(c, UZK, N, invariant = true)
  end

  if gcd(index(N), n) != 1
    # I need to saturate
    @vprint :NormRelation 1 "Saturating at "
    for (p, e) in factor(index(N))
      @vprint :NormRelation 1 "$p "
      b = Hecke.saturate!(c, UZK, Int(p))
      while b
        b = Hecke.saturate!(c, UZK, Int(p))
      end
    end
  end
  @vprint :NormRelation 1 "\n"

  # This makes c.R.gen be a basis of the S-units (modulo torsion)
  c, _ = simplify(c)

  if invariant
    sunitsmodunits = c.R_gen # These are generators for the S-units (mod units, mod n)
  else
    # I need to extract the S-units from the Sclosed-units
    # Now I need to find the correct indices in the c.FB.ideals
    sunitsmodunits = typeof(c.R_gen)()
    ind = Int[]
    for P in S
      for i in 1:length(c.FB.ideals)
        if P == c.FB.ideals[i]
          push!(ind, i)
          break
        end
      end
    end
    sort!(ind)
    # ind = indices of S inside c.FB.ideals
    @assert length(Sclosed) == length(c.FB.ideals)
    @assert length(ind) == length(S)
    z = zero_matrix(FlintZZ, length(c.R_gen), length(Sclosed) - length(S))
    for i in 1:length(c.R_gen)
      k = 1
      for j in 1:length(Sclosed)
        if !(j in ind)
          z[i, k] = c.M.bas_gens[i, j]
          if k == ncols(z)
            break
          end
          k = k + 1
        end
      end
    end
    k, K = left_kernel(z)
    for i in 1:nrows(K)
      if iszero_row(K, i)
        continue
      end
      push!(sunitsmodunits, FacElem(c.R_gen, fmpz[K[i, j] for j in 1:ncols(K)]))
    end
  end

  unitsmodtorsion = UZK.units # These are generators for the units (mod n)
  T, mT = torsion_unit_group(O)
  Q, mQ = quo(T, n, false)
  @assert issnf(Q)
  @assert ngens(Q) == 1
  m = order(Q)

  if !isone(m)
    tomodn = FacElem(elem_in_nf(mT(mQ\Q[1])))
    res_group = abelian_group(append!(fmpz[m], [fmpz(n) for i in 1:(length(sunitsmodunits) + length(unitsmodtorsion))]))

    local exp
    let res_group = res_group, unitsmodtorsion = unitsmodtorsion, tomodn = tomodn
      function exp(a::GrpAbFinGenElem)
        @assert parent(a) == res_group
        zz = FacElem(convert(Vector{Hecke.nf_elem_or_fac_elem}, unitsmodtorsion), fmpz[a[1 + i] for i in 1:length(unitsmodtorsion)])
        z = mul!(zz, zz, tomodn^a[1])
        zzz = FacElem(convert(Vector{Hecke.nf_elem_or_fac_elem}, sunitsmodunits), fmpz[a[1 + length(unitsmodtorsion) + i] for i in 1:length(sunitsmodunits)])
        mul!(z, z, zzz)

        for (k, v) in z.fac
          if iszero(v)
            delete!(z.fac, k)
          end
        end
        return z
      end
    end

    disclog = function(a)
      throw(NotImplemented())
    end
  else # torsion part is one
    res_group = abelian_group([fmpz(n) for i in 1:(length(sunitsmodunits) + length(unitsmodtorsion))])

    exp = function(a::GrpAbFinGenElem)
      @assert parent(a) == res_group
      z = FacElem(convert(Vector{nf_elem_or_fac_elem}, unitsmodtorsion), fmpz[a[i] for i in 1:length(unitsmodtorsion)])
      # This is madness
      zz = FacElem(convert(Vector{nf_elem_or_fac_elem}, sunitsmodunits), fmpz[a[length(unitsmodtorsion) + i] for i in 1:length(sunitsmodunits)])
      z = mul!(z, z, zz)

      for (k, v) in z.fac
        if iszero(v)
          delete!(z.fac, k)
        end
      end

      return z
    end

    disclog = function(a)
      throw(NotImplemented())
    end
  end

  r = Hecke.MapSUnitGrpFacElem()
  r.idl = S
  r.isquotientmap = n

  r.header = Hecke.MapHeader(res_group, FacElemMon(nf(O)), exp, disclog)
  @hassert :NormRelation 9000 begin
    _S, _mS = sunit_group_fac_elem(S)
    _Q, _mQ = quo(_S, n)
    V = quo(_Q, [_mQ(_mS\(r(res_group[i]))) for i in 1:ngens(res_group)])
    order(_Q) == order(res_group) && order(V[1]) == 1
  end
  return res_group, r
end


text:
function solve()
    n = parse(Int, readline())
    as = [parse(Int, x) for x in split(readline())]

    ret = 0
    pos = 0
    m = 0
    cs = 0

    for a in as
        cs += a
        m = max(m, cs)
        ret = max(ret, pos+m)
        pos = pos + cs
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module CIAOAlgorithms #how we use this module? always paths are saved in a same format/ any other formats? how to switch to dev mode? in the dev mode, what happens to the other packages? how to add a package to the dev mode?

const RealOrComplex{R} = Union{R,Complex{R}}
const Maybe{T} = Union{T, Nothing}

# utilities
include("utilities/indexingUtilities.jl")
include("utilities/IndexIterator.jl")
include("utilities/lbfgs.jl")
# include("utilities/indNonnegativeBallL2.jl")
# algorithms
include("algorithms/Finito/Finito.jl")
include("algorithms/ProShI/ProShI.jl")
include("algorithms/SVRG/SVRG.jl")
include("algorithms/SGD/SGD.jl")
include("algorithms/SARAH/SARAH.jl")
include("algorithms/SAGA_SAG/SAGA.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "alsa_plugins"
version = v"1.2.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("ftp://ftp.alsa-project.org/pub/plugins/alsa-plugins-1.2.2.tar.bz2", "1c0f06450c928d711719686c9dbece2d480184f36fab11b8f0534cb7b41e337d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd alsa-plugins-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libasound_module_pcm_a52", :libasound_module_pcm_a52, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_pulse", :libasound_module_pcm_pulse, "lib/alsa-lib"),
    LibraryProduct("libasound_module_rate_lavrate", :libasound_module_rate_lavrate, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_speex", :libasound_module_pcm_speex, "lib/alsa-lib"),
    LibraryProduct("libasound_module_conf_pulse", :libasound_module_conf_pulse, "lib/alsa-lib"),
    LibraryProduct("libasound_module_ctl_pulse", :libasound_module_ctl_pulse, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_oss", :libasound_module_pcm_oss, "lib/alsa-lib"),
    LibraryProduct("libasound_module_rate_samplerate", :libasound_module_rate_samplerate, "lib/alsa-lib"),
    LibraryProduct("libasound_module_ctl_arcam_av", :libasound_module_ctl_arcam_av, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_usb_stream", :libasound_module_pcm_usb_stream, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_vdownmix", :libasound_module_pcm_vdownmix, "lib/alsa-lib"),
    LibraryProduct("libasound_module_pcm_upmix", :libasound_module_pcm_upmix, "lib/alsa-lib"),
    LibraryProduct("libasound_module_rate_speexrate", :libasound_module_rate_speexrate, "lib/alsa-lib"),
    LibraryProduct("libasound_module_ctl_oss", :libasound_module_ctl_oss, "lib/alsa-lib")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFMPEG_jll", uuid="b22a6f82-2f65-5046-a5b2-351ab43fb4e5"))
    Dependency(PackageSpec(name="alsa_jll", uuid="45378030-f8ea-5b20-a7c7-1a9d95efb90e"))
    Dependency(PackageSpec(name="libsamplerate_jll", uuid="9427e74d-4e05-59c1-8ff3-7d74b6e52ac8"))
    Dependency(PackageSpec(name="PulseAudio_jll", uuid="02771fc1-bdb7-5db5-8d11-300768e00fbd"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Fitting" begin
    ε = 1.0e-8

    @testset "Fitting-curve_R" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        a_org = [Point(i1, rand()) for i1 in 1:n1]
        M = BSplineCurve([P1], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)

        M′ = refinement(M, [P1′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′])
        @test norm(a_fit - a_ref) < ε

        a_fit = fittingcontrolpoints(M, [P1′], domain=:R)
        @test norm(a_fit - a_ref) < ε
    end

    @testset "Fitting-curve_I" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        a_org = [Point(i1, rand()) for i1 in 1:n1]
        M = BSplineCurve([P1], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)

        M′ = refinement(M, [P1′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′])

        @test norm(a_fit - a_ref) < ε
    end

    @testset "Fitting-surface_R" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        p2 = 1
        k2 = Knots(rand(4)) + (p2 + 1) * Knots([0, 1])
        P2 = FastBSplineSpace(p2, k2)
        n2 = dim(P2)
        a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
        M = BSplineSurface([P1, P2], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)
        p2′ = p2 + 1
        k2′ = k2 + unique(k2) + Knots(rand(2))
        P2′ = FastBSplineSpace(p2′, k2′)

        M′ = refinement(M, [P1′, P2′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′, P2′])
        @test norm(a_fit - a_ref) < ε

        a_fit = fittingcontrolpoints(M, [P1′, P2′], domain=:R)
        @test norm(a_fit - a_ref) < ε
    end

    @testset "Fitting-surface_I" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        p2 = 1
        k2 = Knots(rand(4)) + Knots([0, 1]) + p2 * Knots([-rand(), 1 + rand()])
        P2 = FastBSplineSpace(p2, k2)
        n2 = dim(P2)
        a_org = [Point(i1, i2, rand()) for i1 in 1:n1, i2 in 1:n2]
        M = BSplineSurface([P1, P2], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)
        p2′ = p2 + 1
        k2′ = k2 + unique(k2[1+p2:end-p2]) + Knots(rand(2))
        P2′ = FastBSplineSpace(p2′, k2′)

        M′ = refinement(M, [P1′, P2′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′, P2′])
        @test norm(a_fit - a_ref) < ε
    end

    @testset "Fitting-solid_R" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        p2 = 1
        k2 = Knots(rand(4)) + (p2 + 1) * Knots([0, 1])
        P2 = FastBSplineSpace(p2, k2)
        n2 = dim(P2)
        p3 = 2
        k3 = Knots(rand(5)) + (p3 + 1) * Knots([0, 1])
        P3 = FastBSplineSpace(p3, k3)
        n3 = dim(P3)
        a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
        M = BSplineSolid([P1, P2, P3], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)
        p2′ = p2 + 1
        k2′ = k2 + unique(k2) + Knots(rand(2))
        P2′ = FastBSplineSpace(p2′, k2′)
        p3′ = p3 + 1
        k3′ = k3 + unique(k3) + Knots(rand(2))
        P3′ = FastBSplineSpace(p3′, k3′)

        M′ = refinement(M, [P1′, P2′, P3′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′])
        @test norm(a_fit - a_ref) < ε

        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′], domain=:R)
        @test norm(a_fit - a_ref) < ε
    end

    @testset "Fitting-solid_I" begin
        Random.seed!(42)

        p1 = 2
        k1 = Knots(rand(3)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
        P1 = FastBSplineSpace(p1, k1)
        n1 = dim(P1)
        p2 = 1
        k2 = Knots(rand(4)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
        P2 = FastBSplineSpace(p2, k2)
        n2 = dim(P2)
        p3 = 2
        k3 = Knots(rand(5)) + Knots([0, 1]) + p1 * Knots([-rand(), 1 + rand()])
        P3 = FastBSplineSpace(p3, k3)
        n3 = dim(P3)
        a_org = [Point(i1, i2, i3, rand()) for i1 in 1:n1, i2 in 1:n2, i3 in 1:n3]
        M = BSplineSolid([P1, P2, P3], a_org)

        p1′ = p1 + 1
        k1′ = k1 + unique(k1[1+p1:end-p1]) + Knots(rand(2))
        P1′ = FastBSplineSpace(p1′, k1′)
        p2′ = p2 + 1
        k2′ = k2 + unique(k2[1+p2:end-p2]) + Knots(rand(2))
        P2′ = FastBSplineSpace(p2′, k2′)
        p3′ = p3 + 1
        k3′ = k3 + unique(k3[1+p3:end-p3]) + Knots(rand(2))
        P3′ = FastBSplineSpace(p3′, k3′)

        M′ = refinement(M, [P1′, P2′, P3′])
        a_ref = M′.controlpoints

        a_fit = fittingcontrolpoints(M, [P1′, P2′, P3′])
        @test norm(a_fit - a_ref) < 1.0e-6
    end
end


text:
## extremes.jl  Test some extreme cases
## (c) 2014--2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

## check some pathological cases where sorting and comparison of equal values plays a role

using Base.Test

function test_both(tar, non, e1, e2)
    r = roc(tar, non)
#    rocplot(r)
    eer1 = eerch(r)
    eer2 = eer(tar, non)
    println("EERch = ", eer1)
    println("EER   = ", eer2)
    @test_approx_eq e1 eer1
    @test_approx_eq e2 eer2
end

test_both([0;2], [1], 1/3, 1/2)
test_both([1], [0;2], 1/3, 1/2)
test_both([0], [1], 1/2, 1)
test_both([1], [0], 0, 0)

test_both(collect(0:2), [1], 2/5, 1/2)
test_both(collect(0.:2), [0.5], 1/4, 1/3)
test_both(collect(0.:2), [1.5], 2/5, 2/3)

test_both([1], collect(0:2), 2/5, 1/2)
test_both([0.5], collect(0.:2), 2/5, 2/3)
test_both([1.5], collect(0.:2), 1/4, 1/3)


text:
export DoubleFloat64, FastDouble, AccurateDouble, sincosh,
double_square, double_add, double_mul, double_div, double_sub, double_sqrt,
square

"""
    DoubleFloat64(x [, mode::ComputeMode]) <: AbstractFloat
"""
struct DoubleFloat64{T<:ComputeMode} <: AbstractFloat
    hi::Float64
    lo::Float64
end

"""
    FastDouble

Shorthand for `DoubleFloat64{ComputeFast}`
"""
const FastDouble = DoubleFloat64{ComputeFast}

"""
    AccurateDouble

Shorthand for `DoubleFloat64{ComputeAccurate}`
"""
const AccurateDouble = DoubleFloat64{ComputeAccurate}

DoubleFloat64(x::DoubleFloat64{T}) where T = DoubleFloat64{T}(x.hi, x.lo)
function DoubleFloat64(x::Float64, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64{compmode}(x, isinf(x) ? Inf : 0.0)
end
function DoubleFloat64(x::Float64, y::Float64, compmode::Type{<:ComputeMode}=ComputeFast)
    FastDouble(x, y)
end
function DoubleFloat64(x::Float32, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64{compmode}(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleFloat64(x::Float16, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64{compmode}(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleFloat64(x::Integer, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64{compmode}(convert(Float64, x), isinf(x) ? Inf : 0.0)
end
function DoubleFloat64(x::BigFloat, compmode::Type{<:ComputeMode}=ComputeFast)
    z = convert(Float64, x)
    DoubleFloat64{compmode}(z, convert(Float64, x-z))
end
function DoubleFloat64(x::Irrational, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64(big(x), compmode)
end
function DoubleFloat64(x::Rational, compmode::Type{<:ComputeMode}=ComputeFast)
    DoubleFloat64(numerator(x), compmode) / DoubleFloat64(denominator(x), compmode)
end

hi(x::DoubleFloat64) = x.hi
lo(x::DoubleFloat64) = x.lo

Base.isbits(::DoubleFloat64) = true

Base.zero(::DoubleFloat64{T}) where T = DoubleFloat64{T}(0.0, 0.0)
Base.zero(::Type{DoubleFloat64{T}}) where T = DoubleFloat64{T}(0.0, 0.0)

Base.one(::DoubleFloat64{T}) where T = DoubleFloat64{T}(1.0, 0.0)
Base.one(::Type{DoubleFloat64{T}}) where T = DoubleFloat64{T}(1.0, 0.0)

Base.convert(::Type{T}, a::DoubleFloat64) where {T<:AbstractFloat} = convert(T, a.hi)
Base.convert(::Type{BigFloat}, a::DoubleFloat64) = big(a.hi) + big(a.lo)
Base.convert(::Type{T}, a::DoubleFloat64) where {T<:Integer} = convert(T, a.hi)
Base.convert(::Type{Integer}, a::DoubleFloat64) = convert(Int64, a.hi)
Base.convert(::Type{BigInt}, a::DoubleFloat64) = convert(BigInt, big(a.hi) + big(a.lo))

Base.convert(::Type{FastDouble}, x::AccurateDouble) = FastDouble(x.hi, x.lo)
Base.convert(::Type{AccurateDouble}, x::FastDouble) = AccurateDouble(x.hi, x.lo)
Base.convert(::Type{FastDouble}, x::FastDouble) = x
Base.convert(::Type{DoubleFloat64}, x::FastDouble) = x
Base.convert(::Type{AccurateDouble}, x::AccurateDouble) = x
Base.convert(::Type{DoubleFloat64}, x::AccurateDouble) = x
Base.convert(::Type{DoubleFloat64{T}}, x::AbstractFloat) where {T<:ComputeMode} = DoubleFloat64(x, T)
Base.convert(::Type{DoubleFloat64{T}}, x::Irrational) where T = DoubleFloat64(x, T)
Base.convert(::Type{DoubleFloat64{T}}, x::Integer) where T = DoubleFloat64(x, T)
Base.convert(::Type{DoubleFloat64}, x::AbstractFloat) = DoubleFloat64(x, ComputeFast)
Base.convert(::Type{DoubleFloat64}, x::Irrational) = DoubleFloat64(x, ComputeFast)
Base.convert(::Type{DoubleFloat64}, x::Integer) = DoubleFloat64(x, ComputeFast)

Base.promote_rule(::Type{FastDouble}, ::Type{DoubleFloat64{ComputeAccurate}}) = DoubleFloat64{ComputeAccurate}
Base.promote_rule(::Type{DoubleFloat64{T}}, ::Type{I}) where {T, I<:Integer} = DoubleFloat64{T}
Base.promote_rule(::Type{DoubleFloat64}, ::Type{<:Integer}) = FastDouble
Base.promote_rule(::Type{DoubleFloat64{T}}, ::Type{BigInt}) where {T} = BigFloat
Base.promote_rule(::Type{DoubleFloat64}, ::Type{BigFloat}) = BigFloat
Base.promote_rule(::Type{DoubleFloat64{T}}, ::Type{Float64}) where T = DoubleFloat64{T}
Base.promote_rule(::Type{DoubleFloat64{T}}, ::Type{Float32}) where T = DoubleFloat64{T}
Base.promote_rule(::Type{DoubleFloat64{T}}, ::Type{Float16}) where T = DoubleFloat64{T}
Base.promote_rule(::Type{DoubleFloat64}, ::Type{Float64}) = FastDouble
Base.promote_rule(::Type{DoubleFloat64}, ::Type{Float32}) = FastDouble
Base.promote_rule(::Type{DoubleFloat64}, ::Type{Float16}) = FastDouble

Base.big(x::DoubleFloat64) = big(x.hi) + big(x.lo)

include("double_const.jl")
include("double_arithmetic.jl")
include("double_misc.jl")
include("double_funcs.jl")

function Base.show(io::IO, x::DoubleFloat64)
    @printf io "%.32g" big(x)  # crude approximation to valid number of digits
end


text:
#========== backward gradient using symbolic derivatives ==========#

using DiffRules

function insert_symbolic_gradient(axislist, store)

    dZ = Symbol(DEL, ZED)
    ∇act! = Symbol(:∇, ACT!)
    gradarrays = map(A -> Symbol(DEL, A), store.arrays)
    # gradscalars = map(A -> Symbol(DEL, A), store.scalars)

    out_ind, in_ind = if store.redfun == :+
        store.sharedind, setdiff(vcat(store.leftind, store.redind), store.sharedind)
    elseif store.redfun in [:min, :max] # :*,
        store.leftind, store.redind
    else
        error("can't take gradients with reduction $(store.redfun)")
    end

    targets = []
    MacroTools_postwalk(symbwalk(targets), store.right)
    # append!(targets, scalars)

    inbody, prebody = [], []
    for (dt, t) in unique(targets)
        drdt = leibnitz(store.right, t)
        deltar = simplitimes(drdt, :(conj($dZ[$(store.leftraw...)])))
        if store.redfun == :+
            push!(inbody, :($dt = $dt + conj($deltar)))
        # elseif store.redfun == :*
        #     push!(inbody, :($dt = conj($deltar) * $ZED[$(store.leftraw...)] * inv($(store.right))))
        #     push!(prebody, :($dt = conj($deltar) * $ACC))
        elseif store.redfun in [:min, :max]
            push!(inbody, :($dt += $deltar)) # only when max attained, i.e. rhs == lhs!
        end
    end
    store.verbose>0 && @info "symbolic gradients" inbody
    ex_body = commonsubex(quote $(inbody...) end)

    ex_pre, ex_post = if store.redfun == :* # then nonzero LHS are handled already, but harder cases here:
        product_grad(prebody, store)
    else
        nothing, nothing
    end
    if store.redfun in [:min, :max] # this case really wants sparse 𝛥x!
        ex_body = :(
            if $ZED[$(store.leftraw...)] == $(store.right)
                $ex_body
            end
        )
    end

    make_many_actors(∇act!,
        vcat(gradarrays, :($dZ::AbstractArray{$TYP}), ZED, store.arrays, store.scalars, axislist),
        # vcat(gradarrays, gradscalars, :($dZ::AbstractArray{$TYP}), store.arrays, store.scalars, axislist),
        nothing, out_ind, ex_pre, in_ind, ex_body, ex_post, store, " (symbolic gradient)")

    if isdefined(store.mod, :Zygote) # special case for FillArrays
        ex_body2 = fillarrayreplace(ex_body, dZ)
        ex_pre2 = fillarrayreplace(ex_pre, dZ)
        ex_value = :($(Symbol(dZ, :_value)) = $dZ.value) # @avx likes this outside the loop

        make_many_actors(∇act!,
            vcat(gradarrays, :($dZ::Zygote.Fill{$TYP}), ZED, store.arrays, store.scalars, axislist),
            ex_value, out_ind, ex_pre2, in_ind, ex_body2, ex_post, store, " (method for FillArrays)")
    end

end

# This works for simple cases, but the general case is more complicatd.
#=
product_grad(prebody, store) = begin
    cnt = Symbol(DEL,:𝒸ℴ𝓊𝓃𝓉,0)

    inds_orig = :(($(store.redind...),))
    inds_prime = :(($(map(i -> Symbol(i,'′',DEL), store.redind)...),))
    inds_zero = :(($(map(i -> 0, store.redind)...),))

    loop_search = recurseloops(:(
        # find and save the index at which RHS is zero
        if iszero($(store.right))
            $cnt += 1
            $inds_prime = $inds_orig
        end
    ), copy(store.redind))

    loop_accum = recurseloops(:(
        # product of RHS at all redind except the one which gives zero
        $ACC = $ACC * ifelse($inds_orig == $inds_prime, 1, $(store.right))
    ), copy(store.redind))

    store.verbose>0 && @info "symbolic gradients extra..." prebody
    ex_prebody = commonsubex(quote $(prebody...) end)

    ex_pre = quote
        if iszero($ZED[$(store.leftraw...)])
            local $cnt = 0
            local $inds_prime = $inds_zero
            $loop_search
            if $cnt == 1
                local $ACC = one($TYP)
                $loop_accum
                let $inds_orig = $inds_prime
                    $ex_prebody
                end
            end # elseif more than one zero, then leave 𝛥x .== 0
            # continue # i.e. skip the ordinary routine, which divides
            @goto JUMP
        end
    end

    ex_post = quote
        @label JUMP
    end

    push!(store.notfree, cnt) # hack to disable @inbounds, avoids ERROR: syntax: misplaced label

    ex_pre, ex_post
end
=#

#========== symbolic differentiation ==========#

# This could probably use https://github.com/dfdx/XGrad.jl
# or https://github.com/SciML/ModelingToolkit.jl
# or https://github.com/JuliaMath/Calculus.jl/blob/master/src/differentiate.jl
# or now I found this: https://github.com/HarrisonGrodin/Simplify.jl
# but seemed simple enough to just write out, using rules from:
# http://www.juliadiff.org/DiffRules.jl/latest/

symbwalk(targets) = ex -> begin
        @capture_(ex, A_[inds__]) && A isa Symbol || return ex
        deltaex = :($(Symbol(DEL, A))[$(inds...)])
        push!(targets, (deltaex, ex))
        return ex
    end

leibnitz(s::Number, target) = 0
leibnitz(s::Symbol, target) = s == target ? 1 : 0
leibnitz(ex::Expr, target) = begin
    ex == target && return 1
    @capture_(ex, B_[ijk__]) && return 0
    if ex.head == Symbol("'")
        ex.head = :call
        pushfirst!(ex.args, :adjoint)
    end
    ex.head == :call || error("expected a functionn call, got $ex.")
    fun = ex.args[1]
    if fun == :log # catch log(a*b) and especially log(a/b)
        arg = ex.args[2]
        if arg isa Expr && arg.args[1] == :* && length(arg.args) == 3
            newex = :(log($(arg.args[2])) + log($(arg.args[3])))
            return leibnitz(newex, target)
        elseif arg isa Expr && arg.args[1] == :/
            newex = :(log($(arg.args[2])) - log($(arg.args[3])))
            return leibnitz(newex, target)
        end
    end
    if length(ex.args) == 2 # one-arg function
        fx = mydiffrule(fun, ex.args[2])
        dx = leibnitz(ex.args[2], target)
        return simplitimes(fx, dx)
    elseif length(ex.args) == 3  # two-arg function
        fx, fy = mydiffrule(fun, ex.args[2:end]...)
        dx = leibnitz(ex.args[2], target)
        dy = leibnitz(ex.args[3], target)
        return simpliplus(simplitimes(fx, dx), simplitimes(fy, dy))
    elseif fun in [:+, :*]
        fun == :* && return leibnitz(:(*($(ex.args[2]), *($(ex.args[3:end]...)))), target)
        dxs = [leibnitz(x, target) for x in ex.args[2:end]]
        fun == :+ && return simpliplus(dxs...)
    end
    error("don't know how to handle $ex.")
end

simplitimes(x::Number, y::Number) = x*y
simplitimes(x::Number, y) = x==0 ? 0 : x==1 ? y : x==-1 ? :(-$y) : :($x * $y)
simplitimes(x, y::Number) = y==0 ? 0 : y==1 ? x : y==-1 ? :(-$x) : :($y * $x)
simplitimes(x, y) = :($y * $x)

simpliplus(x::Number, y::Number) = x + y
simpliplus(x::Number, y) = x==0 ? y : :($x + $y)
simpliplus(x, y::Number) = y==0 ? x : :($x + $y)
simpliplus(x, y) = :($x + $y)
simpliplus(x, y, zs...) = simpliplus(simpliplus(x, y), zs...)

mydiffrule(f, xs...) = begin
    f == :+ && return map(_->1, xs)
    f == :- && return length(xs)==1 ? -1 : (1,-1)
    f == :^ && return mypowrule(xs...)
    f == :/ && return mydivrule(xs...)
    f == :// && return mydivrule(xs...)
    f == :inv && return mydivrule(1, xs...)[2]
    f == :log && return simpliinv(xs...)
    f == :sqrt && return mysqrtrule(xs...)
    f in BASE_NOGRAD && return map(_->0, xs)
    DiffRules.hasdiffrule(:Base, f, length(xs)) &&
        return DiffRules.diffrule(:Base, f, xs...)
    DiffRules.hasdiffrule(:SpecialFunctions, f, length(xs)) &&
        return DiffRules.diffrule(:SpecialFunctions, f, xs...)
    error("no diffrule found for function $f($(join(map(_->"_",xs),", "))).")
end

BASE_NOGRAD = [:(==), :(!=), :(<), :(<=), :(>), :(>=), :trunc, :round]

# Goals of these rules, besides correctness, are:
# 1. don't cause promotion of Float32, e.g. by factors (1/2)
# 2. make it easy for commonsubex(), e.g. by re-using inv(x)

mydivrule(x, y) = begin # (:(one(x) / y), :(-((x / y) / y)))
    invy = simpliinv(y)
    invy, :( -($x) * $invy * $invy )
end
mydivrule(x, y::Integer) = (y==1 ? 1 : 1//y), 0
mydivrule(x, y::Number) = (y==1 ? 1 : :(one($TYP)/$y)), 0

mydivrule(x::Number, y) = 0, :((-$x)*inv($y)*inv($y))
mydivrule(x::Number, y::Number) = 0, 0
mydivrule(x::Number, y::Integer) = 0, 0

mysqrtrule(x::Number) = sqrt(x)
mysqrtrule(x) = :(inv(sqrt($x))/2)

simpliinv(x) = :(inv($x))
simpliinv(x::Integer) = :(1//$x)
simpliinv(x::Expr) = if x.head == :call && x.args[1] == :/
        :($(x.args[3]) / $(x.args[2]))
    else
        :(inv($x))
    end

mypowrule(x, p) = begin
    dx = simplitimes(p, simplipow(x, simpliplus(p, -1)))
    dp = simplitimes(simplipow(x,p), :(log($x)))
    dx, dp
end

simplipow(x::Number, p::Number) = x^p
simplipow(x, p::Number) = p==1 ? x : p==2 ? :($x*$x) : :($x^$p)
simplipow(x, p) = :($x^$p)

#========== CSE ==========#

# My approach was to look for things occuring twice, biggest first.
# Then I found https://github.com/rdeits/CommonSubexpressions.jl
# which just pulls everything out, but doesn't like indexing expressions.

function commonsubex(expr::Expr)
    dict, defs, nope = Dict(), [], Set()
    if expr.head == :block
        args = [csewalk(ex, dict, defs, nope) for ex in copy(expr).args]
        quote
            $(defs...)
            $(args...)
        end
    else
        res = csewalk(copy(expr), dict, defs, nope)
        quote
            $(defs...)
            $res
        end
    end
end

csewalk(ex, dict, defs, nope) = ex
csewalk(ex::Expr, dict::Dict, defs::Vector, nope::Set) =
    # The goal is to alter RHS of assignments,
    # this mess is the most common case, A = A + stuff
    if ex.head == :(=) && ex.args[2] isa Expr && ex.args[2].head == :call &&
        ex.args[2].args[1] == :+ && ex.args[2].args[2] == ex.args[1]
        for n in 3:length(ex.args[2].args)
            ex.args[2].args[n] = csewalk(ex.args[2].args[n], dict, defs, nope)
        end
        push!(nope, ex.args[1]) # new Ex3 = ... cannot have this on RHS
        ex
    elseif ex.head in (:(=), :(+=)) # easier case of A = stuff
        push!(nope, ex.args[1])
        ex.args[2] = csewalk(ex.args[2], dict, defs, nope)
        ex

    # Then we work on sub-expressions, replace those we're seen immediately,
    # and don't look inside A[i,j] at all:
    elseif haskey(dict, ex)
        dict[ex]
    elseif ex.head == :ref
        ex

    # Simplest case is the last one, replace a whole expression with Ex5 & work inwards.
    # Can't replace "illegal" expressions, but can look for parts which are safe:
    elseif illegal(ex, nope)
        args = Any[x in nope ? x : csewalk(x, dict, defs, nope) for x in ex.args]
        Expr(ex.head, args...)

    elseif ex.head == :call && ex.args[1] in (:*, :+) && length(ex.args) >= 4 # e.g. 1*2*3
        inner = []
        while length(ex.args) >= 3
            pushfirst!(inner, pop!(ex.args))
        end
        binary = Expr(:call, ex.args..., Expr(:call, ex.args[1], inner...))
        csewalk(binary, dict, defs, nope)

    else
        args = Any[csewalk(x, dict, defs, nope) for x in ex.args]
        sy = Symbol(EXPR, length(defs)+1)
        dict[ex] = sy
        # add defn for the outermost operation:
        push!(defs, Expr(:(=), sy, Expr(ex.head, args...)))
        # and return the name for caller:
        sy
    end

illegal(ex, nope) = ex in nope
illegal(ex::Expr, nope) = ex in nope || any(illegal(x, nope) for x in ex.args)

#========== examination ==========#

"""
    Tullio.@printgrad log(x/y) x y

Prints the symbolic gradient, showing `∂f/∂x` and `∂f/∂y` for `f=log(x/y)`.
Useful to check that simplifications, and common subexpression elimination,
are working OK for a given RHS.
"""
macro printgrad(exs...)
    printgrad(exs...)
end

function printgrad(ex::Expr, ts::Symbol...)
    out = quote end
    for t in ts
        df = leibnitz(ex, t)
        dt = Symbol(:δ, t) # Symbol("∂f_∂", t)
        push!(out.args, :($dt = $df))
    end
    print("Initial:\n   ")
    println(join(filter(x -> !(x isa LineNumberNode), out.args), "\n   "))
    print("After CSE:\n   ")
    done = filter(x -> !(x isa LineNumberNode), commonsubex(out).args)
    println(join(done, "\n   "))
    nothing
end

#=

using Tullio: @printgrad

@printgrad  x * y * z   x y z
@printgrad  x * (y * z)   x y z
@printgrad  x + y * z   x y z

@printgrad  1/x   x
@printgrad  x^-1   x   # could make inv(x) for CSE
@printgrad  inv(x)   x
@printgrad  sqrt(x)   x
@printgrad  1/sqrt(x)   x
@printgrad  inv(sqrt(x))   x
@printgrad  x/sqrt(y)   x y

@printgrad  sqrt(x*y)   x y
@printgrad  sqrt(x) * sqrt(y)   x y # worse than line above

@printgrad  1/sqrt(x*y)   x y       # could use repeated CSE

@printgrad  x/sqrt(y*z)   x y z     # could use repeated CSE
@printgrad  x/(sqrt(y)*sqrt(z))   x y z
@printgrad  x*inv(sqrt(y))*inv(sqrt(z))   x y z

@printgrad  x/2   x
@printgrad  x/y   x y

@printgrad  x^2   x
@printgrad  (x*y)^2   x y
@printgrad  (x+y)^3   x y
@printgrad  x^y   x y
@printgrad  log(x)^2   x

@printgrad  log(x)   x
@printgrad  log(x/2)   x
@printgrad  log(2x)   x
@printgrad  log(k*x)   x

@printgrad  x*log(y)   x y

@printgrad  log(x*y)   x y
@printgrad  log(x) + log(y)   x y  # better, now used for log(x*y)

@printgrad  log(x/y)   x y
@printgrad  log(x*inv(y))   x y
@printgrad  log(x)-log(y)   x y    # much better, now used for log(x/y)

@printgrad  log(x/y) * z   x y z
@printgrad  (log(x) - log(y)) * z   x y z
@printgrad  log(x)*z - log(y)* z   x y z

@printgrad  exp(2x)   x
@printgrad  exp(x/y)   x y
@printgrad  exp((x-y)^2/2)   x y

@printgrad  exp(x) * y   x y
@printgrad  exp(x) / 2y   x y

@printgrad a * b / sqrt(d * e)  a b d e
@printgrad x * z / sqrt(y * z)  x y z

=#

#========== the end ==========#