text:
function FieldUnivariateDerivative2(arg0::FieldDerivativeStructure)
    return FieldUnivariateDerivative2((FieldDerivativeStructure,), arg0)
end

function FieldUnivariateDerivative2(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return FieldUnivariateDerivative2((RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function abs(obj::FieldUnivariateDerivative2)
    return jcall(obj, "abs", FieldUnivariateDerivative2, ())
end

function acos(obj::FieldUnivariateDerivative2)
    return jcall(obj, "acos", FieldUnivariateDerivative2, ())
end

function acosh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "acosh", FieldUnivariateDerivative2, ())
end

function add(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "add", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function add(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "add", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function add(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "add", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function asin(obj::FieldUnivariateDerivative2)
    return jcall(obj, "asin", FieldUnivariateDerivative2, ())
end

function asinh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "asinh", FieldUnivariateDerivative2, ())
end

function atan(obj::FieldUnivariateDerivative2)
    return jcall(obj, "atan", FieldUnivariateDerivative2, ())
end

function atan2(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "atan2", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function atanh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "atanh", FieldUnivariateDerivative2, ())
end

function cbrt(obj::FieldUnivariateDerivative2)
    return jcall(obj, "cbrt", FieldUnivariateDerivative2, ())
end

function ceil(obj::FieldUnivariateDerivative2)
    return jcall(obj, "ceil", FieldUnivariateDerivative2, ())
end

function compose(obj::FieldUnivariateDerivative2, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(obj, "compose", FieldUnivariateDerivative2, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function copy_sign(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "copySign", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function copy_sign(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "copySign", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function copy_sign(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "copySign", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function cos(obj::FieldUnivariateDerivative2)
    return jcall(obj, "cos", FieldUnivariateDerivative2, ())
end

function cosh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "cosh", FieldUnivariateDerivative2, ())
end

function divide(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "divide", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function divide(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "divide", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function divide(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "divide", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function equals(obj::FieldUnivariateDerivative2, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::FieldUnivariateDerivative2)
    return jcall(obj, "exp", FieldUnivariateDerivative2, ())
end

function expm1(obj::FieldUnivariateDerivative2)
    return jcall(obj, "expm1", FieldUnivariateDerivative2, ())
end

function floor(obj::FieldUnivariateDerivative2)
    return jcall(obj, "floor", FieldUnivariateDerivative2, ())
end

function get_derivative(obj::FieldUnivariateDerivative2, arg0::jint)
    return jcall(obj, "getDerivative", RealFieldElement, (jint,), arg0)
end

function get_exponent(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getField", FieldUnivariateDerivative2Field, ())
end

function get_first_derivative(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getFirstDerivative", RealFieldElement, ())
end

function get_order(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getOrder", jint, ())
end

function get_real(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getReal", jdouble, ())
end

function get_second_derivative(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getSecondDerivative", RealFieldElement, ())
end

function get_value(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getValue", RealFieldElement, ())
end

function get_value_field(obj::FieldUnivariateDerivative2)
    return jcall(obj, "getValueField", Field, ())
end

function hash_code(obj::FieldUnivariateDerivative2)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "hypot", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2, arg1::FieldUnivariateDerivative2, arg2::FieldUnivariateDerivative2, arg3::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2, arg1::FieldUnivariateDerivative2, arg2::FieldUnivariateDerivative2, arg3::FieldUnivariateDerivative2, arg4::FieldUnivariateDerivative2, arg5::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2, arg1::FieldUnivariateDerivative2, arg2::FieldUnivariateDerivative2, arg3::FieldUnivariateDerivative2, arg4::FieldUnivariateDerivative2, arg5::FieldUnivariateDerivative2, arg6::FieldUnivariateDerivative2, arg7::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::RealFieldElement, arg1::FieldUnivariateDerivative2, arg2::RealFieldElement, arg3::FieldUnivariateDerivative2, arg4::RealFieldElement, arg5::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (RealFieldElement, FieldUnivariateDerivative2, RealFieldElement, FieldUnivariateDerivative2, RealFieldElement, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::Vector{FieldUnivariateDerivative2}, arg1::Vector{FieldUnivariateDerivative2})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (Vector{FieldUnivariateDerivative2}, Vector{FieldUnivariateDerivative2}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::Vector{RealFieldElement}, arg1::Vector{FieldUnivariateDerivative2})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (Vector{RealFieldElement}, Vector{FieldUnivariateDerivative2}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::Vector{jdouble}, arg1::Vector{FieldUnivariateDerivative2})
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (Vector{jdouble}, Vector{FieldUnivariateDerivative2}), arg0, arg1)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::jdouble, arg1::FieldUnivariateDerivative2, arg2::jdouble, arg3::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::jdouble, arg1::FieldUnivariateDerivative2, arg2::jdouble, arg3::FieldUnivariateDerivative2, arg4::jdouble, arg5::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::FieldUnivariateDerivative2, arg0::jdouble, arg1::FieldUnivariateDerivative2, arg2::jdouble, arg3::FieldUnivariateDerivative2, arg4::jdouble, arg5::FieldUnivariateDerivative2, arg6::jdouble, arg7::FieldUnivariateDerivative2)
    return jcall(obj, "linearCombination", FieldUnivariateDerivative2, (jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2, jdouble, FieldUnivariateDerivative2), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::FieldUnivariateDerivative2)
    return jcall(obj, "log", FieldUnivariateDerivative2, ())
end

function log10(obj::FieldUnivariateDerivative2)
    return jcall(obj, "log10", FieldUnivariateDerivative2, ())
end

function log1p(obj::FieldUnivariateDerivative2)
    return jcall(obj, "log1p", FieldUnivariateDerivative2, ())
end

function multiply(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "multiply", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function multiply(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "multiply", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function multiply(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "multiply", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function multiply(obj::FieldUnivariateDerivative2, arg0::jint)
    return jcall(obj, "multiply", FieldUnivariateDerivative2, (jint,), arg0)
end

function negate(obj::FieldUnivariateDerivative2)
    return jcall(obj, "negate", FieldUnivariateDerivative2, ())
end

function new_instance(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "newInstance", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function pow(::Type{FieldUnivariateDerivative2}, arg0::jdouble, arg1::FieldUnivariateDerivative2)
    return jcall(FieldUnivariateDerivative2, "pow", FieldUnivariateDerivative2, (jdouble, FieldUnivariateDerivative2), arg0, arg1)
end

function pow(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "pow", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function pow(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "pow", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function pow(obj::FieldUnivariateDerivative2, arg0::jint)
    return jcall(obj, "pow", FieldUnivariateDerivative2, (jint,), arg0)
end

function reciprocal(obj::FieldUnivariateDerivative2)
    return jcall(obj, "reciprocal", FieldUnivariateDerivative2, ())
end

function remainder(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "remainder", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function remainder(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "remainder", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function remainder(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "remainder", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function rint(obj::FieldUnivariateDerivative2)
    return jcall(obj, "rint", FieldUnivariateDerivative2, ())
end

function root_n(obj::FieldUnivariateDerivative2, arg0::jint)
    return jcall(obj, "rootN", FieldUnivariateDerivative2, (jint,), arg0)
end

function scalb(obj::FieldUnivariateDerivative2, arg0::jint)
    return jcall(obj, "scalb", FieldUnivariateDerivative2, (jint,), arg0)
end

function signum(obj::FieldUnivariateDerivative2)
    return jcall(obj, "signum", FieldUnivariateDerivative2, ())
end

function sin(obj::FieldUnivariateDerivative2)
    return jcall(obj, "sin", FieldUnivariateDerivative2, ())
end

function sin_cos(obj::FieldUnivariateDerivative2)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "sinh", FieldUnivariateDerivative2, ())
end

function sqrt(obj::FieldUnivariateDerivative2)
    return jcall(obj, "sqrt", FieldUnivariateDerivative2, ())
end

function subtract(obj::FieldUnivariateDerivative2, arg0::FieldUnivariateDerivative2)
    return jcall(obj, "subtract", FieldUnivariateDerivative2, (FieldUnivariateDerivative2,), arg0)
end

function subtract(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "subtract", FieldUnivariateDerivative2, (RealFieldElement,), arg0)
end

function subtract(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "subtract", FieldUnivariateDerivative2, (jdouble,), arg0)
end

function tan(obj::FieldUnivariateDerivative2)
    return jcall(obj, "tan", FieldUnivariateDerivative2, ())
end

function tanh(obj::FieldUnivariateDerivative2)
    return jcall(obj, "tanh", FieldUnivariateDerivative2, ())
end

function taylor(obj::FieldUnivariateDerivative2, arg0::RealFieldElement)
    return jcall(obj, "taylor", RealFieldElement, (RealFieldElement,), arg0)
end

function taylor(obj::FieldUnivariateDerivative2, arg0::jdouble)
    return jcall(obj, "taylor", RealFieldElement, (jdouble,), arg0)
end

function to_degrees(obj::FieldUnivariateDerivative2)
    return jcall(obj, "toDegrees", FieldUnivariateDerivative2, ())
end

function to_derivative_structure(obj::FieldUnivariateDerivative2)
    return jcall(obj, "toDerivativeStructure", FieldDerivativeStructure, ())
end

function to_radians(obj::FieldUnivariateDerivative2)
    return jcall(obj, "toRadians", FieldUnivariateDerivative2, ())
end


text:
# Import the functions for which we define methods to be called when using PeriodicNumbers
import Base.show
import Base: +, -, *, /, abs, abs2
using LinearAlgebra

#######################################################################################################################
#######################################################################################################################
### Definition of 'PeriodicNumber' with its respective methods for basic operations
#######################################################################################################################
#######################################################################################################################
"""
    PeriodicNumber{T}(value::T, L::T)

Define a number whose value is always contained between 0 and L; more precisely ∈ [0,L).

When it is initialized with value>L or value<0 the value is transformed to be its L-modulus.

"""
struct PeriodicNumber{T<:Real} <:Number
    value::T # Value to be converted if not ∈ [0,L)
    L::T # upper bound
    function PeriodicNumber{T}(value, L) where {T<:Real}
        zero_T = convert(T, 0)
        if L<=zero_T
            error("The period should be a positive number.")
        end
        if (value<zero_T)||(value>=L)
            value = mod(value, L)
            if value==L
                value=zero_T
            end
        end
    new(value, L)
    end
end

# Simple constructor for periodic number
"Constructor for 'PeriodicNumber', with L=1 by default."
PeriodicNumber(x::T, L::T=1.0) where {T<:Real} = PeriodicNumber{T}(x, L)

#### Format the output of periodic numbers
function Base.show(io::IO, PN::PeriodicNumber)
    L = PN.L
    if L==1.0
        print(io, "PN(", PN.value, ")" )
    else
        print(io, "PN(", PN.value, " L=", L, ")" )
    end
end
PeriodicNumber(1.2)

#################
### Definition of basic operations
#################

# Define sum of periodic values among themselves or with real numbers. The result is always a PeriodicNumber type
function Base.:+(PN1::PeriodicNumber, PN2::PeriodicNumber)::PeriodicNumber
    x1=PN1.value; x2=PN2.value
    L1=PN1.L; L2=PN2.L
    if (L1 !=L2)
        error("Comparing PeriodicNumbers of different periodicity")
    end

    return PeriodicNumber(x1+x2, L1)
end

function Base.:+(PN1::PeriodicNumber, x2::Real)::PeriodicNumber
    x1=PN1.value;  L=PN1.L
    return PeriodicNumber(x1+x2, L)
end

Base.:+(x1::Real, PN2::PeriodicNumber)::PeriodicNumber = PN2+x1

# Define product periodic numbers among themselves or with real numbers. The result is always a PeriodicNumber type
function Base.:*(PN1::PeriodicNumber, PN2::PeriodicNumber)::PeriodicNumber
    x1=PN1.value; x2=PN2.value
    L1=PN1.L; L2=PN2.L
    if (L1 !=L2)
        error("Comparing quantities of different periodicity")
    end
    return PeriodicNumber(x1*x2, L1)
end

function Base.:*(PN1::PeriodicNumber, x2::Real)::PeriodicNumber
    x1=PN1.value; L=PN1.L
    return PeriodicNumber(x1*x2, L)
end

Base.:*(x1::Real, PN2::PeriodicNumber)::PeriodicNumber = PN2*x1

#= Define difference of periodic values among themselves or with real numbers. The result is always a PeriodicNumber type, and hence positive.
The difference between periodic values is somewhat tricky, because we need to consider that two particles far apart in the real space might be very close in the periodic one.
For instance, if x1=0 and x2=L, the real distance is L, but in the periodic space it's 0
# NB: The following function is essentially computing a distance in the periodic space, since the value is always positive and
does not depend on the order, i.e. x1-x2 = x2-x1.
=#
"""
    -(PN1::PeriodicNumber, PN2::PeriodicNumber)

Return the difference of two periodic numbers, also as a 'PeriodicNumber' type.

Given that the output is a periodic number, and hence its value is positive, this 'difference
function' is actually commutative.
Moreover, the output's value can never be larger than L/2.
"""
function Base.:-(PN1::PeriodicNumber, PN2::PeriodicNumber)::PeriodicNumber
    x1=PN1.value; x2=PN2.value
    L1=PN1.L; L2=PN2.L
    if (L1 !=L2)
        error("Comparing quantities of different periodicity")
    end
    arg=x1-x2
    arg=min(PeriodicNumber(arg, L1).value, PeriodicNumber(-arg, L1).value )
    return PeriodicNumber(arg, L1)
end

#=
NB: This is also an UNsigned difference.
Thus, it is made to work only with PeriodicNumber's whose lower limit is equal to 0
=#
function Base.:-(PN1::PeriodicNumber, x2::Real)::PeriodicNumber
    x1=PN1.value; L=PN1.L
    dif=abs(x1-x2)
    while dif>0.5*L
        dif -= L
    end
    return PeriodicNumber(dif, L)
end

#NB: Note also that since the difference of a periodic number and any real number is UNsigned, the '-' operator is commutative
Base.:-(x2::Real, PN1::PeriodicNumber)::PeriodicNumber = PN1-x2

#Simple first calls for compilation purposes
PeriodicNumber(0.1, 3.0)+ PeriodicNumber(2.5, 3.0); PeriodicNumber(2.5, 3.0) - PeriodicNumber(0.1, 3.0)
PeriodicNumber(0.1, 3.0) +2.0; PeriodicNumber(0.1, 3.0) - 2.0
rand()*PeriodicNumber(2.5, 3.0); PeriodicNumber(2.5, 3.0)*rand()

# Define other useful methods for computing absolute values, norm, etc.
value(PN::PeriodicNumber{<:Real}) = PN.value
abs(PN::PeriodicNumber{<:Real}) = Base.abs(value(PN))
abs2(PN::PeriodicNumber{<:Real}) = Base.abs2(value(PN))

# Define a 'distance' function to work with numbers with periodic BC
"""
    periodic_distance(PN1::PeriodicNumber, PN2::PeriodicNumber)

Compute the *un*-signed distance between a pair of periodic numbers.

It is implemented simply as the 'value' field of '-' method for periodic numbers.
Since such method always returns a non-negative value and is commutative, this function
is a well defined distance.
"""
periodic_distance(P1::PeriodicNumber, P2::PeriodicNumber)::Float64=(P1-P2).value

ext: jl

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
using Printf: @sprintf

push!(LOAD_PATH, @__DIR__)
using PCREO

const IMAGE_URL_PREFIX = "http://web.stanford.edu/~dkzhang/PCREOImages/"
const DATABASE_URL_PREFIX = "http://web.stanford.edu/~dkzhang/PCREODatabase/"
const HTML_OUTPUT_DIRECTORY = "C:\\Users\\Zhang\\Documents\\GitHub\\my-madoko-docs\\PersonalWebsite\\PCREOTables"

function point_group_string(points)
    group_name = point_group(points)
    if group_name == "C_1"
        return ""
    elseif '_' in group_name
        a, b = split(group_name, '_')
        return "$a$b"
    end
    return group_name
end

function sigfigs_string(x::Float64)
    @assert x >= 0.0
    if x == 0.0
        return "0.000000000000"
    elseif x < 10.0
        return @sprintf("%.12f", x)
    elseif x < 100.0
        return @sprintf("%.11f", x)
    elseif x < 1000.0
        return @sprintf("%.10f", x)
    elseif x < 10000.0
        return @sprintf("%.9f", x)
    elseif x < 100000.0
        return @sprintf("%.8f", x)
    elseif x < 1000000.0
        return @sprintf("%.7f", x)
    elseif x < 10000000.0
        return @sprintf("%.6f", x)
    elseif x < 100000000.0
        return @sprintf("%.5f", x)
    elseif x < 1000000000.0
        return @sprintf("%.4f", x)
    elseif x < 10000000000.0
        return @sprintf("%.3f", x)
    elseif x < 100000000000.0
        return @sprintf("%.2f", x)
    elseif x < 1000000000000.0
        return @sprintf("%.1f", x)
    else
        @assert false
    end
end

function percentage_string(x::Float64)
    @assert 0.0 <= x <= 100.0
    if x == 0.0
        return "0.000%"
    elseif x == 100.0
        return "100.0%"
    elseif x < 10.0
        return @sprintf("%.3f", x) * '%'
    else
        return @sprintf("%.2f", x) * '%'
    end
end

const TABLE_HEADER = """

      Image
      Z-score
      Energy
      Sym.
      Frequency
      PCREO Database ID

"""

const TABLE_FOOTER = """

"""

function main()
    for num_points = 50 : 99
        ids = filter(startswith("PCREO-03-" * lpad(num_points, 4, '0')),
                     readdir(PCREO_DATABASE_DIRECTORY))
        @assert !isempty(ids)
        records = Dict(id => PCREORecord(id) for id in ids)
        energies = sort!([
            (records[id].energy, id, length(readdir(
                joinpath(PCREO_DATABASE_DIRECTORY, id))))
            for id in ids])
        total_count = sum(n for (_, _, n) in energies)
        if length(energies) > 1
            total_energy = sum(n * e for (e, _, n) in energies)
            mean_energy = total_energy / total_count
            total_deviation = sum(n * (e - mean_energy)^2
                                  for (e, _, n) in energies)
            stdev_energy = sqrt(total_deviation / (total_count - 1))
        else
            mean_energy = first(energies)[1]
            stdev_energy = 1.0
        end
        table_entries = String[]
        for (energy, id, freq) in energies
            z_score = (energy - mean_energy) / stdev_energy
            z_score_string = @sprintf("%+.3f", z_score)
            defects = defect_classes(records[id].facets)
            defect_strings = [
                string(num) * " × " * html_defect_string(shapes)
                for (num, shapes) in defects]
            push!(table_entries, """

                  $z_score_string
                  $(sigfigs_string(energy))
                  $(point_group_string(records[id].points))
                  $freq ($(percentage_string(100 * freq / total_count)))

                    $(id[15:end])

                  $(join(defect_strings, " "))

            """)
        end
        open(joinpath(HTML_OUTPUT_DIRECTORY, "$num_points.html"), "w+") do io
            println(io, "
$num_points Points
")
            print(io, TABLE_HEADER)
            for entry in table_entries
                print(io, entry)
            end
            print(io, TABLE_FOOTER)
        end
    end
end

main()


text:
module Utils

    import ..Lib

    const MAPPING = Lib.CEnum.namemap(Lib.helics_error_types)

    function snakecase_to_camelcase(s::String)
        s = replace(s, r"['`]"=>"")
        s = split(s, r"[a-z]\K(?=[A-Z][a-z])|\W|_")
        s = join(map(uppercasefirst, s))
    end

    function snakecase_to_camelcase(s::Symbol)
        s = String(s)
        s = snakecase_to_camelcase(s)
        s = Symbol(s)
    end

    for (sym, ans) in Lib.CEnum.name_value_pairs(Lib.helics_error_types)
        sym = snakecase_to_camelcase(sym)
        eval(:(struct $sym <: Exception
                   msg::String
               end))
    end

    function invoke_and_check(f::Function, args...)
        err = Lib.helicsErrorInitialize()
        r = f(args..., Ref(err))
        if err.error_code == 0
            return r
        elseif err.error_code == -8
            func = String(nameof(f))
            c = Lib.helics_error_types(err.error_code)
            @warn "`$func` returned warning `$(unsafe_string(err.message))`"
            return r
        else
            func = String(nameof(f))
            s = snakecase_to_camelcase(MAPPING[err.error_code])
            S = getfield(Utils, Symbol(s))
            throw(S("`$func` errored with message `$(unsafe_string(err.message))`"))
        end
    end

    macro invoke_and_check(expr)
        pushfirst!(expr.args, :(Utils.invoke_and_check))
        return esc(expr)
    end

end


text:
import Combinatorics: combinations

## ---------- DISCRIMINANT MEASURES ----------
"""
Discriminant measure for Local Discriminant Basis. Current available subtypes
are:
- [`ProbabilityDensityDM`](@ref)
- [`SignaturesDM`](@ref)
"""
abstract type DiscriminantMeasure end

"""
Discriminant measure for Probability Density and Time Frequency based energy
maps. Current available measures are:
- [`AsymmetricRelativeEntropy`](@ref)
- [`SymmetricRelativeEntropy`](@ref)
- [`LpDistance`](@ref)
- [`HellingerDistance`](@ref)
"""
abstract type ProbabilityDensityDM <: DiscriminantMeasure end

"""
Discriminant measure for Signatures based energy maps. Current available
measures are:
- [`EarthMoverDistance`](@ref)
"""
abstract type SignaturesDM <: DiscriminantMeasure end

@doc raw"""
    AsymmetricRelativeEntropy <: ProbabilityDensityDM

Asymmetric Relative Entropy discriminant measure for the Probability Density and Time
Frequency based energy maps. This measure is also known as cross entropy and
Kullback-Leibler divergence.

Equation:

``D(p,q) = \sum p(x) \log \frac{p(x)}{q(x)}``

where ``\int_{-\infty}^{\infty} p(t) dt = \int_{-\infty}^{\infty} q(t) dt = 1`` (``p(t)``
and ``q(t)`` are probability density functions.)
"""
struct AsymmetricRelativeEntropy <: ProbabilityDensityDM end

@doc raw"""
    SymmetricRelativeEntropy <: ProbabilityDensityDM

Symmetric Relative Entropy discriminant measure for the Probability Density and Time
Frequency energy maps. Similar idea to the Asymmetric Relative Entropy, but this aims to
make the measure more symmetric.

Equation: Denote the Asymmetric Relative Entropy as ``D_A(p,q)``, then

``D(p,q) = D_A(p,q) + D_A(q,p) = \sum p(x) \log \frac{p(x)}{q(x)} + q(x) \log
\frac{q(x)}{p(x)}``

where ``\int_{-\infty}^{\infty} p(t) dt = \int_{-\infty}^{\infty} q(t) dt = 1`` (``p(t)``
and ``q(t)`` are probability density functions.)

**See also:** [`AsymmetricRelativeEntropy`](@ref)
"""
struct SymmetricRelativeEntropy <: ProbabilityDensityDM end

@doc raw"""
    LpDistance <: ProbabilityDensityDM

``\ell^p`` Distance discriminant measure for the Probability Density and Time
Frequency based energy maps. The default ``p`` value is set to 2.

Equation:

``W(q,r) = ||q-r||_p^p = \sum_{i=1}^n (q_i - r_i)^p``
"""
@with_kw struct LpDistance{T<:Real} <: ProbabilityDensityDM
    p::T = 2
end

@doc raw"""
    HellingerDistance <: ProbabilityDensityDM

Hellinger Distance discriminant measure for the Probability Density energy
map.

Equation:

``H(p,q) = \sum_{i=1}^n (\sqrt{p_i} - \sqrt{q_i})^2``
"""
struct HellingerDistance <: ProbabilityDensityDM end

@doc raw"""
    EarthMoverDistance <: SignaturesDM

Earth Mover Distance discriminant measure for the Signatures energy map.

Equation:

``E(P,Q) = \frac{\sum_{k=1}^{m+n+1} |\hat p_k - \hat q_k| (r_{k+1} - r_k)}{w_\Sigma}``

where
- ``r_1, r_2, \ldots, r_{m+n}`` is the sorted list of ``p_1, \ldots, p_m, q_1, \ldots, q_n``
- ``\hat p_k = \sum_{p_i \leq r_k} w_{p_i}``
- ``\hat q_k = \sum_{q_i \leq r_k} w_{q_i}``
"""
struct EarthMoverDistance <: SignaturesDM end

"""
    discriminant_measure(Γ, dm)

Computes the discriminant measure of each subspace calculated from the energy maps.

# Arguments
- `Γ`: Energy map computed from `energy_map` function. The data structures of `Γ`, depending
  on their corresponding energy map, should be:
    - `TimeFrequency()`: `AbstractArray{T,3}` for 1D signals or `AbstractArray{T,4}` for 2D
      signals.
    - `ProbabilityDensity()`: `AbstractArray{T,4}` for 1D signals or `AbstractArray{T,5}`
      for 2D signals.
    - `Signatures()`: `AbstractVector{NamedTuple{(:coef, :weight), Tuple{S₁,S₂}}}`.
- `dm::DiscriminantMeasure`: Type of Discriminant Measure. The type of `dm` must match the
  type of `Γ`.

# Returns
- `D::Array{T}`: Discriminant measure at each coefficient of the decomposed signals.

# Examples
```julia
using Wavelets, WaveletsExt

X, y = generateclassdata(ClassData(:tri, 5, 5, 5))
Xw = wpdall(X, wavelet(WT.haar))

Γ = energy_map(Xw, y, TimeFrequency()); discriminant_measure(Γ, AsymmetricRelativeEntropy())
Γ = energy_map(Xw, y, ProbabilityDensity()); discriminant_measure(Γ, LpDistance())
Γ = energy_map(Xw, y, Signatures()); discriminant_measure(Γ, EarthMoverDistance())
```

**See also:** [`pairwise_discriminant_measure`](@ref)
"""
function discriminant_measure(Γ::AbstractArray{T}, dm::ProbabilityDensityDM) where
                              T<:AbstractFloat
    # Basic summary of data
    N = ndims(Γ)
    nc = size(Γ,N)
    @assert 3 ≤ N ≤ 5
    @assert nc > 1
    #=======================================================================================
    Classifying the type of signals and energy maps based on Γ

    1. ndims(Γ)=3                           : Time frequency energy map on 1D signals
    2. ndims(Γ)=4 and size(Γ,3) small (<100): Time frequency energy map on 2D signals
    3. ndims(Γ)=4 and size(Γ,3) large (≥100): Probability density energy map on 1D signals
    4. ndims(Γ)=5                           : Probability density energy map on 2D signals
    =======================================================================================#
    if (N==3) || (N==4 && size(Γ,3)≥100)
        sz = [size(Γ,1)]
        L = size(Γ,2)
        pdf_len = N==3 ? 1 : size(Γ,3)
    elseif (N==4 && size(Γ,3)<100) || (N==5)
        sz = size(Γ)[1:2]
        L = size(Γ,3)
        pdf_len = N==5 ? size(Γ,4) : 1
    else
        throw(ArgumentError("Γ has uninterpretable dimensions/unknown size."))
    end

    # --- Compute pairwise discriminant measure ---
    D = zeros(T, (sz...,L))
    for (i,j) in combinations(1:nc,2)
        if N==3
            Γᵢ = @view Γ[:,:,i]         # Energy map for i-th class
            Γⱼ = @view Γ[:,:,j]         # Energy map for j-th class
        elseif N==4
            Γᵢ = @view Γ[:,:,:,i]       # Energy map for i-th class
            Γⱼ = @view Γ[:,:,:,j]       # Energy map for j-th class
        else # N==5
            Γᵢ = @view Γ[:,:,:,:,i]     # Energy map for i-th class
            Γⱼ = @view Γ[:,:,:,:,j]     # Energy map for j-th class
        end
        D += pairwise_discriminant_measure(Γᵢ, Γⱼ, dm)
    end
    return D
end

# discriminant measure for EMD
function discriminant_measure(Γ::AbstractVector{NamedTuple{(:coef, :weight), Tuple{S₁,S₂}}},
                              dm::SignaturesDM) where
                             {S₁<:Array{T} where T<:AbstractFloat,
                              S₂<:Union{T,Array{T}} where T<:AbstractFloat}
    # Basic summary of data
    nc = length(Γ)
    @assert nc > 1
    T = eltype(Γ[1][:coef])
    sz = size(Γ[1][:coef])[1:end-2]
    L = size(Γ[1][:coef])[end-1]

    D = zeros(T, (sz...,L))
    for (Γᵢ,Γⱼ) in combinations(Γ,2)
        D += pairwise_discriminant_measure(Γᵢ, Γⱼ, dm)
    end
    return D
end

# discriminant measure between 2 energy maps
"""
    pairwise_discriminant_measure(Γ₁, Γ₂, dm)

Computes the discriminant measure between 2 classes based on their energy maps.

# Arguments
- `Γ₁::AbstractArray{T} where T<:AbstractFloat`: Energy map for class 1.
- `Γ₂::AbstractArray{T} where T<:AbstractFloat`: Energy map for class 2.
- `dm::DiscriminantMeasure`: Type of discriminant measure.

# Returns
- `::Array{T}`: Discriminant measure between `Γ₁` and `Γ₂`.
"""
function pairwise_discriminant_measure(Γ₁::AbstractArray{T}, Γ₂::AbstractArray{T},
                                       dm::ProbabilityDensityDM) where T<:AbstractFloat
    # parameter checking and basic summary
    N = ndims(Γ₁)
    @assert 2 ≤ ndims(Γ₁) ≤ 4
    @assert size(Γ₁) == size(Γ₂)
    #=======================================================================================
    Classifying the type of signals and energy maps based on Γ

    1. ndims(Γ)=2                           : Time frequency energy map on 1D signals
    2. ndims(Γ)=3 and size(Γ,3) small (<100): Time frequency energy map on 2D signals
    3. ndims(Γ)=3 and size(Γ,3) large (≥100): Probability density energy map on 1D signals
    4. ndims(Γ)=4                           : Probability density energy map on 2D signals
    =======================================================================================#
    if (N==2) || (N==3 && size(Γ₁,3)≥100)
        sz = [size(Γ₁,1)]
        L = size(Γ₁,2)
    elseif (N==3 && size(Γ₁,3)<100) || (N==4)
        sz = size(Γ₁)[1:2]
        L = size(Γ₁,3)
    else
        throw(ArgumentError("Γ has uninterpretable dimensions/unknown size."))
    end

    # --- Pairwise discriminant measure for each element ---
    D = zeros(T, (sz...,L))
    slice_size = prod([sz...,L])        # Number of elements in each slice of the pdf
    map_size = prod(size(Γ₁))           # Number of elements in entire energy map
    for i in 1:slice_size
        for j in i:slice_size:map_size
            D[i] += pairwise_discriminant_measure(Γ₁[j], Γ₂[j], dm)
        end
    end
    return D
end

# discriminant measure between 2 nergy maps for EMD
function pairwise_discriminant_measure(Γ₁::NamedTuple{(:coef, :weight), Tuple{S₁,S₂}},
                                       Γ₂::NamedTuple{(:coef, :weight), Tuple{S₁,S₂}},
                                       dm::SignaturesDM) where
                                      {S₁<:Array{T} where T<:AbstractFloat,
                                       S₂<:Union{T,Array{T}} where T<:AbstractFloat}

    # parameter checking and basic summary
    N = ndims(Γ₁.coef)
    T = eltype(Γ₁.coef)
    sz = size(Γ₁.coef)[1:end-2]
    L = size(Γ₁.coef)[end-1]
    @assert isa(Γ₁.weight, AbstractFloat) || size(Γ₁.coef) == size(Γ₁.weight)
    @assert isa(Γ₂.weight, AbstractFloat) || size(Γ₂.coef) == size(Γ₂.weight)
    @assert typeof(Γ₁.weight) == typeof(Γ₂.weight)
    @assert sz == size(Γ₂.coef)[1:end-2]
    @assert L == size(Γ₂.coef)[end-1]

    D = Array{T,N-1}(undef, (sz...,L))
    slice_size = prod([sz...,L])                # Number of elements for each signal's coefficients
    for i in 1:slice_size
        # Signatures
        if isa(Γ₁.weight, AbstractFloat)        # Equal weight
            P = (coef = Γ₁.coef[i:slice_size:end], weight = Γ₁.weight)
            Q = (coef = Γ₂.coef[i:slice_size:end], weight = Γ₂.weight)
        else                                    # Probability density weight
            P = (coef = Γ₁.coef[i:slice_size:end], weight = Γ₁.weight[i:slice_size:end])
            Q = (coef = Γ₂.coef[i:slice_size:end], weight = Γ₂.weight[i:slice_size:end])
        end
        D[i] = pairwise_discriminant_measure(P, Q, dm)
    end
    return D
end

# Asymmetric Relative Entropy
"""
    pairwise_discriminant_measure(p, q, dm)

Computes the discriminant measure between 2 classes at an index ``i``, where `Γ₁[i] = p` and
`Γ₂[i] = q`.

# Arguments
- `p::T where T<:AbstractFloat`: Coefficient at index ``i`` from `Γ₁`.
- `q::T where T<:AbstractFloat`: Coefficient at index ``i`` from `Γ₂`.
- `dm::DiscriminantMeasure`: Type of discriminant measure.

# Returns
- `::T`: Discriminant measure between `p` and `q`.
"""
function pairwise_discriminant_measure(p::T, q::T, dm::AsymmetricRelativeEntropy) where
                                       T<:AbstractFloat
    @assert p ≥ 0 && q ≥ 0
    return (p==0 || q==0) ? 0 : p * log(p/q)
end

# Symmetric Relative Entropy
function pairwise_discriminant_measure(p::T, q::T, dm::SymmetricRelativeEntropy) where
                                       T<:AbstractFloat
    return pairwise_discriminant_measure(p, q, AsymmetricRelativeEntropy()) +
           pairwise_discriminant_measure(q, p, AsymmetricRelativeEntropy())
end

# Hellinger Distance
function pairwise_discriminant_measure(p::T, q::T, dm::HellingerDistance) where
                                       T<:AbstractFloat
    return (sqrt(p) - sqrt(q))^2
end

# Lᵖ Distance
function pairwise_discriminant_measure(p::T, q::T, dm::LpDistance) where T<:AbstractFloat
    return (p - q)^dm.p
end

# Earth Mover Distance
function pairwise_discriminant_measure(P::NamedTuple{(:coef, :weight), Tuple{S1,S2}},
                                       Q::NamedTuple{(:coef, :weight), Tuple{S1, S2}},
                                       dm::EarthMoverDistance) where
                                      {S1<:Vector{T} where T<:AbstractFloat,
                                       S2<:Union{T,Vector{T}} where T<:AbstractFloat}
    # assigning tuple signatures into coef and weight
    p, w_p = P
    q, w_q = Q

    # sort signature values
    p_order = sortperm(p)
    q_order = sortperm(q)
    p = p[p_order]
    q = q[q_order]
    w_p = isa(w_p, AbstractFloat) ? repeat([w_p], length(p)) : w_p[p_order]
    w_q = isa(w_q, AbstractFloat) ? repeat([w_q], length(q)) : w_q[q_order]

    # concatenate p and q, then sort them
    r = [p; q]
    sort!(r)

    # compute emd
    n = length(r)
    emd = 0
    for i in 1:(n-1)
        # get total weight of p and q less than or equal to r[i]
        ∑w_p = sum(w_p[p .≤ r[i]])
        ∑w_q = sum(w_q[q .≤ r[i]])
        # add into emd
        @inbounds emd += abs(∑w_p - ∑w_q) * (r[i+1] - r[i])
    end
    emd /= (sum(w_p) + sum(w_q))
    return emd
end

## ---------- DISCRIMINATION POWER ----------
"""
Discriminant Power measure for the Local Discriminant Basis. Current available
measures are
- [`BasisDiscriminantMeasure`](@ref)
- [`FishersClassSeparability`](@ref)
- [`RobustFishersClassSeparability`](@ref)
"""
abstract type DiscriminantPower end

"""
    BasisDiscriminantMeasure <: DiscriminantPower

This is the discriminant measure of a single basis function computed in a previous step to
construct the energy maps.
"""
struct BasisDiscriminantMeasure <: DiscriminantPower end

@doc raw"""
    FishersClassSeparability <: DiscriminantPower

The Fisher's class separability of the expansion coefficients in the basis function.

Equation: ``\frac{\sum_{c=1}^C \pi_c\{{\rm mean}_i(\alpha_{\lambda,i}^{(c)}) - {\rm mean}_c
\cdot {\rm mean}_i(\alpha_{\lambda,i}^{(c)})\}^2}{\sum_{c=1}^C \pi_c {\rm
var}_i(\alpha_{\lambda,i}^{(c)})}``
"""
struct FishersClassSeparability <: DiscriminantPower end

@doc raw"""
    RobustFishersClassSeparability <: DiscriminantPower

The robust version of Fisher's class separability of the expansion coefficients in the basis
function.

Equation: ``\frac{\sum_{c=1}^C \pi_c\{{\rm med}_i(\alpha_{\lambda,i}^{(c)}) - {\rm med}_c
\cdot {\rm med}_i(\alpha_{\lambda,i}^{(c)})\}^2}{\sum_{c=1}^C \pi_c {\rm
mad}_i(\alpha_{\lambda,i}^{(c)})}``
"""
struct RobustFishersClassSeparability <: DiscriminantPower end

"""
    discriminant_power(D, tree, dp)
    discriminant_power(coefs, y, dp)

Returns the discriminant power of each leaf from the local discriminant basis (LDB) tree.

# Arguments
- `D::AbstractArray{T} where T<:AbstractFloat`: Discriminant measures.
- `tree::BitVector`: Best basis tree for selecting coefficients with largest discriminant
  measures.
- `coefs::AbstractArray{T} where T<:AbstractFloat`: Best basis coefficients for the input
  signals.
- `y::AbstractVector{S} where S`: Labels corresponding to each signal in `coefs`.
- `dp::DiscriminantPower`: The measure for discriminant power.

!!! note
    `discriminant_power(D, tree, dp)` only works for `dp = BasisDiscriminantMeasure()`,
    whereas `discriminant_power(coefs, y, dp)` works for `dp = FishersClassSeparability()`
    and `dp = RobustFishersClassSeparability()`.

# Returns
- `power::Array{T}`: The discriminant power at each index of `D` or `coefs`.
- `order::Vector{T}`: The order of discriminant power in descending order.
"""
function discriminant_power(D::AbstractArray{T}, tree::BitVector,
                            dp::BasisDiscriminantMeasure) where T<:AbstractFloat
    @assert 2 ≤ ndims(D) ≤ 3
    N = ndims(D)
    sz = size(D)[1:end-1]
    tmp = Array{T,N-1}(undef, sz)
    @assert isvalidtree(tmp, tree)

    power = getbasiscoef(D, tree)
    order = sortperm(vec(power), rev = true)

    return (power, order)
end

function discriminant_power(coefs::AbstractArray{T}, y::AbstractVector{S},
                            dp::FishersClassSeparability) where {T<:AbstractFloat, S}
    N = ndims(coefs)
    @assert 2 ≤ N ≤ 3
    sz = size(coefs)[1:end-1]           # signal length

    C = unique(y)
    Nc = length(C)                      # number of classes

    Nᵢ = Array{T,1}(undef, Nc)
    Eαᵢ = Array{T,N}(undef, (sz...,Nc))                # mean of each entry
    Varαᵢ = Array{T,N}(undef, (sz...,Nc))              # variance of each entry
    for (i, c) in enumerate(C)
        idx = findall(yᵢ -> yᵢ==c, y)
        Nᵢ[i] = length(idx)
        if N==2
            coefsᵢ = coefs[:,idx]
            Eαᵢ[:,i] = mean(coefsᵢ, dims = N)
            Varαᵢ[:,i] = var(coefsᵢ, dims = N)
        elseif N==3
            coefsᵢ = coefs[:,:,idx]
            Eαᵢ[:,:,i] = mean(coefsᵢ, dims = N)
            Varαᵢ[:,:,i] = var(coefsᵢ, dims = N)
        end
    end
    Eα = mean(Eαᵢ, dims = N)                      # overall mean of each entry
    pᵢ = Nᵢ / sum(Nᵢ)                             # proportions of each class

    if N==2     # For 1D signals, can be done via matrix multiplication
        power = ((Eαᵢ - (Eα .* Eαᵢ)).^2 * pᵢ) ./ (Varαᵢ * pᵢ)
    else        # For 2D signals, requires some manual work
        pᵢ = reshape(pᵢ,1,1,:)
        power = sum((Eαᵢ-(Eα.*Eαᵢ)).^2 .* pᵢ, dims=N) ./ sum(Varαᵢ.*pᵢ, dims=N)
        power = reshape(power, sz...)
    end
    order = sortperm(vec(power), rev = true)

    return (power, order)
end

function discriminant_power(coefs::AbstractArray{T}, y::AbstractVector{S},
                            dp::RobustFishersClassSeparability) where {T<:AbstractFloat,S}
    N = ndims(coefs)
    @assert 2 ≤ N ≤ 3
    sz = size(coefs)[1:end-1]           # signal length

    C = unique(y)
    Nc = length(C)                      # number of classes

    Nᵢ = Array{T,1}(undef, Nc)
    Medαᵢ = Array{T,N}(undef, (sz...,Nc))              # mean of each entry
    Madαᵢ = Array{T,N}(undef, (sz...,Nc))              # variance of each entry
    for (i, c) in enumerate(C)
        idx = findall(yᵢ -> yᵢ==c, y)
        Nᵢ[i] = length(idx)
        if N==2
            coefsᵢ = coefs[:,idx]
            Medαᵢ[:,i] = median(coefsᵢ, dims = N)
            Madαᵢ[:,i] = mapslices(x -> mad(x, normalize=false), coefsᵢ, dims = N)
        elseif N==3
            coefsᵢ = coefs[:,:,idx]
            Medαᵢ[:,:,i] = median(coefsᵢ, dims = N)
            Madαᵢ[:,:,i] = mapslices(x -> mad(x, normalize=false), coefsᵢ, dims = N)
        end
    end
    Medα = median(Medαᵢ, dims = N)               # overall mean of each entry
    pᵢ = Nᵢ / sum(Nᵢ)                            # proportions of each class

    if N==2     # For 1D signals, can be done via matrix multiplication
        power = ((Medαᵢ - (Medα.*Medαᵢ)).^2 *pᵢ) ./ (Madαᵢ * pᵢ)
    else        # For 2D signals, requires some manual work
        pᵢ = reshape(pᵢ,1,1,:)
        power = sum((Medαᵢ-(Medα.*Medαᵢ)).^2 .* pᵢ, dims=N) ./ sum(Madαᵢ.*pᵢ, dims=N)
        power = reshape(power, sz...)
    end
    order = sortperm(vec(power), rev = true)

    return (power, order)
end


text:
module QJuliaFields

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "../..", "core"))
push!(LOAD_PATH, @__DIR__)

import QJuliaRegisters
import QJuliaEnums
import QJuliaGrid

const complex_length = 2

abstract type QJuliaGenericField_qj end

mutable struct QJuliaLatticeField_qj <: QJuliaGenericField_qj
  # Field Matrix{T}, e.g., a single color spinor, a block color spinor (an eigenvector set), an SU(N) field etc.:
  v::Any

  # Lattice field structure
  field_desc::QJuliaGrid.QJuliaLatticeFieldDescr_qj

  QJuliaLatticeField_qj() = new( missing )

  function QJuliaLatticeField_qj(fdesc::QJuliaGrid.QJuliaLatticeFieldDescr_qj)

    # Since the field descriptor is an external argument (and data are not protected in Julia) we need this check again:
    if fdesc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY && fdesc.nSpin == nothing ; error("Spin dof is not defined for the field with scalar geometry."); end

    if fdesc.grid.location == QJuliaEnums.QJULIA_CPU_FIELD_LOCATION
      # Set the field array total elements:
      tot_elems = complex_length*fdesc.siteSubset*(fdesc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY ? fdesc.nSpin : fdesc.nColor)*fdesc.nColor*fdesc.grid.grid_volumeCB

      v = Matrix{fdesc.grid.register_type}(undef, tot_elems, Int(fdesc.geom)*fdesc.nBlock)
    else
      # Ok, set the field as a "missing" one. It will be inititalized later.
      v = missing
    end

    # call constructor
    new(v, fdesc)

  end #QJuliaLatticeField_qj constructor

  function QJuliaLatticeField_qj(grid::QJuliaGrid.QJuliaGridDescr_qj, geom::QJuliaEnums.QJuliaFieldGeometry_qj, parity::QJuliaEnums.QJuliaParity_qj; NSpin = nothing, NColor = 3, NBlock = 1)

    fdesc = QJuliaLatticeFieldDescr_qj{NSpin, NColor, NBlock}(grid::QJuliaGridDescr_qj, geom::QJuliaEnums.QJuliaFieldGeometry_qj, parity::QJuliaEnums.QJuliaParity_qj)

    if fdesc.grid.location == QJuliaEnums.QJULIA_CPU_FIELD_LOCATION
      # Set the field array total elements:
      tot_elems = complex_length*fdesc.siteSubset*(fdesc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY ? fdesc.nSpin : fdesc.nColor)*fdesc.nColor*fdesc.grid.grid_volumeCB

      v = Matrix{fdesc.grid.register_type}(undef, tot_elems, Int(fdesc.geom)*fdesc.nBlock)
    else
      # Ok, set the field as a "missing" one. It will be inititalized later.
      v = missing
    end

    # call constructor
    new(v, fdesc)

  end #QJuliaLatticeField_qj constructor

  function QJuliaLatticeField_qj(field::QJuliaLatticeField_qj, parity::QJuliaEnums.QJuliaParity_qj)

     if field.v === missing; error("Cannot reference a missing field"); end
     if field.field_desc.siteSubset == 1; error("Cannot reference a parity field from the parity field"); end

     fdesc = deepcopy(field.field_desc)

     griddsc   = fdesc.grid
     griddsc.X = NTuple{QJuliaGrid.QJULIA_MAX_DIMS, Int}(ntuple(i->(i == 1 ? field.field_desc.grid.X[i] >> 1 : field.field_desc.grid.X[i]), QJuliaGrid.QJULIA_MAX_DIMS))

     griddsc.volume        = field.field_desc.grid.volumeCB
     griddsc.volumeCB      = field.field_desc.grid.volumeCB
     griddsc.grid_volume   = field.field_desc.grid.grid_volumeCB
     griddsc.grid_volumeCB = field.field_desc.grid.grid_volumeCB

     fdesc.parity     = parity
     fdesc.siteSubset = 1

     parity_volume = complex_length*griddsc.grid_volumeCB*(fdesc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY ? fdesc.nSpin : fdesc.nColor)*fdesc.nColor

     offset = parity == QJuliaEnums.QJULIA_EVEN_PARITY ? 1 : parity_volume+1

     v = view(field.v, offset:offset+parity_volume-1, :)

     new(v, fdesc)

  end

end #QJuliaLatticeField_qj

# Pointer helpers:
@inline function GetPtr(field::QJuliaGenericField_qj); return field.v; end

@inline function GetParityPtr(field::QJuliaGenericField_qj, parity::QJuliaEnums.QJuliaParity_qj)
   if(field.field_desc.siteSubset == 1); return field.v; end

   parity_volume = complex_length*field.field_desc.real_volumeCB*(field.field_desc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY ? field.field_desc.nSpin : field.field_desc.nColor)*field.field_desc.nColor
   offset        = parity == QJuliaEnums.QJULIA_EVEN_PARITY ? 1 : parity_volume+1

   return view(field.v, offset:offset+parity_volume-1, :)
end

function field_info(field::QJuliaGenericField_qj; vrbs = true)

  field_desc = field.field_desc
  println(" ")
  println("General field info for ", typeof(field), ": ")

  if vrbs == true; QJuliaGrid.field_descr_info(field_desc); end

  println("Field data reference type : ", typeof(field.v))
  println("Field object pointer: ", pointer_from_objref(field))

  if field.v === missing
    println("Warning: data array was not allocated.")
  elseif field_desc.grid.location == QJuliaEnums.QJULIA_CUDA_FIELD_LOCATION
    println("Warning: data array was allocated on GPU")
    println("Physical length : ", length(field.v))
  else
    println("Row data pointer: ", pointer(field.v))
    println("Bytes : ", length(field.v))
  end

  println(" ")

end

# Specialized outer constructors:

# WARNING: transfered into CUDA impl
#CreateGenericField(fdesc::QJuliaGrid.QJuliaLatticeFieldDescr_qj) = QJuliaLatticeField_qj(fdesc)

CreateColorSpinor(grid::QJuliaGrid.QJuliaGridDescr_qj; parity::QJuliaEnums.QJuliaParity_qj = QJuliaEnums.QJULIA_INVALID_PARITY) = QJuliaLatticeField_qj(grid, QJuliaEnums.QJULIA_SCALAR_GEOMETRY, parity)

CreateBlockColorSpinor(grid::QJuliaGrid.QJuliaGridDescr_qj, NBlock; parity::QJuliaEnums.QJuliaParity_qj = QJuliaEnums.QJULIA_INVALID_PARITY) = QJuliaLatticeField_qj(grid, QJuliaEnums.QJULIA_SCALAR_GEOMETRY, parity; nBlock = NBlock)

CreateGaugeField(grid::QJuliaGrid.QJuliaGridDescr_qj) = QJuliaLatticeField_qj(grid, QJuliaEnums.QJULIA_VECTOR_GEOMETRY, QJuliaEnums.QJULIA_INVALID_PARITY)

#Create references to parity fields
Even(field::QJuliaGenericField_qj) = QJuliaLatticeField_qj(field, QJuliaEnums.QJULIA_EVEN_PARITY)
Odd(field::QJuliaGenericField_qj)  = QJuliaLatticeField_qj(field, QJuliaEnums.QJULIA_ODD_PARITY)

function (field::QJuliaLatticeField_qj)(i::Int)
  return field.v[i]
end

# Field Accessors
# input : field object
# out   : tuple of data array and args (which is a tuple of dof and offsets)

struct AccessorArgs
  #
  nParity::Int
  #
  nSpin::Int
  #
  nColor::Int
  #
  nVec::Int
  #
  length::Int
  #
  offset_cb::Int
  #
  field_order_id::Int
  # 1 for matter fields, 4 for gauge fields
  geom::Int
end

function create_field_accessor(field::QJuliaGenericField_qj; verbosity = true)

  nParity = field.field_desc.siteSubset
  nSpin   = field.field_desc.nSpin == nothing ? 0 : field.field_desc.nSpin
  nColor  = field.field_desc.nColor
  nVecs   = field.field_desc.nBlock

  field_order_id = 0 # currently default field order id : color inside spin inside space-time for spinors, color cols inside color rows inside space-time for gauge fields

  stride  = complex_length*field.field_desc.grid.grid_volume #include padding

  if     field.field_desc.geom == QJuliaEnums.QJULIA_SCALAR_GEOMETRY
    field_order_shape  = (stride, nColor, nSpin, nVecs)
    geom = 1
  elseif field.field_desc.geom == QJuliaEnums.QJULIA_VECTOR_GEOMETRY
    field_order_shape  = (stride, nColor, nColor, 4) #4 directions
    geom = 4
  elseif field.field_desc.geom == QJuliaEnums.QJULIA_COARSE_GEOMETRY
    field_order_shape  = (stride, nColor, nColor, nSpin, nSpin, 4) #4 directions
    geom = 8
  else
    error("Field geometry is not supported.")
  end

  if verbosity == true ; println("Create field accessor with dims ", field_order_shape);end

  field_order = reshape(field.v, field_order_shape)
  #
  offset_cb   = Int(length(field.v) / field.field_desc.siteSubset)
  #
  args        = AccessorArgs(nParity, nSpin, nColor, nVecs, stride, offset_cb, field_order_id, geom)
  #
  return (field_order, args)
end

end # module QJuliaFields

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# reach-set interfaces
include("AbstractReachSet.jl")
include("AbstractLazyReachSet.jl")
include("AbstractTaylorModelReachSet.jl")

# concrete implementations
include("ReachSet.jl")
include("SparseReachSet.jl")
include("ShiftedReachSet.jl")
include("TaylorModelReachSet.jl")
include("TemplateReachSet.jl")


text:
# kernels.jl
# ==========
# The @kernel macro and associated machinery
#

load("utils/req.jl")
req("utils/staged.jl")
req("utils/utils.jl")
req("dag/transforms.jl")
req("tangle.jl")
req("midsection.jl")
req("julia_backend.jl")

# -- Front half: code -> general DAG (no argument types available yet) --------

function ast_to_gendag(code)
    # Front end: ast --> dag
    @expect is_expr(code, :let)
    @expect length(code.args) == 1 # no let arguments, just body
    rawdag = tangle(code.args[1])

    # Front midsection: DAG transforms independent of argument types
    kernelargs = collect_arguments(rawdag)
    gendag = general_transform(rawdag)

    (kernelargs, gendag)
end

function collect_arguments(sink::Node)
    symnode_names = collect_symnode_names(sink)
    append(get(symnode_names, :output, Symbol{}),
           get(symnode_names, :input,  Symbol{}))
end

# Transform raw DAG independent of argument types
function general_transform(rawdag::Node)
    dag2 = scattered(rawdag)
end

# -- Back half: (general DAG, argtypes) -> specific kernel body ---------------

function gendag_to_specbody(gendag::Node, argnames::Vector{Symbol},
                            argtypes::Vector{Type})

    # todo: Add middle processing dependent on argtypes

    # Back end: Create julia code and wrap into for loops
    value, flat_code = untangle(gendag)

    indvars = (:_i, :_j, :_k, :_l) # todo: use gensyms instead

    nd = ndims(argtypes[1])
    body = wrap_kernel_body(flat_code, indvars[1:nd])
    body = :(shape = size($(argnames[1])); $body)
end

# -- @kernel ------------------------------------------------------------------

macro kernel(code)
    code_kernel(code)
end
function code_kernel(code)
    # Front end:
    #   process code as far as we can without argument types
    kernelargs, gendag = ast_to_gendag(code)

    #   create staged kernel function to invoke back half
    @gensym kernel
    @eval begin
        @staged function ($kernel)($kernelargs...)
            # Back half
            gendag_to_specbody(($gendag),
                               Symbol[{$quoted_exprs(kernelargs)...}...],
                               Type[{$kernelargs...}...])
        end
    end

    #    finally insert a call to it in the code
    quote
        ($kernel)($kernelargs...)
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
import Base: getindex
import FourierCore.Domain: Torus, Circle
import Base: size

export FourierGrid
export create_grid, getindex, size

struct FourierGrid{𝒢,𝒲,𝒟}
    nodes::𝒢
    wavenumbers::𝒲
    domain::𝒟
end

getindex(f::FourierGrid, i) = f.nodes[i]

function FourierGrid(grid_points, Ω::Circle; arraytype = Array) # perhaps change to match the product domain
    @assert length(grid_points) == 1
    grid = arraytype(nodes(grid_points, a = Ω.a, b = Ω.b))
    k⃗ = arraytype(wavenumbers(grid_points, L = Ω.b - Ω.a))
    return FourierGrid([grid], [k⃗], Ω)
end

"""
FourierGrid(grid_points, Ω::ProductDomain, arraytype=Array)
# Description
Create a numerical grid with grid_points resolution in the domain Ω \n
Only works for fully periodic grids at the moment
# Arguments
- `grid_points`: tuple | a tuple of ints in each direction for product domain
- `Ω`: ProductDomain   | a product domain object
# Keyword Arguments
- arraytype = Array
# Return
A Fourier Grid object
"""
function FourierGrid(grid_points::Tuple, Ω::Torus; arraytype = Array)
    @assert length(grid_points) == length(Ω.domains)
    grid = []
    k⃗ = []
    for (i, domain) in enumerate(Ω.domains)
        L = domain.b - domain.a
        reshape_dims = appropriate_dims(length(Ω.domains), i, grid_points[i])
        push!(grid, arraytype(reshape(nodes(grid_points[i], a = domain.a, b = domain.b), reshape_dims)))
        push!(k⃗, arraytype(reshape(wavenumbers(grid_points[i], L = L), reshape_dims)))
    end
    return FourierGrid(Tuple(grid), Tuple(k⃗), Ω)
end

FourierGrid(grid_points::Int, Ω::Torus; arraytype = Array) = FourierGrid(Tuple(fill(grid_points, ndims(Ω))), Ω, arraytype = arraytype)
FourierGrid(grid_points::Array, Ω::Torus; arraytype = Array) = FourierGrid(Tuple(grid_points), Ω, arraytype = arraytype)

function Base.show(io::IO, F::FourierGrid)
    println("domain:", F.domain)
    print("gridpoints:")
    if typeof(F.nodes[1]) <: AbstractArray
        for (i, grid) in enumerate(F.nodes)
            print(length(grid))
            if i != length(F.nodes)
                printstyled(io, "×")
            end
        end
    else
        print(length(F.nodes))
    end
end

size(f::FourierGrid) = length.(f.nodes)

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
using DiskArrays
using DiskArrays: ReshapedDiskArray, PermutedDiskArray
using Test

#Define a data structure that can be used for testing
struct _DiskArray{T,N,A<:AbstractArray{T,N}} <: AbstractDiskArray{T,N}
  getindex_count::Ref{Int}
  setindex_count::Ref{Int}
  parent::A
  chunksize::NTuple{N,Int}
end
_DiskArray(a;chunksize=size(a)) = _DiskArray(Ref(0),Ref(0),a,chunksize)
Base.size(a::_DiskArray) = size(a.parent)
DiskArrays.haschunks(::_DiskArray) = DiskArrays.Chunked()
DiskArrays.eachchunk(a::_DiskArray) = DiskArrays.GridChunks(a,a.chunksize)
getindex_count(a::_DiskArray) = a.getindex_count[]
setindex_count(a::_DiskArray) = a.setindex_count[]
trueparent(a::_DiskArray) = a.parent
getindex_count(a::ReshapedDiskArray) = getindex_count(a.parent)
setindex_count(a::ReshapedDiskArray) = setindex_count(a.parent)
trueparent(a::ReshapedDiskArray) = trueparent(a.parent)
getindex_count(a::PermutedDiskArray) = getindex_count(a.a.parent)
setindex_count(a::PermutedDiskArray) = setindex_count(a.a.parent)
trueparent(a::PermutedDiskArray{T,N,<:PermutedDimsArray{T,N,perm,iperm}}) where {T,N,perm,iperm} = permutedims(trueparent(a.a.parent),perm)
function DiskArrays.readblock!(a::_DiskArray,aout,i::AbstractUnitRange...)
  ndims(a) == length(i) || error("Number of indices is not correct")
  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
  #println("reading from indices ", join(string.(i)," "))
  a.getindex_count[] += 1
  aout .= a.parent[i...]
end
function DiskArrays.writeblock!(a::_DiskArray,v,i::AbstractUnitRange...)
  ndims(a) == length(i) || error("Number of indices is not correct")
  all(r->isa(r,AbstractUnitRange),i) || error("Not all indices are unit ranges")
  #println("Writing to indices ", join(string.(i)," "))
  a.setindex_count[] += 1
  view(a.parent,i...) .= v
end

function test_getindex(a)
  @test a[2,3,1] == 10
  @test a[2,3] == 10
  @test a[2,3,1,1] == 10
  @test a[:,1] == [1, 2, 3, 4]
  @test a[1:2, 1:2,1,1] == [1 5; 2 6]
  @test a[2:2:4,1:2:5] == [2 10 18; 4 12 20]
  @test a[end:-1:1,1,1] == [4,3,2,1]
  @test a[[1,3,4],[1,3],1] == [1 9; 3 11; 4 12]
  # Test bitmask indexing
  m = falses(4,5,1)
  m[2,:,1] .= true
  @test a[m] == [2,6,10,14,18]
  # Test linear indexing
  @test a[11:15] == 11:15
  @test a[20:-1:9] == 20:-1:9
  @test a[[3, 5, 8]] == [3, 5, 8]
  @test a[2:4:14] == [2,6,10,14]
  #Test that readblock was called exactly onces for every getindex
  @test getindex_count(a) == 13
end

function test_setindex(a)
  a[1,1,1] = 1
  a[1,2] = 2
  a[1,3,1,1] = 3
  a[2,:] = [1, 2, 3, 4, 5]
  a[3, 3:4,1,1] = [3,4]
  # Test bitmask indexing
  m = falses(4,5,1)
  m[4,:,1] .= true
  a[m] = [10,11,12,13,14]
  #Test that readblock was called exactly onces for every getindex
  @test setindex_count(a) == 6
  @test trueparent(a)[1,1:3,1] == [1,2,3]
  @test trueparent(a)[2,:,1] == [1,2,3,4,5]
  @test trueparent(a)[3,3:4,1] == [3,4]
  @test trueparent(a)[4,:,1] == [10,11,12,13,14]
  a[1:2:4,1:2:5,1] = [1 2 3; 5 6 7]
  @test trueparent(a)[1:2:4,1:2:5,1] == [1 2 3; 5 6 7]
  @test setindex_count(a) == 7
  a[[2,4],1:2,1] = [1 2; 5 6]
  @test trueparent(a)[[2,4],1:2,1] == [1 2; 5 6]
  @test setindex_count(a) == 8
end

function test_view(a)
  v = view(a,2:3,2:4,1)

  v[1:2,1] = [1,2]
  v[1:2,2:3] = [4 4; 4 4]
  @test v[1:2,1] == [1,2]
  @test v[1:2,2:3] == [4 4; 4 4]
  @test trueparent(a)[2:3,2] == [1,2]
  @test trueparent(a)[2:3,3:4] == [4 4; 4 4]
  @test getindex_count(a) == 2
  @test setindex_count(a) == 2
end

function test_reductions(af)
  data = rand(10,20,2)
  for f in (minimum,maximum,sum,
            (i,args...;kwargs...)->all(j->j>0.1,i,args...;kwargs...),
            (i,args...;kwargs...)->any(j->j<0.1,i,args...;kwargs...),
            (i,args...;kwargs...)->mapreduce(x->2*x,+,i,args...;kwargs...))
    a = af(data)
    @test isapprox(f(a),f(data))
    @test getindex_count(a) <= 10
    #And test reduction along dimensions
    a = _DiskArray(data,chunksize=(5,4,2))
    @test all(isapprox.(f(a,dims=2),f(data,dims=2)))
    #The minimum and maximum functions do some initialization, which will increase
    #the number of reads
    @test f in (minimum, maximum) || getindex_count(a) <= 12
    a = _DiskArray(data,chunksize=(5,4,2))
    @test all(isapprox.(f(a,dims=(1,3)),f(data,dims=(1,3))))
    @test f in (minimum, maximum) || getindex_count(a) <= 12
  end
end

function test_broadcast(a_disk1)
  a_disk2 = _DiskArray(rand(1:10,1,9), chunksize=(1,3))
  a_mem   = reshape(1:2,1,1,2);

  s = a_disk1 .+ a_disk2 .* Ref(2) ./ (2,)
  #Test lazy broadcasting
  @test s isa DiskArrays.BroadcastDiskArray
  @test s === DiskArrays.BroadcastDiskArray(s.bc)
  @test getindex_count(a_disk1)==0
  @test setindex_count(a_disk1)==0
  @test getindex_count(a_disk2)==0
  @test setindex_count(a_disk2)==0
  @test size(s)==(10,9,2)
  @test eltype(s) == Float64
  #Lets merge another broadcast
  s2 = s ./ a_mem
  @test s isa DiskArrays.BroadcastDiskArray
  @test getindex_count(a_disk1)==0
  @test getindex_count(a_disk2)==0
  @test size(s)==(10,9,2)
  @test eltype(s) == Float64
  #And now do the computation with Array as a sink
  aout = zeros(10,9,2)
  aout .= s2 .* 2 ./ Ref(2)
  #Test if the result is correct
  @test aout == (trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem
  @test getindex_count(a_disk1)==6
  @test getindex_count(a_disk2)==6
  #Now use another DiskArray as the output
  aout = _DiskArray(zeros(10,9,2),chunksize=(5,3,2))
  aout .= s ./ a_mem
  @test trueparent(aout) == (trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem
  @test setindex_count(aout)==6
  @test getindex_count(a_disk1)==12
  @test getindex_count(a_disk2)==12
  #Test reduction of broadcasted expression
  r = sum(s2, dims=(1,2))
  @test all(isapprox.(sum((trueparent(a_disk1) .+ trueparent(a_disk2))./a_mem,dims=(1,2)),r))
  @test getindex_count(a_disk1)==18
  @test getindex_count(a_disk2)==18
end

@testset "Index interpretation" begin
  import DiskArrays: DimsDropper, Reshaper
  a = zeros(3,3,1)
  @test interpret_indices_disk(a, (:,2,:)) == ((Base.OneTo(3), 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int}}((2,)))
  @test interpret_indices_disk(a, (1,2,:)) == ((1:1, 2:2, Base.OneTo(1)), DimsDropper{Tuple{Int,Int}}((1, 2)))
  @test interpret_indices_disk(a, (1,2,2,1)) == ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
  @test interpret_indices_disk(a, (1,2,2,1)) == ((1:1, 2:2, 2:2), DimsDropper{Tuple{Int,Int,Int}}((1, 2, 3)))
  @test interpret_indices_disk(a, (:,1:2)) == ((Base.OneTo(3), 1:2, 1:1), DimsDropper{Tuple{Int}}((3,)))
  @test interpret_indices_disk(a, (:,)) == ((Base.OneTo(3), Base.OneTo(3), Base.OneTo(1)), DiskArrays.Reshaper{Int}(9))
end

@testset "AbstractDiskArray getindex" begin
  a = _DiskArray(reshape(1:20,4,5,1))
  test_getindex(a)
end

@testset "AbstractDiskArray setindex" begin
  a = _DiskArray(zeros(Int,4,5,1))
  test_setindex(a)
end

@testset "Zerodimensional" begin
  a = _DiskArray(zeros(Int))
  @test a[] == 0
  @test a[1] == 0
  a[] = 5
  @test a[] == 5
  a[1] = 6
  @test a[] == 6
end

@testset "Views" begin
  a = _DiskArray(zeros(Int,4,5,1))
  test_view(a)
end

# The remaing tests only work for Julia >= 1.3
if VERSION >= v"1.3.0"
import Statistics: mean
@testset "Reductions" begin
  a = data -> _DiskArray(data,chunksize=(5,4,2))
  test_reductions(a)
end

@testset "Broadcast" begin
  a_disk1 = _DiskArray(rand(10,9,2), chunksize=(5,3,2))
  test_broadcast(a_disk1)
end

@testset "Broadcast with length 1 final dim" begin
  a_disk1 = _DiskArray(rand(10,9,1), chunksize=(5,3,1))
  a_disk2 = _DiskArray(rand(1:10,1,9), chunksize=(1,3))
  s = a_disk1 .+ a_disk2
  @test DiskArrays.eachchunk(s) isa DiskArrays.GridChunks{3}
  @test size(collect(s)) == (10, 9, 1)
end

@testset "Reshape" begin
  a = reshape(_DiskArray(reshape(1:20,4,5)),4,5,1)
  test_getindex(a)
  a = reshape(_DiskArray(zeros(Int,4,5)),4,5,1)
  test_setindex(a)
  a = reshape(_DiskArray(zeros(Int,4,5)),4,5,1)
  test_view(a)
  a = data -> reshape(_DiskArray(data,chunksize=(5,4,2)),10,20,2,1)
  test_reductions(a)
  a = reshape(_DiskArray(reshape(1:20,4,5)),4,5,1)
  @test ReshapedDiskArray(a.parent, a.keepdim, a.newsize) === a
end

import Base.PermutedDimsArrays.invperm
@testset "Permutedims" begin
  p = (3,1,2)
  ip = invperm(p)
  a = permutedims(_DiskArray(permutedims(reshape(1:20,4,5,1),ip)),p)
  test_getindex(a)
  a = permutedims(_DiskArray(zeros(Int,5,1,4)),p)
  test_setindex(a)
  a = permutedims(_DiskArray(zeros(Int,5,1,4)),p)
  test_view(a)
  a = data -> permutedims(_DiskArray(permutedims(data,ip),chunksize=(4,2,5)),p)
  test_reductions(a)
  a_disk1 = permutedims(_DiskArray(rand(9,2,10), chunksize=(3,2,5)),p)
  test_broadcast(a_disk1)
  @test PermutedDiskArray(a_disk1.a) === a_disk1
end

@testset "Unchunked String arrays" begin
  a = reshape(1:200000,200,1000)
  b = string.(a)
  c = collect(Union{Int,Missing},a)

  DiskArrays.default_chunk_size[] = 100
  DiskArrays.fallback_element_size[] = 100
  @test DiskArrays.estimate_chunksize(a) == (200,1000)
  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,1000))
  @test DiskArrays.estimate_chunksize(b) == (200,1000)
  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,1000))
  @test DiskArrays.estimate_chunksize(c) == (200,1000)
  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,1000))
  DiskArrays.default_chunk_size[] = 1
  @test DiskArrays.estimate_chunksize(a) == (200,625)
  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,625))
  @test DiskArrays.estimate_chunksize(b) == (200,50)
  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,50))
  @test DiskArrays.estimate_chunksize(c) == (200,625)
  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,625))
  DiskArrays.fallback_element_size[] = 1000
  @test DiskArrays.estimate_chunksize(a) == (200,625)
  @test DiskArrays.eachchunk(a) == DiskArrays.GridChunks(a,(200,625))
  @test DiskArrays.estimate_chunksize(b) == (200,5)
  @test DiskArrays.eachchunk(b) == DiskArrays.GridChunks(b,(200,5))
  @test DiskArrays.estimate_chunksize(c) == (200,625)
  @test DiskArrays.eachchunk(c) == DiskArrays.GridChunks(c,(200,625))
end

end


text:
using Knet, Random, Distributions

function generate(model, vocab, samples=1; prefix="")
    for n in 1:samples
        reset!(model)
        input = vocab.eos
        for token in split(prefix, "")
            model([input])
            input = vocab.w2i[token]
            print(token)
        end
        while true
            #output = model([input])
            #input = sample(output)
            probs = softmax(vec(model([input])))
            input = rand(Categorical(Array(probs)))
            if input == vocab.eos; break; end
            print(vocab.i2w[input])
        end
        println()
    end
end

function generate_old(model,chars,n)
    x = 1
    reset!(model)
    for i=1:n
        y = model([x])
        x = sample(y)
        print(chars[x])
    end
    println()
end

reset!(m)=(for r in m.layers; r isa RNN && (r.c=r.h=0); end);

function sample(y)
    p = Array(exp.(y)); r = rand()*sum(p)
    for j=1:length(p); (r -= p[j]) < 0 && return j; end
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
#
# objects.jl --
#
# Implement most part of OI-FITS objects behavior.
#
#------------------------------------------------------------------------------

"""
    OIFITS.float_type(obj)

yields the floating-point type used for OI-FITS object `obj`.

"""
float_type(::OIDataBlock{T}) where {T} = T
float_type(::OIMaster{T}) where {T} = T

"""
    OIFITS.contents(obj)

yields the vector of the data-blocks owned by OI-FITS master object `obj`.

"""
contents(obj::OIMaster) = getfield(obj, :all)

"""
    OIFITS.is_attached(db)

yields whether OI-FITS data-block `db` is attached to an OI-FITS master object.

    OIFITS.is_attached(db, owner)

yields whether OI-FITS data-block `db` is attached to the OI-FITS master object
`owner`.

"""
is_attached(db::OIDataBlock) = isdefined(db, :owner)
is_attached(db::OIDataBlock{T}, owner::OIMaster{T}) where {T} =
    (is_attached(db) && db.owner === owner)

# By default, create data in double precision.
OIMaster(args...) = OIMaster{Float64}(args...)
OITarget(; kwds...) = OITarget{Float64}(; kwds...)
OIArray(; kwds...) = OIArray{Float64}(; kwds...)
OIWavelength(; kwds...) = OIWavelength{Float64}(; kwds...)
OICorrelation(; kwds...) = OICorrelation{Float64}(; kwds...)
OIVis(; kwds...) = OIVis{Float64}(; kwds...)
OIVis2(; kwds...) = OIVis2{Float64}()
OIT3(; kwds...) = OIT3{Float64}(; kwds...)
OIFlux(; kwds...) = OIFlux{Float64}(; kwds...)
OIPolarization(; kwds...) = OIPolarization{Float64}(; kwds...)

OIMaster{T}(args::OIDataBlock...) where {T} = OIMaster{T}(args)
function OIMaster{T}(args::Union{AbstractVector{<:OIDataBlock},
                                 Tuple{Vararg{OIDataBlock}}}) where {T}
    push!(OIMaster{T}(), args)
end

"""
    OIFITS.convert_float_type(T, obj)

converts the floating-point type of array elements stored by OI-FITS object
`obj` to be `T`.  If floating-point type of `obj` is already `T`, `obj` itself
is returned; otherwise an unlinked object is returned.

""" convert_float_type

# Constructors for copy/conversion.
for T in (:OITarget, :OIArray, :OIWavelength, :OICorrelation,
          :OIVis, :OIVis2, :OIT3, :OIFlux, :OIPolarization)
    @eval begin
        # Constructor from an object of same kind yields same object if exactly
        # same type, a copy if floating-point type is different.
        $T{T}(obj::$T{T}) where {T} = obj
        $T{T}(obj::$T) where {T} = copyto!($T{T}(), obj)

        # Method to convert the floating-point type.
        convert_float_type(::Type{T}, obj::$T) where {T<:AbstractFloat} =
            $T{T}(obj)

        # Extend `convert` to call constructor.
        Base.convert(::Type{$T{T}}, obj::$T) where {T} = $T{T}(obj)

        # Copy to object of same kind.  Only defined and writable fields are
        # copied.
        function Base.copyto!(dst::$T{T}, src::$T) where {T}
            for sym in fieldnames($T{T})
                if iswritable($T{T}, sym) && isdefined(src, sym)
                    _setproperty!(dst, sym, getfield(src, sym))
                end
            end
            return dst
        end
    end
end

"""
    copy(db)

yields an unlinked copy of OI-FITS data-block `db`.  Only the fields of `db`
that are defined and are not links are initialized in the returned instance.
Defined data fields will be shared between `db` and the returned copy.  Links
such as `db.array`, `db.instr`, etc.  are left undefined in the returned object
which can be safely pushed into another `OIMaster` instance.  The `copyto!`
method has the same behavior of only copying writable and defined fields.

"""
Base.copy(obj::T) where {T<:OIDataBlock} = copyto!(T(), obj)

#------------------------------------------------------------------------------
# INDEXING

@inline Base.iswritable(::T, sym::Symbol) where {T<:OIDataBlock} =
    iswritable(T, sym)
@inline Base.iswritable(::Type{T}, sym::Symbol) where {T<:OIDataBlock} =
    isreadonly(T, sym) == false

@inline Base.isreadonly(::T, sym::Symbol) where {T<:OIDataBlock} =
    isreadonly(T, sym)
@inline Base.isreadonly(::Type{T}, sym::Symbol) where {T<:OIDataBlock} =
    (sym == :extname || sym == :owner)
@inline Base.isreadonly(::Type{T}, sym::Symbol) where {T<:OIData} =
    (sym == :extname || sym == :owner || sym == :array || sym == :instr ||
     sym == :correl)
@inline Base.isreadonly(::Type{T}, sym::Symbol) where {T<:OIPolarization} =
    (sym == :extname || sym == :owner || sym == :array)

# Extend `getproperty` and `setproperty!` to implement consistent `obj.field`
# syntax.  The behavior depends on the type of data-block.  Note that
# `setfield!` does no convert given value, it is the job of `setproperty!` to
# do that if needed.
@inline Base.getproperty(obj::OIDataBlock, sym::Symbol) =
    (sym == :extname ? get_extname(obj) : _getfield(obj, sym))

@inline Base.setproperty!(obj::T, sym::Symbol, val) where {T<:OIDataBlock} =
    (isreadonly(T, sym) ? read_only_field(T, sym) :
     _setproperty!(obj, sym, val))

# Like `getfield` but yields `nothing` if undefined.
@inline function _getfield(obj::T, sym::Symbol) where {T<:OIDataBlock}
    hasfield(T, sym) || type_has_no_field(T, sym)
    isdefined(obj, sym) ? getfield(obj, sym) : nothing
end

# Like `setproperty!` but by-pass checking whether field is writable.
@inline _setproperty!(obj::T, sym::Symbol, val) where {T} =
    setfield!(obj, sym, convert(fieldtype(T, sym), val))

@noinline type_has_no_field(T::Type, sym::Symbol) =
    error("type `", T, "` has no field `", sym, "`")

read_only_field(::T, sym::Symbol) where {T} =
    read_only_field(T, sym)

@noinline read_only_field(::Type{T}, sym::Symbol) where {T} =
    error("read-only field `", sym, "` in type `", T, "`")

# Make OIMaster usable as an interator and as a vector of data-blocks.
Base.length(obj::OIMaster) = length(contents(obj))
Base.size(obj::OIMaster) = (length(obj),)
Base.IndexStyle(::Type{<:OIMaster}) = IndexLinear()
Base.getindex(obj::OIMaster, i::Integer) = getindex(contents(obj), i)
Base.iterate(obj::OIMaster) = iterate(contents(obj))
Base.iterate(obj::OIMaster, state) = iterate(contents(obj), state)

# Iterating over a data-block yields (key,val) pairs where keys are the
# symbolic names of defined fields.

Base.iterate(obj::T) where {T<:OIDataBlock} =
    ((:extname, get_extname(obj)), (0, fieldnames(T)))

function Base.iterate(obj::OIDataBlock,
                      state::Tuple{Int,NTuple{N,Symbol}}) where {N}
    k, fields = state
    while k < N
        k += 1
        sym = fields[k]
        if isdefined(obj, sym)
            return ((sym, getfield(obj, sym)), (k, fields))
        end
    end
    nothing
end

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

# Known OI-FITS extension names.
const EXTNAMES = ("OI_TARGET", "OI_WAVELENGTH", "OI_ARRAY", "OI_VIS",
                  "OI_VIS2", "OI_T3", "OI_FLUX", "OI_CORR", "OI_INSPOL")

"""
    OIFITS.get_datablock_type([T = Float64,] extname)

yields the data-block type `<:OIDataBlock{T}` associated to the OI-FITS
extension `extname`.

"""
get_datablock_type(extname::Union{AbstractString,Symbol}) =
    get_datablock_type(Float64, extname)

get_datablock_type(::Type{T}, extname::Symbol) where {T<:AbstractFloat} =
    (extname === :OI_TARGET     ? OITarget{T}       :
     extname === :OI_WAVELENGTH ? OIWavelength{T}   :
     extname === :OI_ARRAY      ? OIArray{T}        :
     extname === :OI_VIS        ? OIVis{T}          :
     extname === :OI_VIS2       ? OIVis2{T}         :
     extname === :OI_T3         ? OIT3{T}           :
     extname === :OI_FLUX       ? OIFlux{T}         :
     extname === :OI_CORR       ? OICorrelation{T}  :
     extname === :OI_INSPOL     ? OIPolarization{T} :
     bad_extname(extname))

get_datablock_type(::Type{T}, extname::AbstractString) where {T<:AbstractFloat} =
    (extname == "OI_TARGET"     ? OITarget{T}       :
     extname == "OI_WAVELENGTH" ? OIWavelength{T}   :
     extname == "OI_ARRAY"      ? OIArray{T}        :
     extname == "OI_VIS"        ? OIVis{T}          :
     extname == "OI_VIS2"       ? OIVis2{T}         :
     extname == "OI_T3"         ? OIT3{T}           :
     extname == "OI_FLUX"       ? OIFlux{T}         :
     extname == "OI_CORR"       ? OICorrelation{T}  :
     extname == "OI_INSPOL"     ? OIPolarization{T} :
     bad_extname(extname))

@noinline bad_extname(extname::Union{Symbol,AbstractString}) =
    error("bad OI-FITS extension name \"", extname, "\"")

"""
    OIFITS.get_extname(arg) -> str

yields the FITS extension name of an OI-FITS data-block corresponding to `arg`.
Argument can be an OI-FITS data-block instance, an OI-FITS data-block type, or
a FITS Header Data Unit.  If `arg` is an OI-FITS data-block instance,
`arg.extname` yields the same result.

"""
get_extname(::T) where {T<:OIDataBlock} = get_extname(T)
get_extname(::Type{<:OITarget})       = "OI_TARGET"
get_extname(::Type{<:OIWavelength})   = "OI_WAVELENGTH"
get_extname(::Type{<:OIArray})        = "OI_ARRAY"
get_extname(::Type{<:OIVis})          = "OI_VIS"
get_extname(::Type{<:OIVis2})         = "OI_VIS2"
get_extname(::Type{<:OIT3})           = "OI_T3"
get_extname(::Type{<:OIFlux})         = "OI_FLUX"
get_extname(::Type{<:OICorrelation})  = "OI_CORR"
get_extname(::Type{<:OIPolarization}) = "OI_INSPOL"
get_extname(hdr::FITSHeader) =
    (get_hdu_type(hdr) !== :binary_table ? "" :
     fix_name(get_string(hdr, "EXTNAME", "")))

function show(io::IO, db::OITarget)
    print(io, db.extname, ": ")
    if (tgt = db.target) !== nothing
        n = length(tgt)
        if n > 0
            print(io, "target=")
            for i in 1:n
                print(io, (i == 1 ? "[\"" : ", \""),
                      tgt[i], (i == n ? "\"]" : "\""))
            end
        end
    else
        n = 0
    end
    n == 0 && print(io, "<empty>")
end

function show(io::IO, db::OIArray)
    print(io, db.extname, ": ")
    ntels =  (db.sta_index === nothing ? 0 : length(db.sta_index))
    if db.arrname !== nothing
        print(io, "arrname=\"", db.arrname, "\", ", ntels, " telescope(s)")
    else
        print(io, ntels, " telescope(s)")
    end
end

function show(io::IO, db::OIWavelength)
    print(io, db.extname, ": ")
    if db.insname === nothing
        print(io, "<empty>")
    else
        print(io, "insname=\"", db.insname, "\"")
        if db.eff_wave !== nothing
            wave = db.eff_wave
            n = length(db.eff_wave)
            if n < 1
                print(io, " with 0 spectral channels")
            elseif n == 1
                print(io, " with 1 spectral channel at ",
                      round(wave[1]*1e6, digits=3), " µm")
            else
                print(io, " with ", n, " spectral channels from ",
                      round(minimum(wave)*1e6, digits=3), " µm to ",
                      round(maximum(wave)*1e6, digits=3), " µm")
            end
        end
    end
end

function show(io::IO, db::OIFlux)
    print(io, db.extname)
end

function _show(io::IO, db::OIDataBlock, nwaves::Integer, ntimes::Integer)
    print(io, db.extname, ": ")
    if nwaves > 0 && ntimes > 0
        print(io, nwaves*ntimes, " measurements in ", nwaves,
              " spectral channel(s) and ", ntimes, " exposure(s)")
    else
        print(io, "<empty>")
    end
end

function show(io::IO, db::OIVis)
    if db.visamp !== nothing
        dims = size(db.visamp)
        nwaves = dims[1]
        ntimes = dims[2]
    elseif db.visphi !== nothing
        dims = size(db.visphi)
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    _show(io, db, nwaves, ntimes)
end

function show(io::IO, db::OIVis2)
    if db.vis2data !== nothing
        dims = size(db.vis2data)
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    _show(io, db, nwaves, ntimes)
end

function show(io::IO, db::OIT3)
    if db.t3amp !== nothing
        dims = size(db.t3amp)
        nwaves = dims[1]
        ntimes = dims[2]
    elseif db.t3phi !== nothing
        dims = size(db.visphi)
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    _show(io, db, nwaves, ntimes)
end

function show(io::IO, master::OIMaster)
    print(io, "OI_MASTER: (", length(master), " data-block(s))")
    for db in master
        print(io, "\n    ")
        show(io, db)
    end
end

##############################
# METHODS FOR DATA SELECTION #
##############################

"""
     OIFITS.select(obj, args)

yields an iterable over the data-blocks of an OI-FITS master object `obj` whose
extension names are specifed in `args`.

"""
function select(obj::OIMaster{T}, args::AbstractString...) where {T}
    datablocks = Array{OIDataBlock{T}}(undef, 0)
    for db in obj
        if db.extname ∈ args
            push!(datablocks, db)
        end
    end
    return datablocks
end

"""
    OIFITS.select_target(src, id)

selects a subset of data from source `src` corresponding to the target
identified by `id` (an integer or a name).  The source can be an instance of
`OIMaster` or of any `OIDataBlock` sub-types.  The result is of the same type
as the source or `nothing` if the source contains no data for the given target.

The result may share part of its contents with the source `src`.

"""
select_target(obj::Union{OIMaster,OIDataBlock}, id::AbstractString) =
    _select_target(obj, fix_name(id))
select_target(obj::Union{OIMaster,OIDataBlock}, id::Integer) =
    _select_target(obj, to_integer(id))

function _select_target(obj::Union{OIMaster,OIDataBlock}, str::String)
    tgt = obj.target
    tgt === nothing && return nothing
    idx = findfirst(x -> same_name(x, str), tgt.target)
    idx === nothing && return nothing
    select_target(obj, tgt.target_id[idx])
end

function _select_target(src::OIMaster, id::Int)
    dst = OIMaster()
    for db in src
        Builder._push!(dst, select_target(db, id))
    end
    Builder._update_links!(dst)
end

_select_target(src::OIDataBlock, id::Int) = copy(src)

function _select_target(src::T, id::Int) where {T<:OITarget}
    idx = findfirst(x -> x == id, src.target_id)
    idx == nothing && return nothing
    (name, revn, defn) = Builder.get_description(src)
    dst = T()
    for (key, val) in src
        let spec = get(defn.spec, key, nothing)
            spec === nothing && continue
            if spec.iskeyword
                setfield!(dst, key, val)
            else
                setfield!(dst, key, val[idx:idx])
            end
        end
    end
    return dst
end

function _select_target(src::T, id::Int) where {T<:Union{OIVis,OIVis2,OIT3,
                                                        OIFlux,
                                                        OIPolarization}}
    # Pre-select target.
    target_id = src.target_id
    sel = findall(x -> x == id, target_id)
    length(sel) > 0 || return nothing

    # Build a data-block with selected data starting with an empty structure.
    dst = T()
    if length(sel) == length(target_id)
        # Just copy everythings.
        copyto!(dst, src)
    else
        (name, revn, defn) = Builder.get_description(src)
        for (key, val) in src
            let spec = get(defn.spec, key, nothing)
                spec != nothing || continue
                if spec.iskeyword
                    setfield!(dst, key, val)
                else
                    # Copy a sub-array corresponding to the selection.  (The target
                    # is specified for each last index.)
                    @assert last(size(val)) == length(target_id)
                    setfield!(dst, key, _select_last(arr, sel))
                end
            end
        end
    end
    return dst
end

function _select_last(arr::AbstractArray{T,N}, sel::Vector{Int}) where {T,N}
    siz = ntuple(i -> i == N ? length(sel) : size(arr, i), Val(N))
    sub = Array{T}(undef, siz)
    if N == 1
        for j in 1:length(sel)
            sub[j] = arr[sel[j]]
        end
    elseif N == 2
        for j in 1:length(sel)
            sub[:,j] = arr[:,sel[j]]
        end
    else
        error("unexpected rank ", N)
    end
    return sub
end

function _select_first(arr::AbstractArray{T,N}, sel::Vector{Int}) where {T,N}
    siz = ntuple(i -> i == 1 ? length(sel) : size(val, i), Val(N))
    sub = Array{T}(undef, siz)
    if N == 1
        for j in 1:length(sel)
            sub[j] = val[sel[j]]
        end
    elseif N == 2
        for j in 1:length(sel)
            sub[j,:] = val[sel[j],:]
        end
    else
        error("unexpected rank ", N)
    end
    return sub
end

"""
    OIFITS.select_wavelength(src, wmin, wmax)

selects a subset of data from source `src` corresponding to the measurements at
wavelengths `w` such that `wmin ≤ w ≤ wmax`.  The wavelength bounds are in the
same units as assumed by OI-FITS convention, that is in meters.  The source can
be an instance of `OIMaster` or of any `OIDataBlock` sub-types.

An alternative is to use a selector function:

    OIFITS.select_wavelength(src, sel)

where `sel` is callable and such that `sel(w)` returns whether wavelength `w`
is to be selected.

In any cases, the result is of the same type as the source or `nothing` if the
source contains no data in the given wavelength range.

The result may share part of its contents with the source `src`.

"""
function select_wavelength(src::Union{OIMaster,OIDataBlock},
                           wavemin::Real, wavemax::Real)
    wmin = convert(Float64, wavemin)
    wmax = convert(Float64, wavemax)
    select_wavelength(src, w -> wmin ≤ w ≤ wmax)
end

function select_wavelength(src::OIMaster, select::Function)
    dst = OIMaster()
    for db in src
        Builder._push!(dst, select_wavelength(db, select))
    end
    Builder._update_links!(dst)
end

select_wavelength(src::OIDataBlock, select::Function) = copy(src)

function select_wavelength(src::T,
                           select::Function) where {T<:Union{OIWavelength,
                                                             OIVis,OIVis2,OIT3,
                                                             OIFlux,
                                                             OIPolarization}}
    # Pre-select wavelengths.
    isintr = isa(src, OIWavelength) # source is instrument?
    wave = (isintr ? src.eff_wave : src.instr.eff_wave)
    wave === nothing && return nothing
    sel = findall(select, wave)
    length(sel) > 0 || return nothing

    # Build a data-block with selected data starting with an empty structure.
    dst = T()
    if length(sel) == length(wave)
        # Just copy everythings.
        copyto!(dst, src)
    else
        (name, revn, defn) = Builder.get_description(src)
        for (key, val) in src
            # Get field specifications.  Do nothing if field.
            let spec = get(defn.spec, key, nothing)
                spec === nothing && continue
                if spec.iskeyword || (spec.multiplier ≥ 0 && isintr)
                    # Field is a keyword or is not a wavelength-wise array.
                    setfield!(dst, key, val)
                else
                    # Copy a sub-array corresponding to the selection.  (The
                    # wavelength corresponds to the first index.)
                    @assert first(size(val)) == length(wave)
                    setfield!(dst, key, _select_first(val, sel))
                end
            end
        end
    end
    return dst
end


max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
function Base.findfirst(A::Vector{T}, name::AbstractString) where T <: Union{Dict, NamedTuple}
    seen = Set{String }()
    num = length(A)
    res = Vector{Int}()
    @inbounds for x in 1:num
        if A[x][name] ∉ seen
            push!(seen, A[x][name])
            push!(res, x)
        end
    end
    res
end

Base.unique(A::Vector{T}, name::AbstractString) where T <: Union{Dict, NamedTuple} = findfirst(A, name) |> q -> view(A, q)

Mytype = Union{Symbol , Expr}
extra(d::Expr) = d.head == :(::) ? d.args[1] : d
extra(d::Symbol) = d
extra(d::AbstractString) = d

function rep_func(a::Mytype , b::Mytype)

    tpe  = a.head == :(::) ? a.args[2] : :AbstractString
    ab ,aa = extra(b), extra(a)

    func = Expr(:call , :rep , a, b )
    b1   = Expr(:(=), :res , Expr(:call , Expr(:curly, :Vector , tpe) , :undef, ab) )
    b2   = Expr(:macrocall, Symbol("@inbounds"), "",
                    Expr(:for , Expr(:(=), :i , Expr(:call, :(:), 1, ab ) ),
                       Expr(:block, Expr(:(=) , Expr(:ref, :res, :i) , aa ) ) ))

    Expr(:function , func, Expr(:block,  b1, b2, :res) )

end

macro rep(a, b )
   esc(rep_func(a,b))
end

@rep(data::Dict, num::Integer)
@rep(data::Number, num::Integer)
@rep(data::Tuple, num::Integer)
@rep(data::AbstractString, num::Integer)
@rep(data::AbstractArray, num::Integer)
@rep(data::AbstractChar, num::Integer)

macro set(t...)
           collect_to_set(collect(t))
end

@inline collect_to_set(t::Vector) = Expr(:call, :Set , Expr(:vect, t...))

max_stars_repo_name: frankill/ftest

max_issues_repo_name: frankill/ftest

max_forks_repo_name: frankill/ftest