text:
using XFields
using FFTransforms
using Spectra

using LinearAlgebra
using PyPlot

# Define the transform
# ------------------------------------------

trn = let
    𝕨 = 𝕎(Float32, 256, π) ⊗ 𝕎(256, 4.0)
    ordinary_scale(𝕨)*𝕨
end;

# Signal and noise spectral density
# --------------------------------------------

# Here we use Matern spectral densities

Cf = let trn = trn, ρ = 0.15, ν = 2.1, σ² = 1 / 2
    l  = wavenum(trn)
    d  = ndims(l)
    cl = σ² .* Spectra.matern_spec.(l; rho=ρ, nu=ν, d=d)
    Cf = DiagOp(Xfourier(trn, cl))
    @show pixel_variance = (Cf[:] ./ (2π)^(d/2))[1]
    (σ² / pixel_variance) * Cf
end;

# Noise spectral density

Cn =  let trn=trn, μKarcminT=15, ℓknee=4*minimum(Δfreq(trn)), αknee=2

    Ωx_unit = deg2rad(1/60)^2 ## area [rad²] for 1arcmin×1arcmin pixel
    wvn     = wavenum(trn); wvn[1] = Inf
    knee    = @. 1 + XFields.nan2zero((ℓknee / wvn) ^ αknee)
    cnl     = μKarcminT^2 .* Ωx_unit .* knee
    Cn      =  Xfourier(trn, μKarcminT^2 .* Ωx_unit .* knee) |> DiagOp

    Cn
end;

# Mask and Transfer function linear operators
# --------------------------------------------

# Mask

Ma =  let trn=trn, x1bdry = (0.1, 0.9), x2bdry = (0.2,0.95)

    lbr1, rbr1 = trn.period[1] .* x1bdry
    lbr2, rbr2 = trn.period[2] .* x2bdry
    x1, x2 = pix(trn)
    ma = (lbr1 .< x1 .< rbr1) .* (lbr2 .< x2 .< rbr2)'
    Ma = Xmap(trn, ma) |> DiagOp

    Ma
end;

# Transfer function

Tr =  let trn=trn, beam_npix = 4

    fwhm_rad = beam_npix * min(Δpix(trn)...)
    beam = l -> exp(-abs2(l * fwhm_rad) / (16*log(2)))
    tr   = beam.(wavenum(trn))
    tr .*= wavenum(trn) .< 0.9nyq(trn)[1]
    Tr   = Xfourier(trn, tr) |> DiagOp

    Tr
end;

# White noise simulator
# ------------------------------------

function ωη(trn::T) where T<:Transform
    zx = randn(eltype_in(trn),size_in(trn))
    Xmap(trn, zx ./ √Ωpix(trn))
end

# Field simulation: signal (`fsim`), noise (`nsim`) and data (`dsim`)
# --------------------------------------------------------------

dsim, fsim, nsim = let trn=trn, Cf=Cf, Cn=Cn, Cf=Cf, Ma=Ma, Tr=Tr

    fsim = √Cf * ωη(trn)
    nsim = √Cn * ωη(trn)

    dsim = Ma * Tr * fsim + Ma * nsim

    ## μsim = Xmap(trn,sum(dsim[:]) ./ sum(Ma[:]))
    ## dsim = dsim - Ma * μsim

    dsim, fsim, nsim
end;

# Plots of the signal, noise and data
# --------------------------------------------------------------

# ## Signal `fsim`
let trn=trn, f=fsim
    fig, ax = subplots(1, figsize=(8,4))
    x1, x2 = pix(trn)
    pcm = ax.pcolormesh(x2, x1, f[:])
    ax.set_title("signal")
    fig.colorbar(pcm, ax = ax)
    fig.tight_layout()
end;

# ## Noise `nsim`
let trn=trn, f=nsim
    fig, ax = subplots(1, figsize=(8,4))
    x1, x2 = pix(trn)
    pcm = ax.pcolormesh(x2, x1, f[:])
    ax.set_title("noise")
    fig.colorbar(pcm, ax = ax)
    fig.tight_layout()
end;

# ## Data `dsim`
let trn=trn, f=dsim
    fig, ax = subplots(1, figsize=(8,4))
    x1, x2 = pix(trn)
    pcm = ax.pcolormesh(x2, x1, f[:])
    ax.set_title("Data")
    fig.colorbar(pcm, ax = ax)
    fig.tight_layout()
end;

# ## Mask and transfer
let trn=trn, Ma=Ma, Tr=Tr
    fig, ax = subplots(2,1, figsize=(8,8))

    x1, x2 = pix(trn)
    pcm1 = ax[1].pcolormesh(x2, x1, Ma[:])
    ax[1].set_title("Pixel mask")

    l1, l2 = freq(trn)
    p2 = sortperm(l2)
    tr = real.(Tr[!])
    pcm2 = ax[2].pcolormesh(l2[p2], l1, tr[:,p2])
    ax[2].set_title("Fourier transfer function")

    fig.colorbar(pcm1, ax=ax[1])
    fig.colorbar(pcm2, ax=ax[2])
    fig.tight_layout()
end;

# Bandpowers (i.e. periodogram)
# --------------------------------------------------------------

function power(f::Xfield{F}, g::Xfield{F}; bin::Int=2, kmax=Inf, mult=1) where F<:Transform
    trn     = fieldtransform(f)
    k      = wavenum(trn)
    pwr    = @. mult * real(f[!] * conj(g[!]) + conj(f[!]) * g[!]) / 2
    Δbin   = bin * minimum(Δfreq(trn))
    k_left = 0
    while k_left < min(kmax, maximum(k))
        k_right    = k_left + Δbin
        indx       = findall(k_left .< k .<= k_right)
        pwr[indx] .= sum( pwr[indx] ) / length( pwr[indx] )
        k_left  = k_right
    end
    return pwr
end

#-

function power(f::Xfield{F}; bin::Int=2, kmax=Inf, mult=1) where F<:Transform
    power(f, f; bin=bin, kmax=kmax, mult=mult)
end

# ## Noise and signal bandpowers
let trn=trn, Cn=Cn, Cf=Cf, f=fsim, n=nsim
    l     = wavenum(trn)
    mult  = l .* (l .+ 1)

    fig, ax = subplots(1, figsize=(8,4))

    pwrf = power(f; mult=mult * Ωfreq(trn) )
    pwrn = power(n; mult=mult * Ωfreq(trn) )
    (l[:,1], pwrf[:,1]) |> x->ax.plot(x[1][2:end],x[2][2:end], label="signal")
    (l[:,1], pwrn[:,1]) |> x->ax.plot(x[1][2:end],x[2][2:end], label="noise")

    cf = real.(Cf[!])
    cn = real.(Cn[!])
    (l[:,1], (mult.*cf)[:,1]) |> x->ax.plot(x[1][3:end],x[2][3:end])
    (l[:,1], (mult.*cn)[:,1]) |> x->ax.plot(x[1][3:end],x[2][3:end])

    ax.set_xlabel("wavenumber")
    ax.set_ylabel("power")
    ax.legend()
    fig.tight_layout()
end;

# Set up basic  d = M⋅Tf⋅f + M⋅n
# --------------------------------------------------------------

# custom pcg with function composition (Minv * A \approx I)
function pcg(Minv::Function, A::Function, b, x=0*b; nsteps::Int=75, rel_tol::Float64 = 1e-8)
    r       = b - A(x)
    z       = Minv(r)
    p       = deepcopy(z)
    res     = dot(r,z)
    reshist = Vector{typeof(res)}()
    for i = 1:nsteps
        Ap        = A(p)
        α         = res / dot(p,Ap)
        x         = x + α * p
        r         = r - α * Ap
        z         = Minv(r)
        res′      = dot(r,z)
        p         = z + (res′ / res) * p
        rel_error = XFields.nan2zero(sqrt(dot(r,r)/dot(b,b)))
        if rel_error < rel_tol
            return x, reshist
        end
        push!(reshist, rel_error)
        res = res′
    end
    return x, reshist
end

#-

function LinearAlgebra.dot(f::Xfield{FT},g::Xfield{FT}) where FT<:Transform
    trn = fieldtransform(f)
    Ωpix(trn) * dot(f[:],g[:])
end

#-

wfsim, wfhist, zwf = let trn=trn, Cn=Cn, Cf=Cf, Tr=Tr, Ma=Ma, dsim=dsim

    A  = Ma * Tr / Cn * Tr * Ma
    B  = 1 / Cf
    P  = inv(Tr / Cn * Tr + B)

    wfsim, wfhist = pcg(
            w -> P * w,
            w -> A * w +  B * w,
            Ma * Tr / Cn * dsim,
            nsteps  = 500,
            rel_tol = 1e-15,
    )

    ## compte the chi2 z-score
    dfd   = sum(Ma[:]) # degrees of freedom of the data
    Δ1    = dsim - Ma * wfsim
    zwf  = - dot(Δ1, Cn \ Δ1) - dot(wfsim, Cf \ wfsim)
    zwf -= -dfd
    zwf /= sqrt(2*dfd)

    wfsim, wfhist, zwf

end;

# the "residual" per iteration

let wfhist=wfhist
    fig, ax = subplots(1, figsize=(8,4))
    semilogy(wfhist)
end;

# The Wiener filter

let trn=trn, f=wfsim
    fig, ax = subplots(1, figsize=(8,4))
    x1, x2 = pix(trn)
    vm = extrema(f[:]) .|> abs |> x->max(x...)
    pcm = ax.pcolormesh(x2, x1, f[:],vmin=-vm, vmax=vm)
    ax.set_title("Wiener filter")
    fig.colorbar(pcm, ax = ax)
    fig.tight_layout()
end;


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rocm_cmake_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rocm_cmake")
JLLWrappers.@generate_main_file("rocm_cmake", UUID("c3eb98f6-eb42-520d-ab3c-506c3b64aac0"))
end  # module rocm_cmake_jll

ext: jl


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Stage-decomposition generation #
# ============================== #
function generate!(stochasticprogram::StochasticProgram{N}, structure::StageDecompositionStructure{N}) where N
    # Generate all stages
    for stage in 1:N
        generate!(stochasticprogram, structure, stage)
    end
    return nothing
end

function generate!(stochasticprogram::StochasticProgram{N}, structure::StageDecompositionStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    if stage == 1
        # Check generators
        has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
        # Set the optimizer (if any)
        if has_provided_optimizer(stochasticprogram.optimizer)
            master_opt = master_optimizer(stochasticprogram)
            if master_opt != nothing
                set_optimizer(structure.first_stage, master_opt)
            end
        end
        # Prepare decisions
        structure.first_stage.ext[:decisions] = Decisions((structure.decisions[1],))
        add_decision_bridges!(structure.first_stage)
        # Generate first stage
        generator(stochasticprogram, :stage_1)(structure.first_stage, stage_parameters(stochasticprogram, 1))
    else
        # Check generators
        stage_key = Symbol(:stage_, stage)
        decision_key = Symbol(:stage_, stage - 1, :_decisions)
        has_generator(stochasticprogram, stage_key) || error("Stage problem $stage not defined in stochastic program. Consider @stage $stage.")
        has_generator(stochasticprogram, decision_key) || error("No decision variables defined in stage problem $(stage-1).")
        # Sanity check on scenario probabilities
        if num_scenarios(stochasticprogram, stage) > 0
            p = stage_probability(stochasticprogram, stage)
            abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
        end
        # Generate
        generate_vertical!(scenarioproblems(structure, stage),
                           stage,
                           generator(stochasticprogram, decision_key),
                           generator(stochasticprogram, stage_key),
                           stage_parameters(stochasticprogram, stage - 1),
                           stage_parameters(stochasticprogram, stage),
                           structure.decisions[stage],
                           subproblem_optimizer(stochasticprogram))
    end
    return nothing
end

function generate_vertical!(scenarioproblems::ScenarioProblems,
                            stage::Integer,
                            decision_generator::Function,
                            generator::Function,
                            decision_params::Any,
                            stage_params::Any,
                            decision_map::DecisionMap,
                            optimizer)
    for i in num_subproblems(scenarioproblems)+1:num_scenarios(scenarioproblems)
        # Create subproblem
        subproblem = optimizer == nothing ? Model() : Model(optimizer)
        # Prepare decisions
        decisions = ntuple(Val{stage}()) do s
            if s == stage - 1
                # Known decisions from the previous stages are
                # the same everywhere.
                decision_map
            else
                # Remaining decisions are unique to each subproblem
                DecisionMap()
            end
        end
        subproblem.ext[:decisions] = Decisions(decisions; is_node = true)
        add_decision_bridges!(subproblem)
        # Generate and return the stage model
        decision_generator(subproblem, decision_params)
        generator(subproblem, stage_params, scenario(scenarioproblems, i))
        push!(scenarioproblems.problems, subproblem)
    end
    return nothing
end
function generate_vertical!(scenarioproblems::DistributedScenarioProblems,
                            stage::Integer,
                            decision_generator::Function,
                            generator::Function,
                            decision_params::Any,
                            stage_params::Any,
                            ::DecisionMap,
                            optimizer)
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                stage,
                decision_generator,
                generator,
                decision_params,
                stage_params,
                scenarioproblems.decisions[w-1],
                optimizer) do sp, stage, dgenerator, generator, dparams, params, decisions, opt
                    generate_vertical!(fetch(sp),
                                       stage,
                                       dgenerator,
                                       generator,
                                       dparams,
                                       params,
                                       fetch(decisions),
                                       opt)
                end
        end
    end
    return nothing
end

function clear!(structure::StageDecompositionStructure{N}) where N
    # Clear decisions
    clear!(structure.decisions)
    # Clear all stages
    for stage in 1:N
        clear_stage!(structure, stage)
    end
    return nothing
end

function clear_stage!(structure::StageDecompositionStructure{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    if s == 1
        empty!(structure.first_stage)
    else
        clear!(scenarioproblems(structure, s))
    end
    return nothing
end

# Getters #
# ========================== #
function first_stage(stochasticprogram::StochasticProgram, structure::StageDecompositionStructure; optimizer = nothing)
    if optimizer == nothing
        return structure.first_stage
    end
    stage_one = copy(structure.first_stage)
    set_optimizer(stage_one, optimizer)
    return stage_one
end


text:
# Julia wrapper for header: /usr/include/scip/cons_orbitope.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrOrbitope(scip)
    ccall((:SCIPincludeConshdlrOrbitope, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsOrbitope(scip, cons, name, vars, orbitopetype, nspcons, nblocks, resolveprop, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsOrbitope, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{Ptr{SCIP_VAR}}}, SCIP_ORBITOPETYPE, Cint, Cint, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, vars, orbitopetype, nspcons, nblocks, resolveprop, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicOrbitope(scip, cons, name, vars, orbitopetype, nspcons, nblocks, resolveprop)
    ccall((:SCIPcreateConsBasicOrbitope, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{Ptr{SCIP_VAR}}}, SCIP_ORBITOPETYPE, Cint, Cint, UInt32), scip, cons, name, vars, orbitopetype, nspcons, nblocks, resolveprop)
end


text:
function op_K( pw::PWGrid, psi::Array{ComplexF64,2} )
    out = zeros(ComplexF64,size(psi))
    Ncol = size(psi,2)
    Ω  = pw.Ω
    G2 = pw.G2
    Npoints = pw.Npoints
    for is = 1:Ncol
        for ip = 1:Npoints
            out[ip,is] = psi[ip,is]*G2[ip]
        end
    end
    return 0.5*out # two minus signs
end


text:
using Base.Test
using Zlib
using Compat

data = convert(Vector{UInt8}, rand(1:255, 1000000))
decompressed = decompress(compress(data))
@test data == decompressed

decompressed = decompress(compress(data, false, true), true)
@test data == decompressed

@test length(decompress(compress(""))) == 0

b = IOBuffer()
w = Zlib.Writer(b)
n = 0
while n < length(data)
    n += write(w, data[n+1:n+200000])
end
close(w)
seekstart(b)
@test data == decompress(readbytes(b))

seekstart(b)
r = Zlib.Reader(b)
decompressed = Array(UInt8, 0)
while !eof(r)
    append!(decompressed, read(r, UInt8, 200000))
end
@test data == decompressed

data = Any[
    convert(UInt8, 20),
    42,
    float(3.14),
    "julia",
    rand(5),
    rand(3, 4),
    sub(rand(10,10), 2:8,2:4)
]

b = IOBuffer()
@test_throws ErrorException read(w, UInt8, 1)
w = Zlib.Writer(b)
for x in data
    write(w, x)
end
close(w)

seekstart(b)
r = Zlib.Reader(b)
@test_throws ErrorException write(r, convert(UInt8, 20))
for x in data
    if typeof(x) == ASCIIString
        @test x == ASCIIString(read(r, UInt8, length(x)))
    elseif typeof(x) <: Array
        y = similar(x)
        y[:] = 0
        @test x == read!(r, y)
        @test x == y
    elseif typeof(x) <: SubArray
        continue # Base knows how to write, but not read
    else
        @test x == read(r, typeof(x))
    end
end
close(r)

@test crc32("") == 0
@test crc32("hello") == 0x3610a686
@test crc32("Julia programming language") == 0xfc485364
crc = crc32("Julia programming")
@test crc == 0xc7db4271
@test crc32(" language", crc) == 0xfc485364

@test_throws ErrorException decompress(compress("abcdefghijklmnopqrstuvwxyz")[1:10])


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    # rootdet
    ((2, 2, MatNegGeom()),),
    ((2, 2, MatNegGeomEFExp()),),
    ((2, 2, MatNegGeomEFPow()),),
    # tr inv
    ((2, 2, MatInv()),),
    ((2, 2, MatInvEigOrd()),),
    ((2, 2, MatInvDirect()),),
    # tr neg2sqrt
    ((1, 2, MatNeg2Sqrt()),),
    ((1, 2, MatNeg2SqrtEigOrd()),),
    # tr negexp1 (domain not positive)
    ((1, 2, MatNegExp1()),),
    ((1, 2, MatNegExp1EigOrd()),),
    # tr power12conj (domain not positive)
    ((1, 2, MatPower12Conj(1.5)),),
    ((1, 2, MatPower12ConjEigOrd(1.5)),),
    ]
insts["fast"] = [
    ((1, 10, MatNegGeomEFExp()),),
    ((1, 15, MatNegGeom()),),
    ((2, 3, MatNegGeomEFPow()),),
    ((2, 3, MatInvEigOrd()),),
    ((2, 5, MatInvDirect()),),
    ((2, 6, MatInv()),),
    ((2, 6, MatNegEntropy()),),
    ((2, 7, MatNegLog()),),
    ((3, 2, MatNegEntropyEigOrd()),),
    ((3, 2, MatNegExp1()),),
    ((3, 2, MatNegExp1EigOrd()),),
    ((3, 2, MatPower12Conj(1.5)),),
    ((3, 2, MatPower12ConjEigOrd(1.5)),),
    ((3, 2, MatPower12EigOrd(1.5)),),
    ((3, 4, MatPower12(1.5)),),
    ((3, 4, MatNegGeom()),),
    ((7, 2, MatNegLog()),),
    ((7, 2, MatInv()),),
    ((7, 2, MatNegEntropy()),),
    ((7, 2, MatPower12(1.5)),),
    ]
insts["various"] = [
    ((2, 5, MatInvEigOrd()),),
    ((2, 5, MatNegLogEigOrd()),),
    ((2, 4, MatPower12EigOrd(1.5)), nothing, relaxed_tols),
    ((2, 4, MatNegEntropyEigOrd()),),
    ((2, 10, MatNegGeomEFExp()),),
    ((2, 8, MatNegGeomEFPow()),),
    ((2, 14, MatNegGeom()),),
    ((2, 12, MatInv()),),
    ((2, 10, MatNegLog()),),
    ((2, 8, MatNegEntropy()),),
    ((2, 6, MatPower12(1.5)),),
    ((3, 6, MatInvDirect()),),
    ((3, 5, MatNegLogDirect()),),
    ((3, 6, MatNegEntropy()),),
    ((3, 6, MatNegLog()),),
    ((3, 5, MatPower12(1.5)),),
    ((3, 5, MatNegExp1()),),
    ((3, 4, MatPower12Conj(1.5)),),
    ((5, 2, MatNegExp1EigOrd()), nothing, relaxed_tols),
    ((5, 2, MatPower12ConjEigOrd(1.5)),),
    ((5, 3, MatInvDirect()),),
    ((5, 3, MatNegLogDirect()),),
    ((6, 3, MatNegGeom()),),
    ((6, 3, MatInv()),),
    ]
return (CentralPolyMatJuMP, insts)


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: BoundingBox, get_overlapping_fields,
                            one_node_joint_infer, infer_init, find_neighbors
import Celeste.SDSSIO: RunCamcolField, load_field_images, PlainFITSStrategy

const rcfs = [
    RunCamcolField(4264,6,160),
    RunCamcolField(4264,6,161),
    RunCamcolField(4264,5,158),
    RunCamcolField(4264,5,159),
    RunCamcolField(4264,5,160),
    RunCamcolField(4264,5,161),
    RunCamcolField(4264,5,162),
    RunCamcolField(4294,6,133),
    RunCamcolField(4294,6,134),
    RunCamcolField(4294,6,135),
    RunCamcolField(4294,6,136),
    RunCamcolField(4294,6,137),
    RunCamcolField(4294,6,138),
    RunCamcolField(4294,5,135)]

const datadir = joinpath(dirname(@__FILE__), "..", "..", "test", "data")
wd = pwd()
cd(datadir)
for rcf in rcfs
    run(`make RUN=$(rcf.run) CAMCOL=$(rcf.camcol) FIELD=$(rcf.field)`)
end
cd(wd)

function benchmark_quarter_degree()
    box = BoundingBox(124.0, 124.5, 58.5, 59.0)
    rcfs = get_overlapping_fields(box, datadir)

    # ctni = (catalogs, target, neighbor_map, images)
    ctni = infer_init(rcfs, PlainFITSStrategy(datadir); box=box)[1:4]

    # Warm up---this compiles the code
    ctni2 = (ctni[1], ctni[2][1:1], ctni[3][1:1], ctni[4][1:1])
    @time one_node_joint_infer(ctni2...)
    println("Done with warm up")

    # resets runtime profiler *and* count for --track-allocation
    Profile.clear_malloc_data()

    if isempty(ARGS)
        @time  one_node_joint_infer(ctni...)
    elseif ARGS[1] == "--profile"
        Profile.init(delay=1.0)
        @profile one_node_joint_infer(ctni...)
        Profile.print(format=:flat, sortedby=:count)
    end
end

benchmark_quarter_degree()


text:
using SourceWalk: textmap, textwalk, @capture, isexpr
using Test

testrep(ex, text) = textmap(_ -> Expr(:file, ex), text)

# Replacement

@test testrep(:(b*2), "a*2") == "b*2"
@test testrep(:(b*2), "a * 2") == "b * 2"
@test testrep(:(2*b*c), "a * b * c") == "2 * b * c"
@test testrep(:(f(a, b)), "a * b") == "f(a, b)"

@test testrep(:((a+1)*2), "a * 2") == "(a + 1) * 2"
@test testrep(:((a+1)*2), "*(a, 2)") == "*(a + 1, 2)"

@test testrep(:(bar(uppercase("baz"))), """
  bar("baz")
  """) == """
  bar(uppercase("baz"))
  """

@test testrep(:(using Bar), "using Foo") == "using Bar"

# Insertion

@test testrep(:(f(a,b)), "f(a)") == "f(a, b)"
@test testrep(:(f(a)), "f()") == "f(a)"

@test testrep(:(f(a, b, c)), "f(a, b)") == "f(a, b, c)"
@test testrep(:(f(a, c, b)), "f(a, b)") == "f(a, c, b)"
@test testrep(:(f(c, a, b)), "f(a, b)") == "f(c, a, b)"
@test testrep(:(c(f, a, b)), "f(a, b)") == "c(f, a, b)"

@test testrep(:(f(a, b, c)), "f(a , b)") == "f(a , b , c)"
@test testrep(:(a+b+c), "a + b") == "a + b + c"
@test testrep(:(a+b+c), "a+b") == "a+b+c"

testrep(:(1; 2; 3), """
  begin
    1
    2
  end
  """) == """
  begin
    1
    2
    3
  end
  """

testrep(:(1; 2; 3), """
  begin
    1; 2
  end
  """) == """
  begin
    1; 2; 3
  end
  """

testrep(:(1; 2; 3), """
  begin
    1; 2;
  end
  """) == """
  begin
    1; 2; 3;
  end
  """

testrep(:(1; 2; 3), "(1; 2)") == "(1; 2; 3)"

# Deletion

@test testrep(:(f()), "f(a)") == "f()"
@test testrep(:(f()), "f(a,)") == "f()"

@test testrep(:(f(a, b)), "f(a, b, c)") == "f(a, b)"
@test testrep(:(f(a, c)), "f(a, b, c)") == "f(a, c)"
@test testrep(:(f(b, c)), "f(a, b, c)") == "f(b, c)"
@test testrep(:(a(b, c)), "f(a, b, c)") == "a(b, c)"

# Others

@test testrep(:(@foo a), "@foo a b") == "@foo a"
@test testrep(:(@foo a b c), "@foo a b") == "@foo a b c"
@test testrep(:(@foo a b c), "@foo(a, b)") == "@foo(a, b, c)"

@test testrep(:(a,), "(a, b)") == "(a)"
@test testrep(:(a,b,c), "a, b") == "a, b, c"


text:
using SliceMap
using Test
using ForwardDiff, Tracker, Zygote, TensorCast, JuliennedArrays

Zygote.refresh()

@testset "columns" begin

    mat = rand(1:9, 3,10)
    fun(x) = 2 .+ x.^2
    res = mapslices(fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)
    @test res ≈ MapCols(fun, mat)

    @test res ≈ tmapcols(fun, mat)
    @test res ≈ ThreadMapCols{3}(fun, mat)
    @test res ≈ ThreadMapCols(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols(fun, m)), mat)[1]

    @test grad ≈ Tracker.gradient(m -> sum(sin, tmapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, ThreadMapCols{3}(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, ThreadMapCols(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, tmapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, ThreadMapCols{3}(fun, m)), mat)[1]

    tcm(mat) = @cast out[i,j] := fun(mat[:,j])[i]
    @test res ≈ tcm(mat)
    @test grad ≈ Zygote.gradient(m -> sum(sin, tcm(m)), mat)[1]

    jcols(f,m) = Align(map(f, Slices(m, True(), False())), True(), False())
    @test res ≈ jcols(fun, mat)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jcols(fun, m)), mat)[1]

end
@testset "columns -> scalar" begin

    mat = rand(1:9, 3,10)
    fun(x) = sum(x) # different function!
    res = mapslices(fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    # tcm3(mat) = @cast out[_,j] := fun(mat[:,j]) # changed here too
    # @test res ≈ tcm3(mat)
    # @test grad ≈ Zygote.gradient(m -> sum(sin, tcm3(m)), mat)[1]

end
@testset "columns -> matrix" begin

    mat = rand(1:9, 3,10)
    fun(x) = x .* x' # different function! vector -> matrix
    res = mapslices(vec∘fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(vec∘fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    # tcm4(mat) = @cast out[i⊗i′,j] := fun(mat[:,j])[i,i′]  i:3
    # @test res ≈ tcm4(mat)
    # @test grad ≈ Zygote.gradient(m -> sum(sin, tcm4(m)), mat)[1]

end
@testset "columns w args" begin

    mat = randn(Float32, 3,10)
    fun(x, s) = 1 .+ x .* s
    res = mapslices(x -> vec(fun(x,5)), mat, dims=1)

    @test res ≈ mapcols(fun, mat, 5)
    @test res ≈ MapCols{3}(fun, mat, 5)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(x -> vec(fun(x,5)), m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m, 5)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m, 5)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m, 5)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m, 5)), mat)[1]

    # tcm5(mat) = @cast out[i,j] := fun(mat[:,j], 5)[i]
    # @test res ≈ tcm5(mat)
    # @test grad ≈ Zygote.gradient(m -> sum(sin, tcm5(m)), mat)[1]

end
@testset "rows" begin

    mat = randn(4,5)
    fun(x) = 2 .+ x.^2 ./ sum(x)

    res = mapslices(fun, mat, dims=2)
    @test res ≈ maprows(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=2)), mat)
    @test grad ≈ Tracker.gradient(m -> sum(sin, maprows(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, maprows(fun, m)), mat)[1]

    # tcm2(mat) = @cast out[i,j] := fun(mat[i,:])[j]
    # @test res ≈ tcm2(mat)
    # @test grad ≈ Zygote.gradient(m -> sum(sin, tcm2(m)), mat)[1]

    jrows(f,m) = Align(map(f, Slices(m, False(), True())), False(), True())
    @test res ≈ jrows(fun, mat)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jrows(fun, m)), mat)[1]

end
@testset "slices of a 4-tensor" begin

    ten = randn(3,4,5,2)
    fun(x::AbstractVector) = sqrt(3) .+ x.^3 ./ (sum(x)^2)
    res = mapslices(fun, ten, dims=3)

    @test res ≈ slicemap(fun, ten, dims=3)

    grad = ForwardDiff.gradient(x -> sum(sin, slicemap(fun, x, dims=3)), ten)
    @test grad ≈ Zygote.gradient(x -> sum(sin, slicemap(fun, x, dims=3)), ten)[1]

    jthree(f,m) = Align(map(f,
        Slices(m, False(), False(), True(), False())
        ), False(), False(), True(), False())
    @test res ≈ jthree(fun, ten)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jthree(fun, m)), ten)[1]

    j3(f,m) = Align(map(f, Slices(m, 3)), 3)
    @test res ≈ j3(fun, ten)
    @test grad ≈ Zygote.gradient(m -> sum(sin, j3(fun, m)), ten)[1]

end


text:
"""
    CpG_windows_expected(positions::Array, refseq::NamedTuple; window=5, step=5, noslide=false)

Ouput the expected number of CpG sites for the specified genomic positions using a user-defined window size and step. Positions are an array of 1:123-345 style positions. A genomic feature may be skipped by inputing a 2d array of postions. For example, 1:100-200  1:230-300, skips the feature from 1:201-229.

The refseq input is generated using the readfasta function.
"""
function CpG_windows_expected(positions::Array, refseq::NamedTuple; window=5, step=5, noslide=false)

    ref = refseq
    m = []
    xvals = []

    for i in 1:size(positions,1)

        if size(positions,2) == 2
            s1 = getref(positions[i,1], ref)
            s2 = getref(positions[i,2], ref)
            s = s1 * s2
        else
            s = getref(positions[i], ref)
        end

        if noslide == true
            window = length(s)
            step = length(s)
        end

        gcdata = sliding_window(s, window, step, expected_CpG)

        push!(m, gcdata)

    end

    return m

end


text:
using ACEexperimental
using Test

@testset "ACEexperimental.jl" begin

include("test_combinations.jl")

end


text:
# simd version is defined in simd.jl
@generated function _map(f, xs::Vararg{AbstractTensor, N}) where {N}
    S = promote_space(map(Space, xs)...)
    exps = map(indices(S)) do i
        vals = [:(xs[$j][$i]) for j in 1:N]
        :(f($(vals...)))
    end
    TT = tensortype(S)
    return quote
        @_inline_meta
        @inbounds $TT($(exps...))
    end
end

@inline Base.:+(x::AbstractTensor, y::AbstractTensor) = _map(+, x, y)
@inline Base.:-(x::AbstractTensor, y::AbstractTensor) = _map(-, x, y)
@inline Base.:*(y::Real, x::AbstractTensor) = _map(x -> x*y, x)
@inline Base.:*(x::AbstractTensor, y::Real) = _map(x -> x*y, x)
@inline Base.:/(x::AbstractTensor, y::Real) = _map(x -> x/y, x)
@inline Base.:+(x::AbstractTensor) = x
@inline Base.:-(x::AbstractTensor) = _map(-, x)

# with AbstractArray
@generated Base.:+(x::AbstractTensor, y::AbstractArray) = :(@_inline_meta; x + convert($(tensortype(Space(size(x)))), y))
@generated Base.:+(x::AbstractArray, y::AbstractTensor) = :(@_inline_meta; convert($(tensortype(Space(size(y)))), x) + y)
@generated Base.:-(x::AbstractTensor, y::AbstractArray) = :(@_inline_meta; x - convert($(tensortype(Space(size(x)))), y))
@generated Base.:-(x::AbstractArray, y::AbstractTensor) = :(@_inline_meta; convert($(tensortype(Space(size(y)))), x) - y)

@inline _add_uniform(x::AbstractSquareTensor, λ::Real) = x + λ * one(x)
@inline Base.:+(x::AbstractTensor, y::UniformScaling) = _add_uniform( x,  y.λ)
@inline Base.:-(x::AbstractTensor, y::UniformScaling) = _add_uniform( x, -y.λ)
@inline Base.:+(x::UniformScaling, y::AbstractTensor) = _add_uniform( y,  x.λ)
@inline Base.:-(x::UniformScaling, y::AbstractTensor) = _add_uniform(-y,  x.λ)
@inline dot(x::Union{AbstractVec, AbstractMat, AbstractSquareTensor}, y::UniformScaling) = x * y.λ
@inline dot(x::UniformScaling, y::Union{AbstractVec, AbstractMat, AbstractSquareTensor}) = x.λ * y
@inline double_contraction(x::AbstractSquareTensor, y::UniformScaling) = tr(x) * y.λ
@inline double_contraction(x::UniformScaling, y::AbstractSquareTensor) = x.λ * tr(y)

# error for standard multiplications
error_multiply() = error("use `⋅` (`\\cdot`) for single contraction and `⊡` (`\\boxdot`) for double contraction instead of `*`")
Base.:*(::AbstractTensor, ::AbstractTensor) = error_multiply()
Base.:*(::AbstractTensor, ::UniformScaling) = error_multiply()
Base.:*(::UniformScaling, ::AbstractTensor) = error_multiply()

function contraction_exprs(x::Type{<: AbstractTensor}, y::Type{<: AbstractTensor}, N::Int)
    nx = ndims(x)
    ny = ndims(y)
    ij_x = UnitRange(1, nx)
    ij_y = UnitRange(nx + 1 - N, nx + ny - N)
    freeinds = find_freeindices([ij_x; ij_y])
    einsum_contraction_expr(Tuple(freeinds), (x, y), (ij_x, ij_y))
end

"""
    contraction(::AbstractTensor, ::AbstractTensor, ::Val{N})

Conduct contraction of `N` inner indices.
For example, `N=2` contraction for third-order tensors ``A_{ij} = B_{ikl} C_{klj}``
can be computed in Tensorial.jl as

```jldoctest
julia> B = rand(Tensor{Tuple{3,3,3}});

julia> C = rand(Tensor{Tuple{3,3,3}});

julia> A = contraction(B, C, Val(2))
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 1.36912   1.86751  1.32531
 1.61744   2.34426  1.94101
 0.929252  1.89656  1.79015
```

Following symbols are also available for specific contractions:

- `x ⊗ y` (where `⊗` can be typed by `\\otimes<tab>`): `contraction(x, y, Val(0))`
- `x ⋅ y` (where `⋅` can be typed by `\\cdot<tab>`): `contraction(x, y, Val(1))`
- `x ⊡ y` (where `⊡` can be typed by `\\boxdot<tab>`): `contraction(x, y, Val(2))`
"""
@generated function contraction(x::AbstractTensor, y::AbstractTensor, ::Val{N}) where {N}
    TT, exps = contraction_exprs(x, y, N)
    quote
        @_inline_meta
        tensors = (x, y)
        @inbounds $TT($(exps...))
    end
end

"""
    otimes(x::AbstractTensor, y::AbstractTensor)
    x ⊗ y

Compute tensor product such as ``A_{ij} = x_i y_j``.
`x ⊗ y` (where `⊗` can be typed by `\\otimes<tab>`) is a synonym for `otimes(x, y)`.

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

julia> y = rand(Vec{3})
3-element Vec{3, Float64}:
 0.4600853424625171
 0.7940257103317943
 0.8541465903790502

julia> A = x ⊗ y
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.271839  0.469146  0.504668
 0.352792  0.608857  0.654957
 0.260518  0.449607  0.48365
```
"""
@inline otimes(x1::AbstractTensor, x2::AbstractTensor) = contraction(x1, x2, Val(0))

"""
    dot(x::AbstractTensor, y::AbstractTensor)
    x ⋅ y

Compute dot product such as ``a = x_i y_i``.
This is equivalent to [`contraction(::AbstractTensor, ::AbstractTensor, Val(1))`](@ref).
`x ⋅ y` (where `⋅` can be typed by `\\cdot<tab>`) is a synonym for `dot(x, y)`.

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

julia> y = rand(Vec{3})
3-element Vec{3, Float64}:
 0.4600853424625171
 0.7940257103317943
 0.8541465903790502

julia> a = x ⋅ y
1.3643452781654772
```
"""
@inline dot(x1::AbstractTensor, x2::AbstractTensor) = contraction(x1, x2, Val(1))
@inline double_contraction(x1::AbstractTensor, x2::AbstractTensor) = contraction(x1, x2, Val(2))

"""
    norm(::AbstractTensor)

Compute norm of a tensor.

# Examples
```jldoctest
julia> x = rand(Mat{3, 3})
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> norm(x)
1.7377443667834922
```
"""
@inline norm(x::AbstractTensor) = sqrt(contraction(x, x, Val(ndims(x))))

@inline normalize(x::AbstractTensor) = x / norm(x)

# v_k * S_ikjl * u_l
@inline function dotdot(v1::Vec{dim}, S::SymmetricFourthOrderTensor{dim}, v2::Vec{dim}) where {dim}
    S′ = SymmetricFourthOrderTensor{dim}((i,j,k,l) -> @inbounds S[j,i,l,k])
    v1 ⋅ S′ ⋅ v2
end

"""
    tr(::AbstractSecondOrderTensor)
    tr(::AbstractSymmetricSecondOrderTensor)

Compute the trace of a square tensor.

# Examples
```jldoctest
julia> x = rand(Mat{3,3})
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> tr(x)
1.6317075356075135
```
"""
@inline tr(x::AbstractSquareTensor) = @einsum x[i,i]

# symmetric
@inline symmetric(x::AbstractSymmetricSecondOrderTensor) = x
@inline function symmetric(x::AbstractSymmetricSecondOrderTensor, uplo::Symbol)
    if uplo == :U || uplo == :L
        return x
    end
    throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
end
## :U or :L
for (uplo, filter) in ((:U, (inds,i,j) -> i ≤ j ? inds[i,j] : inds[j,i]),
                       (:L, (inds,i,j) -> i ≥ j ? inds[i,j] : inds[j,i]))
    @eval @generated function $(Symbol(:symmetric, uplo))(x::AbstractSecondOrderTensor{dim}) where {dim}
        exps = [:(Tuple(x)[$($filter(independent_indices(x), Tuple(I)...))]) for I in CartesianIndices(x)]
        TT = SymmetricSecondOrderTensor{dim}
        quote
            @_inline_meta
            @inbounds $TT($(exps[indices(TT)]...))
        end
    end
end
@inline function symmetric(x::AbstractSecondOrderTensor{dim}, uplo::Symbol) where {dim}
    uplo == :U && return symmetricU(x)
    uplo == :L && return symmetricL(x)
    throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
end
## no uplo
@inline symmetric(x::AbstractSecondOrderTensor{dim}) where {dim} =
    SymmetricSecondOrderTensor{dim}((i,j) -> @inbounds i == j ? x[i,j] : (x[i,j] + x[j,i]) / 2)
@inline symmetric(x::AbstractSecondOrderTensor{1}) =
    @inbounds SymmetricSecondOrderTensor{1}(x[1,1])
@inline symmetric(x::AbstractSecondOrderTensor{2}) =
    @inbounds SymmetricSecondOrderTensor{2}(x[1,1], (x[2,1]+x[1,2])/2, x[2,2])
@inline symmetric(x::AbstractSecondOrderTensor{3}) =
    @inbounds SymmetricSecondOrderTensor{3}(x[1,1], (x[2,1]+x[1,2])/2, (x[3,1]+x[1,3])/2, x[2,2], (x[3,2]+x[2,3])/2, x[3,3])

"""
    skew(::AbstractSecondOrderTensor)
    skew(::AbstractSymmetricSecondOrderTensor)

Compute skew-symmetric (anti-symmetric) part of a second order tensor.
"""
@inline skew(x::AbstractSecondOrderTensor) = (x - x') / 2
@inline skew(x::AbstractSymmetricSecondOrderTensor{dim, T}) where {dim, T} = zero(SecondOrderTensor{dim, T})

"""
    skew(ω::Vec{3})

Construct a skew-symmetric (anti-symmetric) tensor `W` from a vector `ω` as

```math
\\bm{\\omega} = \\begin{Bmatrix}
    \\omega_1 \\\\
    \\omega_2 \\\\
    \\omega_3
\\end{Bmatrix}, \\quad
\\bm{W} = \\begin{bmatrix}
     0         & -\\omega_3 &  \\omega_2 \\\\
     \\omega_3 & 0          & -\\omega_1 \\\\
    -\\omega_2 &  \\omega_1 &  0
\\end{bmatrix}
```

# Examples
```jldoctest
julia> skew(Vec(1,2,3))
3×3 Tensor{Tuple{3, 3}, Int64, 2, 9}:
  0  -3   2
  3   0  -1
 -2   1   0
```
"""
@inline function skew(ω::Vec{3})
    z = zero(eltype(ω))
    @inbounds @Mat [ z    -ω[3]  ω[2]
                     ω[3]     z -ω[1]
                    -ω[2]  ω[1]     z]
end

# transpose/adjoint
@inline transpose(x::AbstractTensor{Tuple{@Symmetry{dim, dim}}}) where {dim} = x
@inline transpose(x::AbstractTensor{Tuple{m, n}}) where {m, n} = Tensor{Tuple{n, m}}((i,j) -> @inbounds x[j,i])
@inline adjoint(x::AbstractTensor) = transpose(x)

# det
@generated function extract_vecs(x::AbstractSquareTensor{dim}) where {dim}
    exps = map(1:dim) do j
        :(Vec($([getindex_expr(x, :x, i, j) for i in 1:dim]...)))
    end
    quote
        @_inline_meta
        @inbounds tuple($(exps...))
    end
end
@inline function det(x::AbstractSquareTensor{1})
    @inbounds x[1,1]
end
@inline function det(x::AbstractSquareTensor{2})
    @inbounds x[1,1] * x[2,2] - x[1,2] * x[2,1]
end
@inline function det(A::AbstractSquareTensor{3})
    a1, a2, a3 = extract_vecs(A)
    (a1 × a2) ⋅ a3
end
@inline function det(x::AbstractSquareTensor{dim}) where {dim}
    det(SMatrix{dim, dim}(x))
end

"""
    cross(x::Vec{3}, y::Vec{3}) -> Vec{3}
    cross(x::Vec{2}, y::Vec{2}) -> Vec{3}
    cross(x::Vec{1}, y::Vec{1}) -> Vec{3}
    x × y

Compute the cross product between two vectors.
The vectors are expanded to 3D frist for dimensions 1 and 2.
The infix operator `×` (written `\\times`) can also be used.
`x × y` (where `×` can be typed by `\\times<tab>`) is a synonym for `cross(x, y)`.

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

julia> y = rand(Vec{3})
3-element Vec{3, Float64}:
 0.4600853424625171
 0.7940257103317943
 0.8541465903790502

julia> x × y
3-element Vec{3, Float64}:
  0.20535000738340053
 -0.24415039787171888
  0.11635375677388776
```
"""
@inline cross(x::Vec{1, T1}, y::Vec{1, T2}) where {T1, T2} = zero(Vec{3, promote_type(T1, T2)})
@inline function cross(x::Vec{2, T1}, y::Vec{2, T2}) where {T1, T2}
    z = zero(promote_type(T1, T2))
    @inbounds Vec(z, z, x[1]*y[2] - x[2]*y[1])
end
@inline function cross(x::Vec{3}, y::Vec{3})
    @inbounds Vec(x[2]*y[3] - x[3]*y[2],
                  x[3]*y[1] - x[1]*y[3],
                  x[1]*y[2] - x[2]*y[1])
end

# power
@inline Base.literal_pow(::typeof(^), x::AbstractSquareTensor, ::Val{-1}) = inv(x)
@inline Base.literal_pow(::typeof(^), x::AbstractSquareTensor, ::Val{0})  = one(x)
@inline Base.literal_pow(::typeof(^), x::AbstractSquareTensor, ::Val{1})  = x
@inline function Base.literal_pow(::typeof(^), x::AbstractSquareTensor, ::Val{p}) where {p}
    p > 0 ? (y = x; q = p) : (y = inv(x); q = -p)
    z = y
    for i in 2:q
        y = _powdot(y, z)
    end
    y
end
## helper functions
@inline _powdot(x::AbstractSecondOrderTensor, y::AbstractSecondOrderTensor) = dot(x, y)
@generated function _powdot(x::AbstractSymmetricSecondOrderTensor{dim}, y::AbstractSymmetricSecondOrderTensor{dim}) where {dim}
    _, exps = contraction_exprs(x, y, 1)
    quote
        @_inline_meta
        tensors = (x, y)
        @inbounds SymmetricSecondOrderTensor{dim}($(exps[indices(x)]...))
    end
end

# rotate
"""
    rotmat(θ::Real; degree::Bool = false)

Construct 2D rotation matrix.

# Examples
```jldoctest
julia> rotmat(30, degree = true)
2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 0.866025  -0.5
 0.5        0.866025
```
"""
@inline function rotmat(θ::Real; degree::Bool = false)
    if degree
        θ = deg2rad(θ)
    end
    sinθ = sin(θ)
    cosθ = cos(θ)
    @Mat [cosθ -sinθ
          sinθ  cosθ]
end

"""
    rotmat(θ::Vec{3}; sequence::Symbol, degree::Bool = false)
    rotmatx(θ::Real)
    rotmaty(θ::Real)
    rotmatz(θ::Real)

Convert Euler angles to rotation matrix.
Use 3 characters belonging to the set (X, Y, Z) for intrinsic rotations,
or (x, y, z) for extrinsic rotations.

# Examples
```jldoctest
julia> α, β, γ = rand(Vec{3});

julia> rotmat(Vec(α,β,γ), sequence = :XYZ) ≈ rotmatx(α) ⋅ rotmaty(β) ⋅ rotmatz(γ)
true

julia> rotmat(Vec(α,β,γ), sequence = :xyz) ≈ rotmatz(γ) ⋅ rotmaty(β) ⋅ rotmatx(α)
true

julia> rotmat(Vec(α,β,γ), sequence = :XYZ) ≈ rotmat(Vec(γ,β,α), sequence = :zyx)
true
```
"""
function rotmat(θ::Vec{3}; sequence::Symbol, degree::Bool = false)
    @inbounds α, β, γ = θ[1], θ[2], θ[3]
    if degree
        α, β, γ = map(deg2rad, (α, β, γ))
    end
    # intrinsic
    sequence == :XZX && return rotmatx(α) ⋅ rotmaty(β) ⋅ rotmatz(γ)
    sequence == :XYX && return rotmatx(α) ⋅ rotmaty(β) ⋅ rotmatx(γ)
    sequence == :YXY && return rotmaty(α) ⋅ rotmatx(β) ⋅ rotmaty(γ)
    sequence == :YZY && return rotmaty(α) ⋅ rotmatz(β) ⋅ rotmaty(γ)
    sequence == :ZYZ && return rotmatz(α) ⋅ rotmaty(β) ⋅ rotmatz(γ)
    sequence == :ZXZ && return rotmatz(α) ⋅ rotmatx(β) ⋅ rotmatz(γ)
    sequence == :XZY && return rotmatx(α) ⋅ rotmatz(β) ⋅ rotmaty(γ)
    sequence == :XYZ && return rotmatx(α) ⋅ rotmaty(β) ⋅ rotmatz(γ)
    sequence == :YXZ && return rotmaty(α) ⋅ rotmatx(β) ⋅ rotmatz(γ)
    sequence == :YZX && return rotmaty(α) ⋅ rotmatz(β) ⋅ rotmatx(γ)
    sequence == :ZYX && return rotmatz(α) ⋅ rotmaty(β) ⋅ rotmatx(γ)
    sequence == :ZXY && return rotmatz(α) ⋅ rotmatx(β) ⋅ rotmaty(γ)
    # extrinsic
    sequence == :xzx && return rotmatx(γ) ⋅ rotmaty(β) ⋅ rotmatz(α)
    sequence == :xyx && return rotmatx(γ) ⋅ rotmaty(β) ⋅ rotmatx(α)
    sequence == :yxy && return rotmaty(γ) ⋅ rotmatx(β) ⋅ rotmaty(α)
    sequence == :yzy && return rotmaty(γ) ⋅ rotmatz(β) ⋅ rotmaty(α)
    sequence == :zyz && return rotmatz(γ) ⋅ rotmaty(β) ⋅ rotmatz(α)
    sequence == :zxz && return rotmatz(γ) ⋅ rotmatx(β) ⋅ rotmatz(α)
    sequence == :xzy && return rotmaty(γ) ⋅ rotmatz(β) ⋅ rotmatx(α)
    sequence == :xyz && return rotmatz(γ) ⋅ rotmaty(β) ⋅ rotmatx(α)
    sequence == :yxz && return rotmatz(γ) ⋅ rotmatx(β) ⋅ rotmaty(α)
    sequence == :yzx && return rotmatx(γ) ⋅ rotmatz(β) ⋅ rotmaty(α)
    sequence == :zyx && return rotmatx(γ) ⋅ rotmaty(β) ⋅ rotmatz(α)
    sequence == :zxy && return rotmaty(γ) ⋅ rotmatx(β) ⋅ rotmatz(α)
    throw(ArgumentError("sequence $sequence is not supported"))
end
@inline function rotmatx(θ::Real; degree::Bool = false)
    if degree
        θ = deg2rad(θ)
    end
    o = one(θ)
    z = zero(θ)
    sinθ = sin(θ)
    cosθ = cos(θ)
    @Mat [o z     z
          z cosθ -sinθ
          z sinθ  cosθ]
end
@inline function rotmaty(θ::Real; degree::Bool = false)
    if degree
        θ = deg2rad(θ)
    end
    o = one(θ)
    z = zero(θ)
    sinθ = sin(θ)
    cosθ = cos(θ)
    @Mat [ cosθ z sinθ
           z    o z
          -sinθ z cosθ]
end
@inline function rotmatz(θ::Real; degree::Bool = false)
    if degree
        θ = deg2rad(θ)
    end
    o = one(θ)
    z = zero(θ)
    sinθ = sin(θ)
    cosθ = cos(θ)
    @Mat [cosθ -sinθ z
          sinθ  cosθ z
          z     z    o]
end

"""
    rotmat(a => b)

Construct rotation matrix rotating vector `a` to `b`.
The norms of two vectors must be the same.

# Examples
```jldoctest
julia> a = normalize(rand(Vec{3}))
3-element Vec{3, Float64}:
 0.526847334217759
 0.683741457787621
 0.5049054419691867

julia> b = normalize(rand(Vec{3}))
3-element Vec{3, Float64}:
 0.36698690362212083
 0.6333543148133657
 0.6813097125956302

julia> R = rotmat(a => b)
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 -0.594528   0.597477   0.538106
  0.597477  -0.119597   0.792917
  0.538106   0.792917  -0.285875

julia> R ⋅ a ≈ b
true
```
"""
function rotmat(pair::Pair{Vec{dim, T}, Vec{dim, T}})::Mat{dim, dim, T} where {dim, T}
    # https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/2672702#2672702
    a = pair.first
    b = pair.second
    dot(a, a) ≈ dot(b, b) || throw(ArgumentError("the norms of two vectors must be the same"))
    a ==  b && return  one(Mat{dim, dim, T})
    a == -b && return -one(Mat{dim, dim, T})
    c = a + b
    2 * (c ⊗ c) / (c ⋅ c) - one(Mat{dim, dim, T})
end

"""
    rotmat(θ, n::Vec; degree::Bool = false)

Construct rotation matrix from angle `θ` and direction `n`.

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

julia> n = Vec(0.0, 0.0, 1.0)
3-element Vec{3, Float64}:
 0.0
 0.0
 1.0

julia> rotmat(π/2, n) ⋅ x
3-element Vec{3, Float64}:
 1.1102230246251565e-16
 1.0
 0.0
```
"""
function rotmat(θ::Real, x::Vec{3}; degree::Bool = false)
    if degree
        θ = deg2rad(θ)
    end
    n = normalize(x)
    W = skew(n)
    I + W*sin(θ) + W^2*(1-cos(θ))
end

"""
    rotate(x::Vec, R::SecondOrderTensor)
    rotate(x::SecondOrderTensor, R::SecondOrderTensor)
    rotate(x::SymmetricSecondOrderTensor, R::SecondOrderTensor)

Rotate `x` by rotation matrix `R`.
This function can hold the symmetry of `SymmetricSecondOrderTensor`.

# Examples
```jldoctest
julia> A = rand(SymmetricSecondOrderTensor{3})
3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 0.590845  0.766797  0.566237
 0.766797  0.460085  0.794026
 0.566237  0.794026  0.854147

julia> R = rotmatz(π/4)
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

julia> rotate(A, R)
3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 -0.241332   0.0653796  -0.161071
  0.0653796  1.29226     0.961851
 -0.161071   0.961851    0.854147

julia> R ⋅ A ⋅ R'
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 -0.241332   0.0653796  -0.161071
  0.0653796  1.29226     0.961851
 -0.161071   0.961851    0.854147
```
"""
@inline rotate(v::Vec, R::SecondOrderTensor) = R ⋅ v
@inline rotate(v::Vec{2}, R::SecondOrderTensor{3}) = (R2x2 = Mat{2, 2}((i,j) -> @inbounds R[i,j]); R2x2 ⋅ v)
@inline rotate(A::SecondOrderTensor, R::SecondOrderTensor) = @einsum R[i,j] * A[j,k] * R[l,k]
@generated function rotate(A::SymmetricSecondOrderTensor{dim}, R::SecondOrderTensor{dim}) where {dim}
    _, exps = contraction_exprs(SecondOrderTensor{dim}, SecondOrderTensor{dim}, 1)
    TT = SymmetricSecondOrderTensor{dim}
    quote
        @_inline_meta
        ARᵀ = @einsum A[i,j] * R[k,j]
        tensors = (R, ARᵀ)
        @inbounds $TT($(exps[indices(TT)]...))
    end
end

function angleaxis(R::SecondOrderTensor{3})
    # https://math.stackexchange.com/questions/893984/conversion-of-rotation-matrix-to-quaternion
    θ = acos((tr(R)-1) / 2)
    n = Tensor(Real.(eigvecs(SArray(R))[:,3]))
    θ, n
end

# ----------------------------------------------#
# operations calling methods in StaticArrays.jl #
# ----------------------------------------------#

# eigvals/eigen
@inline function eigvals(x::AbstractSymmetricSecondOrderTensor; permute::Bool = true, scale::Bool = true)
    Tensor(eigvals(Symmetric(SArray(x)); permute, scale))
end
@inline function eigen(x::AbstractSymmetricSecondOrderTensor)
    eig = eigen(Symmetric(SArray(x)))
    Eigen(Tensor(eig.values), Tensor(eig.vectors))
end

# exp
@inline Base.exp(x::AbstractSecondOrderTensor) = typeof(x)(exp(SArray(x)))
@inline Base.exp(x::AbstractSymmetricSecondOrderTensor) = typeof(x)(exp(Symmetric(SArray(x))))

# diag/diagm
@inline diag(x::Union{AbstractMat, AbstractSymmetricSecondOrderTensor}, ::Val{k} = Val(0)) where {k} = Tensor(diag(SArray(x), Val{k}))
@inline diagm(kvs::Pair{<: Val, <: AbstractVec}...) = Tensor(diagm(map(kv -> kv.first => SArray(kv.second), kvs)...))
@inline diagm(x::Vec) = diagm(Val(0) => x)

# qr
struct QR{Q, R, P}
    Q::Q
    R::R
    p::P
end

# iteration for destructuring into components
Base.iterate(S::QR) = (S.Q, Val(:R))
Base.iterate(S::QR, ::Val{:R}) = (S.R, Val(:p))
Base.iterate(S::QR, ::Val{:p}) = (S.p, Val(:done))
Base.iterate(S::QR, ::Val{:done}) = nothing

for pv in (:true, :false)
    @eval function qr(A::AbstractMat, pivot::Val{$pv})
        F = qr(SArray(A), pivot)
        QR(Tensor(F.Q), Tensor(F.R), Tensor(F.p))
    end
end
qr(A::AbstractMat) = qr(A, Val(false))

# lu
struct LU{L, U, p}
    L::L
    U::U
    p::p
end

# iteration for destructuring into components
Base.iterate(S::LU) = (S.L, Val(:U))
Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p))
Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done))
Base.iterate(S::LU, ::Val{:done}) = nothing

for pv in (:true, :false)
    @eval function lu(A::AbstractMat, pivot::Val{$pv}; check = true)
        F = lu(SArray(A), pivot; check)
        LU(LowerTriangular(Tensor(parent(F.L))), UpperTriangular(Tensor(parent(F.U))), Tensor(F.p))
    end
end
lu(A::AbstractMat; check = true) = lu(A, Val(true); check)

# svd
struct SVD{T, TU, TS, TVt} <: Factorization{T}
    U::TU
    S::TS
    Vt::TVt
end
SVD(U::AbstractArray{T}, S::AbstractVector, Vt::AbstractArray{T}) where {T} = SVD{T, typeof(U), typeof(S), typeof(Vt)}(U, S, Vt)

@inline function Base.getproperty(F::SVD, s::Symbol)
    if s === :V
        return getfield(F, :Vt)'
    else
        return getfield(F, s)
    end
end
Base.propertynames(::SVD) = (:U, :S, :V, :Vt)

# iteration for destructuring into components
Base.iterate(S::SVD) = (S.U, Val(:S))
Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V))
Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done))
Base.iterate(S::SVD, ::Val{:done}) = nothing

function svd(A::AbstractMat; full=Val(false))
    F = svd(SArray(A); full)
    SVD(Tensor(F.U), Tensor(F.S), Tensor(F.Vt))
end


text:
module Characters
using LinearAlgebra
using SparseArrays
using Primes
using GroupsCore
using Cyclotomics

import PermutationGroups: AbstractOrbit, degree
import PermutationGroups

export AbstractClassFunction, Character, CharacterTable

export conjugacy_classes,
    constituents,
    degree,
    irreducible_characters,
    isirreducible,
    table

include("gf.jl")

include("echelon_form.jl")
include("eigenspacedecomposition.jl")
include("cmmatrix.jl")

include("powermap.jl")
include("character_tables.jl")
include("class_functions.jl")
include("dixon.jl")
end


text:
function UndoShift(AShot,datainfo)
	### Input:
	### AShot: one shot gather (nt,nrx,nry) with padding zeros
	### datainfo: a Dictionary contains info for one shot gather
	### Output:
	### A shot gather w/o linear shift
	nrx = datainfo["nrx"];
	nry = datainfo["nry"];
	Rgrid = datainfo["Rgrid"];
	dt = datainfo["dt"];

	vel = 1415; #velocity sand(dry) 0.2-1 Km/s, sand(water saturated 1.5-2 Km/s)

	# find the center position of receiver = source position
	CentX = (nrx+1)/2;
	CentY = (nry+1)/2;

	# calculate the shift grid for each receiver and save them in the matrix shiftT
	a = 1:nrx;
	b = 1:nry;
	row = reshape(repeat(a,outer=nry),nrx,nry);
	col = reshape(repeat(b,inner=nry),nrx,nry);
	shiftT = round.(sqrt.((row .-CentX).^2 .+ (col .-CentY).^2) .*(Rgrid/vel/(dt*1e-3)));

	nt1 =  datainfo["nt"] + 300 + 100;
	Shift_Back_Data = zeros(size(AShot));
	for i = 1:nrx
	    for j = 1:nry
		Shift_Back_Data[Int(shiftT[i,j]+1):nt1,i,j] = AShot[1:Int(nt1-shiftT[i,j]),i,j];
	    end
	end

	return Shift_Back_Data[300+1:300+datainfo["nt"],:,:]
end


text:
export unscale, unscale!

@doc """
    unscale!(S::GphysData[, chans=CC, irr=false])

Divide out the gains of all channels `i` where `S.fs[i] > 0.0`. Specify
`irr=true` to also remove the gains of irregularly-sampled channels. Use keyword
`chans=CC` to only resample channel numbers `CC`.

""" unscale!
function unscale!(S::GphysData;
  chans::ChanSpec=Int64[],
  irr::Bool=false)

  if chans == Int64[]
    chans = 1:S.n
  end

  proc_str = string("unscale!(S, chans=", chans, ")")

  @inbounds for i = 1:S.n
    (irr==false && S.fs[i]<=0.0) && continue
    if (S.gain[i] != 1.0) && (i in chans)
      T = eltype(S.x[i])
      rmul!(S.x[i], T(1.0/S.gain[i]))
      proc_note!(S, i, proc_str, string("divided out gain = ", repr(S.gain[i], context=:compact=>true)))
      S.gain[i] = 1.0
    end
  end
  return nothing
end
function unscale!(C::GphysChannel)
  rmul!(C.x, eltype(C.x)(1.0/C.gain))
  proc_note!(C, "unscale!(C)", string("divided out gain = ", repr(C.gain, context=:compact=>true)))
  C.gain = 1.0
  return nothing
end

@doc (@doc unscale!)
function unscale(S::GphysData;
  chans::ChanSpec=Int64[],
  irr::Bool=false)

  U = deepcopy(S)
  unscale!(U, chans=chans, irr=irr)
  return U
end
function unscale(C::GphysChannel; irr::Bool=false)
  U = deepcopy(C)
  rmul!(U.x, eltype(C.x)(1.0/U.gain))
  U.gain = 1.0
  return U
end


text:
module PerFlagFuncs

export split_by_consecutives
export count_flags
export find_next_iter
export merge_even_odd
export get_alternating_consecutive_vector

function split_by_consecutives(A::Vector{Int})
  N = length(A)
  if N>1
    sig = vcat(diff(A), [1])
    N_groups = count(map(x->x==1, sig))+1
    A_sub = Vector{Int}()
    A_split = Vector{Vector{Int}}([])
    for (a, s) in zip(A, sig)
      push!(A_sub, a)
      if !(s==1)
        push!(A_split, A_sub)
        A_sub = Vector{Int}([])
      end
      if a==A[N] && !(s==1)==false
        push!(A_split, A_sub)
      end
    end
    return A_split
  else
    return [A]
  end
end

"""
count_flags(s::String, flag::String)

counts the number of flags `flag` in string `s`.
"""
function count_flags(s::String, flag::String)
  cnt = 0
  LN = length(flag)
  LS = length(s)
  if LN <= LS
    for i in 1:LS-LN+1
      substr = s[i:i+LN-1]
      if flag==substr
        cnt += 1
      end
    end
  end
  return cnt
end

function find_next_iter(s::String, flag::String)
  ind = []
  LN = length(flag)
  LS = length(s)
  if LN <= LS
    for i in 1:LS-LN+1
      substr = s[i:i+LN-1]
      if flag==substr
        push!(ind, i)
      end
    end
  end
  return ind
end

function merge_even_odd(odd::Vector{Int64}, even::Vector{Int64})
  return [odd even]'[:]
end

function get_alternating_consecutive_vector(A::Vector{Int64},
                                            B::Vector{Int64},
                                            level_total_in=nothing,
                                            level_outer_in=nothing,
                                            s=nothing)
  e = Vector{Int64}(undef, 0)
  (C, D) = Tuple([e, e])
  if isempty(A)
    return (C, D)
  end
  if isempty(B)
    return (C, D)
  end
  N_AB = max(A..., B...)
  N_s = s == nothing ? N_AB : length(s)
  level_outer = level_outer_in==nothing ? zeros(N_s) : level_outer_in
  level_total = level_total_in==nothing ? zeros(N_s) : level_total_in
  L = Vector{Int64}(undef, 0)
  B_available = B
  if length(A) > 0 && length(B) > 0
    b_previous = B[1]
    for (i, a) in enumerate(A)
      found = false
      for (j, b) in enumerate(B_available)
        if !found
          start_match = level_outer[a]==1 && level_outer[a-1]==0
          stop_match  = level_outer[b]==1 && level_outer[b+1]==0
          cond_outer = level_outer_in == nothing ? true : start_match && stop_match
          cond_total = level_total[a]==level_total[b]
          cond_i1 = (b > a)
          cond_i2 = ( true || a == A[1] )
          cond_incr = cond_i1 && cond_i2
          cond = cond_incr && cond_total && cond_outer
          if cond
            push!(L, a)
            push!(L, b)
            filter!(x -> x > b_previous, B_available)
            b_previous = b
            found = true
          end
        end
      end
    end
    C = L[1:2:end]
    D = L[2:2:end]
  end
  return Tuple([C, D])
end

end


text:
module Pointy

using MinkowskiReduction, LinearAlgebra

#u = [2,1,0]; v = [3,4,6]; w = [0,-1,-2]

export pointgroup_basic, threeDrotation, pointgroup
function threeDrotation(u,v,w,α,β,γ)
    """
    threeDrotation(u,v,w,α,β,γ)

    Rotate a basis by three angles to any orientation

    """
A = [u v w]
R = [[cos(α)cos(β) cos(α)sin(β)sin(γ)-sin(α)cos(γ) cos(α)sin(β)cos(γ)+sin(α)sin(γ)];
     [sin(α)cos(β) sin(α)sin(β)sin(γ)+cos(α)cos(γ) sin(α)sin(β)cos(γ)-cos(α)sin(γ)];
     [-sin(β)      cos(β)sin(γ)                    cos(β)cos(γ)                   ]]
A = R*A
return A[:,1],A[:,2],A[:,3]
end

function pointgroup_basic(a1,a2,a3,debug=false)
    """
    Generate the symmetry operations of a lattice, defined by three 3-vectors
    (This function is _very_ simple but not efficient. Used to test more efficient algorithms.)

    ```juliadoctest
    julia> u = [1,0,0]; v = [.5,√3/2,0]; w = [0,0,√(8/3)];
    julia> pointgroup_basic(u,v,w)
    24-element Array{Array{Float64,2},1}:
    [-1.0 0.0 0.0; -1.0 1.0 0.0; 0.0 0.0 -0.9999999999999999]
    ...
    ```
    """
u,v,w = minkReduce(a1,a2,a3)
A = [u v w] # Put the lattice vectors as columns in matrix A
B = inv(A)*transpose(inv(A)) # Use this for checking for orthogonality
# A list of all possible lattice vectors in a rotated basis
c = [A*[i;j;k] for i ∈ (-1,0,1) for j ∈ (-1,0,1) for k ∈ (-1,0,1)]
# A list of all possible bases, (i.e., all combinations of c vectors)
R = [[i j k] for i ∈ c for j ∈ c for k ∈ c]
RT = [transpose(i) for i ∈ R]
# This is the U^T*U, where U transforms original basis to candidate basis
T = [R[i]*B*RT[i] for i ∈ 1:length(R)]
if debug return T end
# If T==identity then the U was a symmetry of the lattice
idx = findall([t≈I(3) for t ∈ T].==true)
Ai = inv(A)
ops = [Ai*R[i] for i in idx]
return ops
end

function pointgroup(a1,a2,a3)
    """
    Generate the symmetry operations of a lattice, defined by three 3-vectors.
    This function aims to be more efficient than `pointgroup_basic` and so is more complex

    ```juliadoctest
    julia> u = [1,0,0]; v = [.5,√3/2,0]; w = [0,0,√(8/3)];
    julia> pointgroup(u,v,w)
    24-element Array{Array{Float64,2},1}:
    [-1.0 0.0 0.0; -1.0 1.0 0.0; 0.0 0.0 -0.9999999999999999]
    ...
    ```
    """
u,v,w = minkReduce(a1,a2,a3) # Always do this first
A = [u v w] # Define a matrix with input vectors as columns
norms=norm.([u,v,w]) # Compute the norms of the three input vectors
vol = abs(u×v⋅w) # Volume of the parallelipiped formed by the basis vectors
AiAiT = inv(A)*transpose(inv(A)) # Use this for checking for orthogonality
# A list of all possible lattice vectors in a rotated basis
# These are lattice points from the vertices of the 8 cells with a corner at the origin)
# There are 27 of these (==3^3)
c = [A*[i;j;k] for i ∈ (-1,0,1) for j ∈ (-1,0,1) for k ∈ (-1,0,1)]
# Now keep only those vectors that have a norm matching one of the input vectors
c = c[findall([any(norm(i).≈norms) for i ∈ c])]
# Construct all possible bases, (i.e., all combinations of c vectors), skip duplicate vectors
R = [[i j k] for i ∈ c for j ∈ c if i ≠ j for k ∈ c if j ≠ k && i ≠ k]
R = R[findall([abs(det(r))≈vol for r in R])] # Delete candidate bases with the wrong volume
RT = [transpose(i) for i ∈ R]
# This is the U^T*U, where U transforms original basis to candidate basis
# If T==identity then the U was a symmetry of the lattice
T = [R[i]*AiAiT*RT[i] for i ∈ 1:length(R)]
# Indices of candidate T's that match the identity
idx = findall([t≈I(3) for t ∈ T].==true)
Ai = inv(A)
ops = [round.(Int,Ai*R[i]) for i in idx]
return ops
end
end


text:
# Keeps track of the variables that are in the current scope
struct Scope
    values::Dict{String, LLVM.AllocaInst}
    parent::Scope

    Scope() = new(Dict{String, LLVM.AllocaInst}())
    Scope(parent::Scope) = new(Dict{String, LLVM.AllocaInst}(), parent)
end

# Look for variable in `scope` and if not found, look recursively in parent
function _get(scope::Scope, var::String, default)
    v = get(scope.values, var, nothing)
    if v != nothing
        return v
    else
        return isdefined(scope, :parent) ? _get(scope.parent, var, default) : default
    end
end

mutable struct CurrentScope
    scope::Scope
end

CurrentScope() = CurrentScope(Scope())
open_scope!(currentscope::CurrentScope) = (currentscope.scope = Scope(currentscope.scope); currentscope)
Base.pop!(currentscope::CurrentScope) = (currentscope.scope = currentscope.scope.parent)
Base.setindex!(currentscope::CurrentScope, v::LLVM.AllocaInst, var::String) = currentscope.scope.values[var] = v
Base.get(currentscope::CurrentScope, var::String, default) = _get(currentscope.scope, var, default)


text:
using Base: Float64
"""
File with definitions of functions for structural analysis of aircraft configurations using beam elements
"""

Fmax = 1e15

"""
Get elasticity matrix for a single beam
"""
beam_get_K(
	L::Fg,
	EA::Fg,
	GJ::Fg,
	EIy::Fg,
	EIz::Fg
) where {Fg <: Real} = - Fg[
	(EA / L) 0.0 0.0 0.0 0.0 0.0 (- EA / L) 0.0 0.0 0.0 0.0 0.0;
	0.0 (12.0 * EIz / L ^ 3) 0.0 0.0 0.0 (6.0 * EIz / L ^ 2) 0.0 (- 12.0 * EIz / L ^ 3) 0.0 0.0 0.0 (6.0 * EIz / L ^ 2);
	0.0 0.0 (12.0 * EIy / L ^ 3) 0.0 (- 6.0 * EIy / L ^ 2) 0.0 0.0 0.0 (- 12.0 * EIy / L ^ 3) 0.0 (- 6.0 * EIy / L ^ 2) 0.0;
	0.0 0.0 0.0 (GJ / L) 0.0 0.0 0.0 0.0 0.0 (- GJ / L) 0.0 0.0;
	0.0 0.0 (- 6.0 * EIy / L ^ 2) 0.0 (4.0 * EIy / L) 0.0 0.0 0.0 (6.0 * EIy / L ^ 2) 0.0 (2.0 * EIy / L) 0.0;
	0.0 (6.0 * EIz / L ^ 2) 0.0 0.0 0.0 (4.0 * EIy / L) 0.0 (- 6.0 * EIz / L ^ 2) 0.0 0.0 0.0 (2.0 * EIz / L);
	(- EA / L) 0.0 0.0 0.0 0.0 0.0 (EA / L) 0.0 0.0 0.0 0.0 0.0;
	0.0 (- 12.0 * EIz / L ^ 3) 0.0 0.0 0.0 (- 6.0 * EIz / L ^ 2) 0.0 (12.0 * EIz / L ^ 3) 0.0 0.0 0.0 (- 6.0 * EIz / L ^ 2);
	0.0 0.0 (- 12.0 * EIy / L ^ 3) 0.0 (6.0 * EIy / L ^ 2) 0.0 0.0 0.0 (12.0 * EIy / L ^ 3) 0.0 (6.0 * EIy / L ^ 2) 0.0;
	0.0 0.0 0.0 ( - GJ / L) 0.0 0.0 0.0 0.0 0.0 (GJ / L) 0.0 0.0;
	0.0 0.0 (- 6.0 * EIy / L ^ 2) 0.0 (2.0 * EIy / L) 0.0 0.0 0.0 (6.0 * EIy / L ^ 2) 0.0 (4.0 * EIy / L) 0.0;
	0.0 (6.0 * EIz / L ^ 2) 0.0 0.0 0.0 (2.0 * EIz / L) 0.0 (- 6.0 * EIz / L ^ 2) 0.0 0.0 0.0 (4.0 * EIz / L)
]

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    GeoMeantoRelEntrBridge{T}

The `geometric mean cone` is representable with a relative entropy constraint
and a nonnegative auxiliary variable.

This is because ``u \\le \\prod_{i=1}^n w_i^{1/n}`` is equivalent to
``y \\ge 0`` and ``0 \\le u + y \\le \\prod_{i=1}^n w_i^{1/n}``,
and the latter inequality is equivalent to
``1 \\le \\prod_{i=1}^n (\\frac{w_i}{u + y})^{1/n}``, which is equivalent to
``0 \\le \\sum_{i=1}^n \\log (\\frac{w_i}{u + y})^{1/n}``, which is
equivalent to ``0 \\ge \\sum_{i=1}^n (u + y) \\log (\\frac{u + y}{w_i})``.

Thus ``(u, w) \\in GeometricMeanCone(1 + n)`` is representable as ``y \\ge 0``,
``(0, w, (u + y) e) \\in RelativeEntropyCone(1 + 2n)``, where ``e`` is a vector
of ones.
"""
struct GeoMeantoRelEntrBridge{T,F,G,H} <: AbstractBridge
    y::MOI.VariableIndex
    nn_index::MOI.ConstraintIndex{F,MOI.Nonnegatives} # for y >= 0
    relentr_index::MOI.ConstraintIndex{G,MOI.RelativeEntropyCone}
end

const GeoMeantoRelEntr{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GeoMeantoRelEntrBridge{T},OT}

function bridge_constraint(
    ::Type{GeoMeantoRelEntrBridge{T,F,G,H}},
    model::MOI.ModelLike,
    f::H,
    s::MOI.GeometricMeanCone,
) where {T,F,G,H}
    f_scalars = MOI.Utilities.eachscalar(f)
    (y, nn_index) = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    w_func = MOI.Utilities.vectorize(
        fill(
            MOI.Utilities.operate(+, T, f_scalars[1], y[1]),
            MOI.dimension(s) - 1,
        ),
    )
    relentr_func = MOI.Utilities.operate(
        vcat,
        T,
        zero(MOI.ScalarAffineFunction{Float64}),
        f_scalars[2:end],
        w_func,
    )
    relentr_index = MOI.add_constraint(
        model,
        relentr_func,
        MOI.RelativeEntropyCone(MOI.output_dimension(relentr_func)),
    )
    return GeoMeantoRelEntrBridge{T,F,G,H}(y[1], nn_index, relentr_index)
end

function MOI.supports_constraint(
    ::Type{<:GeoMeantoRelEntrBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    return true
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:GeoMeantoRelEntrBridge},
)
    return Tuple{Type}[(MOI.Nonnegatives,)]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:GeoMeantoRelEntrBridge{T,F,G}},
) where {T,F,G}
    return Tuple{Type,Type}[(G, MOI.RelativeEntropyCone)]
end

function concrete_bridge_type(
    ::Type{<:GeoMeantoRelEntrBridge{T}},
    H::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.GeometricMeanCone},
) where {T}
    F = MOI.VectorOfVariables
    S = MOI.Utilities.scalar_type(H)
    G = MOI.Utilities.promote_operation(
        vcat,
        T,
        T,
        S,
        MOI.Utilities.promote_operation(+, T, S, MOI.VariableIndex),
    )
    return GeoMeantoRelEntrBridge{T,F,G,H}
end

MOI.get(::GeoMeantoRelEntrBridge, ::MOI.NumberOfVariables)::Int64 = 1

function MOI.get(bridge::GeoMeantoRelEntrBridge, ::MOI.ListOfVariableIndices)
    return [bridge.y]
end

function MOI.get(
    ::GeoMeantoRelEntrBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.Nonnegatives},
)::Int64 where {T,F}
    return 1
end

function MOI.get(
    ::GeoMeantoRelEntrBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.RelativeEntropyCone},
)::Int64 where {T,F,G}
    return 1
end

function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.Nonnegatives},
) where {T,F}
    return [bridge.nn_index]
end

function MOI.get(
    bridge::GeoMeantoRelEntrBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.RelativeEntropyCone},
) where {T,F,G}
    return [bridge.relentr_index]
end

function MOI.delete(model::MOI.ModelLike, bridge::GeoMeantoRelEntrBridge)
    MOI.delete(model, bridge.relentr_index)
    MOI.delete(model, bridge.nn_index)
    MOI.delete(model, bridge.y)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::GeoMeantoRelEntrBridge{T,F,G,H},
) where {T,F,G,H}
    relentr_func = MOI.Utilities.eachscalar(
        MOI.get(model, MOI.ConstraintFunction(), bridge.relentr_index),
    )
    d = div(length(relentr_func) - 1, 2)
    u_func = MOI.Utilities.remove_variable(
        MOI.Utilities.operate(-, T, relentr_func[end], bridge.y),
        bridge.y,
    )
    w_func = relentr_func[2:(1+d)]
    return MOI.Utilities.convert_approx(
        H,
        MOI.Utilities.operate(vcat, T, u_func, w_func),
    )
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::GeoMeantoRelEntrBridge,
)
    return MOI.GeometricMeanCone(
        1 + div(
            MOI.dimension(
                MOI.get(model, MOI.ConstraintSet(), bridge.relentr_index),
            ) - 1,
            2,
        ),
    )
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{GeoMeantoRelEntrBridge{T,F,G,H}},
) where {T,F,G,H}
    ci_1 = MOI.ConstraintIndex{F,MOI.Nonnegatives}
    ci_2 = MOI.ConstraintIndex{G,MOI.RelativeEntropyCone}
    return MOI.supports(model, attr, ci_1) && MOI.supports(model, attr, ci_2)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::GeoMeantoRelEntrBridge,
)
    relentr_primal = MOI.get(model, attr, bridge.relentr_index)
    d = div(length(relentr_primal) - 1, 2)
    y_val = MOI.get(model, attr, bridge.nn_index)[1]
    u_val = sum(relentr_primal[(2+d):end]) / d - y_val
    return vcat(u_val, relentr_primal[2:(1+d)])
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::GeoMeantoRelEntrBridge,
    value,
)
    u_val = value[1]
    u_pos = -min(zero(u_val), u_val)
    d = length(value) - 1
    MOI.set(model, MOI.VariablePrimalStart(), bridge.y, u_pos)
    MOI.set(model, attr, bridge.nn_index, [u_pos])
    MOI.set(
        model,
        attr,
        bridge.relentr_index,
        vcat(0, value[2:end], fill(u_pos + u_val, d)),
    )
    return
end

# Given a is dual on y >= 0 and (b, c, d) is dual on RelativeEntropyCone
# constraint, dual on (u, w) in GeometricMeanCone is (-a, c). Note that
# sum(d) = -a, so we could instead use (sum(d), c).
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::GeoMeantoRelEntrBridge,
)
    u_dual = -MOI.get(model, attr, bridge.nn_index)[1]
    relentr_dual = MOI.get(model, attr, bridge.relentr_index)
    d = div(length(relentr_dual) - 1, 2)
    w_dual = relentr_dual[2:(d+1)]
    return vcat(u_dual, w_dual)
end

# Given GeometricMeanCone constraint dual start of (u, w), constraint dual on
# y >= 0 is -u and on RelativeEntropyCone constraint is
# (-u/n, w, u/n * (log.(w/geomean(w)) .+ 1)).
# Note log.(w/geomean(w)) = log.(w) .- sum(log, w) / n.
function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintDualStart,
    bridge::GeoMeantoRelEntrBridge,
    value,
)
    d = length(value) - 1
    u = value[1]
    w = value[2:end]
    MOI.set(model, MOI.ConstraintDualStart(), bridge.nn_index, [-u])
    relentr_dual =
        vcat(-u / d, value[2:end], u / d * (log.(w) .+ (1 - sum(log, w) / d)))
    MOI.set(
        model,
        MOI.ConstraintDualStart(),
        bridge.relentr_index,
        relentr_dual,
    )
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
"""
heq_aa()
Duplicate the part of Table 4 in Toth-Kelley for l^2 optimization
with c=.99.

Compare solution with one from nsoli.
"""
function heq_aa()
    c = 0.99
    tol = 1.e-8
    fcount = [11, 10, 10, 11, 12, 12]
    anormref = [4.0, 5.4, 5.4, 5.4, 5.4, 5.4]
    condref = [1.0, 2.e2, 1.9e5, 1.9e7, 5.5e9, 6.5e10]
    n = 500
    u0 = ones(n)
    hdata = heqinit(u0, c)
    FS = zeros(n)
    FPS = zeros(n, 20)
    rtol = tol
    atol = tol
    maxit = 100
    mmax = 6
    Vstore = zeros(n, 3 * mmax + 2)
    itrecords = zeros(6, 4)
    houtn = nsoli(
        heqf!,
        u0,
        FS,
        FPS;
        pdata = hdata,
        rtol = rtol,
        atol = atol,
        lmaxit = 10,
        eta = 0.01,
    )
    #
    # Solve H-equation with Anderson(m) for several values of m
    #
    for m = 1:6
        houta = aasol(
            HeqFix!,
            u0,
            m,
            Vstore;
            maxit = maxit,
            pdata = hdata,
            rtol = rtol,
            atol = atol,
        )
        #
        # Keep the books.
        #
        itrecords[m, 1] = hresults(houta.solution, houtn.solution)
        itrecords[m, 2] = length(houta.history)
        itrecords[m, 3] = norm(houta.stats.condhist, Inf)
        itrecords[m, 4] = norm(houta.stats.alphanorm, Inf)
    end
    #
    # Grade the results. I only use the coefficent norm and the condition
    # numbers for my own research.
    #
    solok = (norm(itrecords[:, 1], Inf) < 1.e-7)
    solok || println("Solution error in Anderson H solve")
    histok = (norm(itrecords[:, 2] - fcount, Inf) < 1.e-5)
    histok || println("History error in Anderson H solve")
    condok = (hresults(itrecords[:, 3], condref) < 1.e-1)
    condok || println("Condition error in Anderson H solve")
    normok = (hresults(itrecords[:, 4], anormref) < 1.e-1)
    normok || println("Coefficient norm error in Anderson H solve")
    #return (itrecords = itrecords)
    return solok && histok && condok && normok
end

function hresults(x, y)
    vdiff = (x - y) ./ abs.(x)
    return norm(vdiff, Inf)
end


text:
# ------------------------------------------------------------------------------
# documentation
# ------------------------------------------------------------------------------
using Documenter, DocumenterTools, Markdown, LFAToolkit
DocMeta.setdocmeta!(LFAToolkit, :DocTestSetup, :(using LFAToolkit); recursive = true)

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

# ------------------------------------------------------------------------------
# make
# ------------------------------------------------------------------------------
makedocs(
    modules = [LFAToolkit],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://jeremylt.github.io/LFAToolkit.jl/stable/",
        highlights = ["yaml"],
    ),
    pages = [
        "Introduction" => "index.md",
        "Mathematical Background" => "background.md",
        "Examples" => "examples.md",
        "Public API" => "public.md",
        "Private API" => "private.md",
        "Release Notes" => "release_notes.md",
        "References" => "references.md",
    ],
    clean = false,
    sitename = "LFAToolkit.jl",
    authors = "Jeremy L Thompson",
    linkcheck = !("skiplinks" in ARGS),
)

# ------------------------------------------------------------------------------
# deploy
# ------------------------------------------------------------------------------
deploydocs(
    repo = "github.com/jeremylt/LFAToolkit.jl.git",
    devbranch = "main",
    target = "build",
    push_preview = true,
)

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


text:
export AbstractMetricReader, CompositMetricReader, MetricReader, PeriodicMetricReader

abstract type AbstractMetricReader end

struct CompositMetricReader <: AbstractMetricReader
    readers::Vector{Any}
    function CompositMetricReader(readers...)
        new([readers...])
    end
end

function (r::CompositMetricReader)()
    for x in r.readers
        x()
    end
end

function shut_down!(r::CompositMetricReader)
    for x in r.readers
        shut_down!(x)
    end
end

#####

struct MetricReader{P,E} <: AbstractMetricReader
    provider::P
    exporter::E
    function MetricReader(provider::P, exporter::E) where {P,E}
        r = new{P,E}(provider, exporter)
        r()
        r
    end
end

MetricReader(; provider, exporter) = MetricReader(provider, exporter)

function (r::MetricReader)()
    for ins in keys(r.provider.async_instruments)
        ins()
    end
    export!(r.exporter, (m for m in values(r.provider.metrics)))
end

# ??? shut_down! provider?
shut_down!(r::MetricReader) = shut_down!(r.exporter)

#####

struct PeriodicMetricReader{R<:AbstractMetricReader} <: AbstractMetricReader
    reader::R
    export_interval_seconds::Int
    export_timeout_seconds::Int
    timer::Timer
    function PeriodicMetricReader(
        reader;
        export_interval_seconds = 60,
        export_timeout_seconds = 30,
    )
        timer = Timer(0; interval = export_interval_seconds) do t
            res = timedwait(export_timeout_seconds; pollint = 1) do
                reader()
                true
            end
            if res == :timed_out
                @warn "timed out when exporting metrics"
                # how to stop exporting?
            end
        end

        new{typeof(reader)}(reader, export_interval_seconds, export_timeout_seconds, timer)
    end
end

function shut_down!(r::PeriodicMetricReader)
    close(r.timer)
    shut_down!(r.reader)
end


text:
# GenBank Writer

struct Writer{S <: TranscodingStream} <: BioGenerics.IO.AbstractWriter
    output::S
end

function BioGenerics.IO.stream(writer::Writer)
    return writer.output
end

"""
    GenBank.Writer(output::IO; width=70)

Create a data writer of the GenBank file format.
```julia
open(GenBank.Writer, outfile) do writer
    write(writer, genome)
end
```
"""
function Writer(output::IO)
    if output isa TranscodingStream
        return Writer{typeof(output)}(stream)
    else
        stream = TranscodingStreams.NoopStream(output)
        return Writer{typeof(stream)}(stream)
    end
end

function Base.write(writer::Writer, record::Record)
    printgbk(writer.output, record)
end

"""
    printgbk(io::IO, chr)
    printgbk(path::AbstractString, chr)

Print `chr` in GenBank format.
"""
function printgbk(path::AbstractString, chrs)
    open(path, "w") do f
        printgbk(f, chrs)
    end
end
function printgbk(chrs::AbstractVector{C}) where {C <: Record}
    io = IOBuffer()
    printgbk(io, chrs)
end
function printgbk(io::IO, chrs::AbstractVector{C}) where {C <: Record}
    for chr in chrs
        printgbk(io, chr)
    end
    return io
end
function printgbk(chr::C) where {C <: Record}
    io = IOBuffer()
    printgbk(io, chr)
end
function printgbk(io::IO, chr::C) where {C <: Record}
    if chr.header[1] != '#'
        println(io, chr.header)
    else
        ### GFF3 header
        println(io, "LOCUS       unknown", lpad(length(chr.sequence), 10, ' '), " bp     DNA")
    end
    println(io, rpad("FEATURES", 21, ' '), "Location/Qualifiers")
    println(io, chr.genes)
    println(io, "ORIGIN")
    formatsequence(chr.sequence, io)
    println(io)
    println(io, "//")
    return io
end

function formatsequence(sequence, io = IOBuffer)
    p = length(string(length(sequence))) + 2
    if length(sequence) > 60
        intervals = [i:i+60 for i in range(1; step = 60, stop = length(sequence)-60)]
        for interval in intervals
            println(io, lpad(string(first(interval)), p, ' '), " ", sequence[interval[1:10]],
                " ", sequence[interval[11:20]], " ", sequence[interval[21:30]],
                " ", sequence[interval[31:40]], " ", sequence[interval[41:50]],
                " ", sequence[interval[51:60]])
        end
    else
        intervals = [1:1]
    end
    i = intervals[end].stop
    if i <= length(sequence)
        print(io, lpad(i, p, ' '), " ")
        j = 0
        while i+j <= length(sequence)
            print(io, sequence[i+j])
            (j+1) % 10 == 0 && print(io, " ")
            j += 1
        end
    end
    return io
end


text:
# Structure for VMLM of OptimPack, see OptimPack.jl
# https://github.com/emmt/OptimPack.jl
ls = OptimPack.MoreThuenteLineSearch(ftol = 1e-8, gtol = 0.95)
gat = 0
grt = 1e-6
vt = false
memsize = 10
mxvl = 1000
mxtr = 1000


text:
abstract type TimeError <: Exception end

"""
    AmbiguousTimeError(local_datetime, timezone)

The provided local datetime is ambiguous within the specified timezone. Typically occurs on
daylight saving time transitions which "fall back" causing duplicate hour long period.
"""
struct AmbiguousTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

function Base.showerror(io::IO, e::AmbiguousTimeError)
    print(io, "AmbiguousTimeError: ")
    print(io, "Local DateTime $(e.local_dt) is ambiguous within $(string(e.timezone))")
end

"""
    NonExistentTimeError(local_datetime, timezone)

The provided local datetime is does not exist within the specified timezone. Typically
occurs on daylight saving time transitions which "spring forward" causing an hour long
period to be skipped.
"""
struct NonExistentTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

function Base.showerror(io::IO, e::NonExistentTimeError)
    print(io, "NonExistentTimeError: ")
    print(io, "Local DateTime $(e.local_dt) does not exist within $(string(e.timezone))")
end

"""
    UnhandledTimeError(timezone)

The timezone calculation occurs beyond the last calculated transition.
"""
struct UnhandledTimeError <: TimeError
    tz::VariableTimeZone
end

function Base.showerror(io::IO, e::UnhandledTimeError)
    print(io, "UnhandledTimeError: ")
    print(io, "TimeZone $(string(e.tz)) does not handle dates on or after $(e.tz.cutoff) UTC")
end

# Note: ParseNextError should avoid raising exceptions when start/end index aren't
# unreasonable to attempt and always show the actual parsing error instead of an internal
# error.

"""
    ParseNextError

An exception which displays the portion of a string which could not be parsed successfully.
"""
struct ParseNextError <: Exception
    msg::AbstractString
    str::AbstractString
    s::Int
    e::Int
end

function ParseNextError(msg, str, s)
    len = lastindex(str)
    ParseNextError(msg, str, s, s >= len ? len : s)
end

function Base.showerror(io::IO, e::ParseNextError)
    print(io, "$ParseNextError: ")
    if !isempty(e.msg)
        print(io, "$(e.msg): ")
    end

    str = e.str
    len = lastindex(str)
    u_start = e.s
    u_end = e.e > len || e.e < e.s ? len : e.e

    print(io, "\"", str[firstindex(str):prevind(str, u_start)])
    printstyled(io, str[u_start:u_end], color=:underline)
    print(io, str[nextind(str, u_end):len])
    if e.s > len
        printstyled(io, "\"", color=:underline)
    else
        print(io, "\"")
    end
end

function Base.show(io::IO, ::MIME"text/plain", e::ParseNextError)
    showerror(io, e)
    print(io, " ($(e.s), $(e.e))")
end


text:
# back propagation
export update!

# update
function update!(v::voidLayer,t::Int64)
end
function update!(d::dataLayer,t::Int64)
end
function update!(h::hiddenLayer,t::Int64)
  ∂W=reshape(mean(h.∂w_L.+λ*regularize(h.W),ndims(h.∂w_L)),size(h.∂W)...)
  ∂b=reshape(mean(h.∂b_L,ndims(h.∂b_L)),size(h.∂b)...)
  h.∂W=momentum*h.∂W+(1-momentum)*clamp!(∂W,-▽max,▽max)
  h.∂b=momentum*h.∂b+(1-momentum)*clamp!(∂b,-▽max,▽max)
  h.W+=-lr/(1+δ*t)*h.∂W;h.b+=-lr/(1+δ*t)*h.∂b;
end
function update!(p::predictLayer,t::Int64)
end
function update!(l::lossLayer,t::Int64)
end
function update!(net::Array{Layer},t::Int64)
  for idx=1:1:length(net)
    update!(net[idx],t::Int64)
  end
end


text:
include("../imports.jl")

@testset ExtendedTestSet "Bounds" begin
    @testset ExtendedTestSet "boundsranges" begin
        keypoints = Keypoints(
            [SVector(1, 1)],
            (50, 50)
        )
        image = Image(rand(RGB, 50, 50))
        @test boundsranges(getbounds(keypoints)) == boundsranges(getbounds(image))
    end
end


text:
using LazIO
using ProgressMeter

function correct!(points::Vector{LazIO.LazPoint})
    sort!(points, by=x->(x.gps_time, LazIO.return_number(x)))
    prevpoint = points[1]
    n = 0
    first = points[1]
    @inbounds @showprogress "Correcting " for i in eachindex(points)
        point = points[i]
        Δt = point.gps_time - prevpoint.gps_time
        nr = LazIO.number_of_returns(point)
        rn = LazIO.return_number(point)
        # println("$rn of $nr with diff $Δt at $i")
        if (nr > 1) && (rn == n+1)
            if rn == 1
                # @info "first $i"
                first = point
                point.classification = 19
            elseif Δt < 1e-7
                # @info "corr $i"
                try
                    correct!(first, point)
                    point.classification = 20
                catch e
                    @warn "Can't correct $i"
                    point.classification = 21
                end
            else
                point.classification = 22
            end

            # point.intensity /= first.intensity

            if nr == rn
                n = 0
            else
                n += 1
            end
        else
            # @info "skip $i"
            n = 0
            point.classification = 10
        end
        points[i] = point
        prevpoint = point
    end
    points
end

function write!(fn, ds, points)
    LazIO.write(fn, ds.header) do io
        for p in points
            LazIO.writepoint(io, p)
        end
    end
end


text:
export image, dockerfile

abstract AbstractImage
type Image <: AbstractImage
  name::String
end

type BuildImage <: AbstractImage
  name::String
  base::AbstractImage

  volumes::Vector{String}
  ppas::Vector
  workdir::String
  env::Dict
  packages::Vector{String}
  pips::Vector{String}
  args::Vector{String}
  rm::Bool
end

const DEFAULT_IMAGE = convert(Dict{Symbol, Any}, {
  :volumes => String[],
  :workdir => "",
  :env => Dict(),
  :packages => String[],
  :pips => String[],
  :args => String[],
  :rm => true
})

function BuildImage(name::String, base::AbstractImage; kwargs...)
  # Treat volume separately
  volumes = Util.getarg(:volumes, kwargs, DEFAULT_IMAGE)
  for (vol, path) in filter(x -> x[1] == :volume, kwargs)
    push!(volumes, string(path))
  end
  filter!(x -> x[1] != :volume, kwargs)

  ppas = [ppa for (key, ppa) in filter(x -> x[1] == :ppa, kwargs)]
  filter!(x -> x[1] != :ppa, kwargs)

  # Remove special keyword arg volumes
  if length(kwargs) !=(map(x->x[1], kwargs) |> unique |> length)
    error("Keyword arguments, other than volume should be unique")
  end

  BuildImage(
    name, base, unique(volumes), ppas,
    map(x -> Util.getarg(x, kwargs, DEFAULT_IMAGE), names(BuildImage)[5:end])...
  )
end

function BuildImage(name::String, base::String="ubuntu:14.04"; kwargs...)
  if length(base) == ""
    error("Base cannot be empty")
  end
  BuildImage(name, Image(base); kwargs...)
end

function image(name::String, args...; kwargs...)
  if length(args) != 0 || length(kwargs) != 0
    BuildImage(name, args...; kwargs...)
  else
    Image(name)
  end
end

function dockerfile(image::BuildImage)
    result = "FROM $(image.base.name)\n"
    result *= "RUN apt-get update\n"
    if length(image.ppas) > 0
      ppas = join(["add-apt-repository ppa:" * n for n in image.ppas], "\\\n    && ")
      packages = "software-properties-common python-software-properties"
      result *= "RUN apt-get install -y $packages \\\n    && $ppas \\\n    && apt-get update\n"
    end
    if length(image.packages) > 0
      packages = ""
      regex = r"([a-z,A-Z,0-9,_,\-,+,\.]*)\s*(?:\(\s*((?:>|<|=)*)\s*([0-9,.,-,+]*)\s*\))?"
      for package in image.packages
        name, constraint, version = match(regex, package).captures
        if is(name, nothing) || length(name) == 0
          error("Could not determine name of package")
        end
        if (!is(version, nothing)) && length(version) > 0
          if is(constraint, nothing)
            packages *= " " * name * "=" * version
          elseif constraint == "==" || constraint == "="
            packages *= " " * name * "=" * version
          else
            packages *= " " * name
          end
        else
          packages *= " " * name
        end
      end
      result = result *"RUN apt-get install -y $packages\n"
    end
    if length(image.pips) > 0
        result = result * "RUN pip install $(join(image.pips, " "))\n"
    end
    if length(image.volumes) > 0
        volumes = join(image.volumes, " ")
        result *= "RUN mkdir -p $volumes \n"
        result *= "VOLUME $volumes\n"
    end
    if length(image.env) > 0
      result *= "ENV " * join(["$(string(k))=\"$(string(v))\"" for (k, v) in image.env], " ")
      result *= "\n"
    end
    if length(image.args) > 0
        result = result * join(["RUN $(chomp(u))" for u in image.args], "\n")
    end
    if length(image.workdir) > 0
        result *= "WORKDIR $(image.workdir)\n"
    end
    result
end

|> (machine::MachineEnv, image::AbstractImage) = (machine, image)
function |> (parent::AbstractImage, child::BuildImage)
  result = deepcopy(child)
  result.base = parent
  result
end
function |> (parents::(MachineEnv, AbstractImage), child::BuildImage)
  parents[1] |> (parents[2] |> child)
end

function command_impl(machine::MachineEnv, image::BuildImage)
  rmargs = if image.rm "--rm=true" else "" end
  name = if length(image.name) > 0 "--tag=\"$(image.name)\"" else "" end
  # File needs to be in context... there is no good way to get a temp file in a specific dir yet, it
  # seems.
  imagefile = Util.heretempfile()
  open(imagefile, "w") do file
    write(file, dockerfile(image))
  end
  filename = "--file=\"$(basename(imagefile))\""
  path = dirname(imagefile)
  cmd = command(machine, "build", "$rmargs  $name $filename $path")
  if isa(image.base, BuildImage)
    cmds, files = command_impl(machine, image.base)
    [cmds..., cmd], [files..., imagefile]
  else
    [cmd], [imagefile]
  end
end

function command_impl(func::Function, machine::MachineEnv, image::BuildImage)
  cmds, files = command_impl(machine, image)
  try
    return cmds |> func
  finally
    for imagefile in files
      (!isfile(imagefile)) || rm(imagefile)
    end
  end
end

run(inputs::(MachineEnv, BuildImage)) = run(inputs[1], inputs[2])
run(machine::MachineEnv, image::BuildImage) = command_impl(run, machine, image)
readchomp(inputs::(MachineEnv, BuildImage)) = readchomp(inputs[1], inputs[2])
readchomp(machine::MachineEnv, image::BuildImage) = command_impl(readchomp, machine, image)
readall(inputs::(MachineEnv, BuildImage)) = readall(inputs[1], inputs[2])
readall(machine::MachineEnv, image::BuildImage) = command_impl(readall, machine, image)


text:
# Mag
set!(res::MagLike, ::Irrational{:π}) = const_pi!(res)

# Arb
function set!(res::ArbLike, x::Union{MagLike,BigFloat})
    set!(radref(res), UInt64(0))
    set!(midref(res), x)
    return res
end

function set!(res::ArbLike, x::Rational; prec::Integer = precision(res))
    set!(res, numerator(x))
    return div!(res, res, denominator(x), prec = prec)
end

for (irr, suffix) in ((:π, "pi"), (:ℯ, "e"), (:γ, "euler"), (:catalan, "catalan"))
    jlf = Symbol("const_$suffix", "!")
    IrrT = Irrational{irr}
    @eval begin
        set!(res::ArbLike, ::$IrrT; prec::Integer = precision(res)) = $jlf(res, prec = prec)
    end
end

function set!(res::ArbLike, ::Irrational{:φ}; prec::Integer = precision(res))
    set!(res, 5)
    sqrt!(res, res, prec)
    add!(res, res, 1, prec)
    return div!(res, res, 2, prec)
end

function set!(
    res::ArbLike,
    (a, b)::NTuple{2,Union{MagLike,ArfLike,BigFloat}};
    prec::Integer = precision(res),
)
    a <= b || throw(ArgumentError("must have a <= b, got a = $a and b = $b"))
    return set_interval!(res, a, b, prec = prec)
end

function set!(res::ArbLike, (a, b)::Tuple{<:Real,<:Real}; prec::Integer = precision(res))
    # TODO: This is not strictly required to check.
    a > b && throw(ArgumentError("must have a <= b, got a = $a and b = $b"))
    # TODO: If we really want to we could avoid one allocation by reusing res
    return union!(res, Arb(a, prec = prec), Arb(b, prec = prec), prec = prec)
end

# Acb
function set!(res::AcbLike, x::Union{Real,arf_struct,mag_struct,Tuple{<:Real,<:Real}})
    set!(realref(res), x)
    set!(imagref(res), 0)
    return res
end

function set!(
    res::AcbLike,
    re::Union{Real,arb_struct,arf_struct,mag_struct,Tuple{<:Real,<:Real}},
    im::Union{Real,arb_struct,arf_struct,mag_struct,Tuple{<:Real,<:Real}},
)
    set!(realref(res), re)
    set!(imagref(res), im)
    return res
end

set!(res::AcbLike, z::Complex) = set!(res, real(z), imag(z))
set!(res::AcbLike, ::Irrational{:π}) = const_pi!(res)


text:
using JsonGrinder, JSON, Test, SparseArrays, Mill

@testset "ExtractMultipleRepresentation" begin
	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"]),
		JsonGrinder.ExtractString(String)))
	e = ex("Olda")

	@test length(e.data) == 2
	@test size(e.data[1].data) == (4, 1)
	@test e.data[1].data[:] ≈ [0, 1, 0, 0]
	@test size(e.data[2].data) == (2053, 1)
	@test findall(x->x > 0, SparseMatrixCSC(e.data[2].data)) .|> Tuple == [(80, 1), (98, 1), (206, 1), (738, 1),
		(1062, 1), (1856, 1)]
end

@testset "show" begin
	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"]),
		JsonGrinder.ExtractString(String)))
	e = ex("Olda")

	buf = IOBuffer()
	Base.show(buf, ex)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
: MultiRepresentation
 : Categorical d = 4
 : String
"""

	buf = IOBuffer()
	Mill.dsprint(buf, e)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
TreeNode
  ├── ArrayNode(4, 1)
  └── ArrayNode(2053, 1)"""
end


text:
# HypothesisTests.jl
# Hypothesis tests in Julia
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

VERSION >= v"0.4.0-dev+6521" && __precompile__()

module HypothesisTests

using Distributions, Roots, StatsBase, Compat
using Combinatorics: combinations
using Rmath: pwilcox, psignrank

export testname, pvalue, ci
abstract HypothesisTest

check_same_length(x::AbstractVector, y::AbstractVector) = if length(x) != length(y)
    throw(DimensionMismatch("Vectors must be the same length"))
end

# Basic function for finding a p-value given a distribution and tail
pvalue(dist::ContinuousUnivariateDistribution, x::Number; tail=:both) =
    if tail == :both
        min(2 * min(cdf(dist, x), ccdf(dist, x)), 1.0)
    elseif tail == :left
        cdf(dist, x)
    elseif tail == :right
        ccdf(dist, x)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end

pvalue(dist::DiscreteUnivariateDistribution, x::Number; tail=:both) =
    if tail == :both
        min(2 * min(ccdf(dist, x-1), cdf(dist, x)), 1.0)
    elseif tail == :left
        cdf(dist, x)
    elseif tail == :right
        ccdf(dist, x-1)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end

function check_alpha(alpha::Float64)
    if alpha <= 0 || alpha >= 0.5
        throw(ArgumentError("alpha $alpha not in range (0, 0.5)"))
    end
end

function significance{T<:Real}(p_val::T)
    if p_val > 0.05
        "not significant"
    elseif p_val < 0.001
        "extremely significant"
    elseif  p_val < 0.01
        "very significant"
    else
        "significant"
    end
end

# Pretty-print
function Base.show{T<:HypothesisTest}(io::IO, test::T)
    println(io, testname(test))
    println(io, repeat("-", length(testname(test))))

    # population details
    has_ci = applicable(ci, test)
    (param_name, param_under_h0, param_estimate) = population_param_of_interest(test)
    println(io, "Population details:")
    println(io, "    parameter of interest:   $param_name")
    println(io, "    value under h_0:         $param_under_h0")
    println(io, "    point estimate:          $param_estimate")
    if has_ci
        println(io, "    95% confidence interval: $(ci(test))")
    end
    println(io)

    # test summary
    p = pvalue(test)
    outcome = if p > 0.05 "fail to reject" else "reject" end
    println(io, "Test summary:")
    println(io, "    outcome with 95% confidence: $outcome h_0")
    println(io, "    two-sided p-value:           $p ($(significance(p)))")
    println(io)

    # further details
    println(io, "Details:")
    show_params(io, test, "    ")
end

# parameter of interest: name, value under h0, point estimate
population_param_of_interest{T<:HypothesisTest}(test::T) = ("not implemented yet", NaN, NaN)

function show_params{T<:HypothesisTest}(io::IO, test::T, ident="")
    fieldidx = find(Bool[t<:Number for t in T.types])
    if !isempty(fieldidx)
        lengths = [length(string(T.names[i])) for i in fieldidx]
        maxlen = maximum(lengths)

        for i = 1:length(fieldidx)
            name = T.names[fieldidx[i]]
            println(io, ident, repeat(" ", maxlen-lengths[i]),
                      replace(string(name), "_", " "),
                      " = $(getfield(test, name))")
        end
    end
end

include("common.jl")

include("binomial.jl")
include("circular.jl")
include("fisher.jl")
include("kolmogorov_smirnov.jl")
include("kruskal_wallis.jl")
include("mann_whitney.jl")
include("t.jl")
include("z.jl")
include("wilcoxon.jl")
include("power_divergence.jl")
include("anderson_darling.jl")
end