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

## 1-dimensional ranges ##

abstract Range{T} <: AbstractArray{T,1}

## ordinal ranges

abstract OrdinalRange{T,S} <: Range{T}
abstract AbstractUnitRange{T} <: OrdinalRange{T,Int}

immutable StepRange{T,S} <: OrdinalRange{T,S}
    start::T
    step::S
    stop::T

    function StepRange(start::T, step::S, stop::T)
        new(start, step, steprange_last(start,step,stop))
    end
end

# to make StepRange constructor inlineable, so optimizer can see `step` value
function steprange_last{T}(start::T, step, stop)
    if isa(start,AbstractFloat) || isa(step,AbstractFloat)
        throw(ArgumentError("StepRange should not be used with floating point"))
    end
    z = zero(step)
    step == z && throw(ArgumentError("step cannot be zero"))

    if stop == start
        last = stop
    else
        if (step > z) != (stop > start)
            last = steprange_last_empty(start, step, stop)
        else
            diff = stop - start
            if T<:Signed && (diff > zero(diff)) != (stop > start)
                # handle overflowed subtraction with unsigned rem
                if diff > zero(diff)
                    remain = -convert(T, unsigned(-diff) % step)
                else
                    remain = convert(T, unsigned(diff) % step)
                end
            else
                remain = steprem(start,stop,step)
            end
            last = stop - remain
        end
    end
    last
end

function steprange_last_empty{T<:Integer}(start::T, step, stop)
    # empty range has a special representation where stop = start-1
    # this is needed to avoid the wrap-around that can happen computing
    # start - step, which leads to a range that looks very large instead
    # of empty.
    if step > zero(step)
        last = start - one(stop-start)
    else
        last = start + one(stop-start)
    end
    last
end
# For types where x+one(x) may not be well-defined
steprange_last_empty(start, step, stop) = start - step

steprem(start,stop,step) = (stop-start) % step

StepRange{T,S}(start::T, step::S, stop::T) = StepRange{T,S}(start, step, stop)

immutable UnitRange{T<:Real} <: AbstractUnitRange{T}
    start::T
    stop::T
    UnitRange(start, stop) = new(start, unitrange_last(start,stop))
end
UnitRange{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop)

unitrange_last(::Bool, stop::Bool) = stop
unitrange_last{T<:Integer}(start::T, stop::T) =
    ifelse(stop >= start, stop, convert(T,start-one(stop-start)))
unitrange_last{T}(start::T, stop::T) =
    ifelse(stop >= start, convert(T,start+floor(stop-start)),
                          convert(T,start-one(stop-start)))

"""
    Base.OneTo(n)

Define an `AbstractUnitRange` that behaves like `1:n`, with the added
distinction that the lower limit is guaranteed (by the type system) to
be 1.
"""
immutable OneTo{T<:Integer} <: AbstractUnitRange{T}
    stop::T
    OneTo(stop) = new(max(zero(T), stop))
end
OneTo{T<:Integer}(stop::T) = OneTo{T}(stop)

colon(a::Real, b::Real) = colon(promote(a,b)...)

colon{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop)

range(a::Real, len::Integer) = UnitRange{typeof(a)}(a, oftype(a, a+len-1))

colon{T}(start::T, stop::T) = StepRange(start, one(stop-start), stop)

range{T}(a::T, len::Integer) =
    StepRange{T, typeof(a-a)}(a, one(a-a), a+oftype(a-a,(len-1)))

# first promote start and stop, leaving step alone
# this is for non-numeric ranges where step can be quite different
colon{A<:Real,C<:Real}(a::A, b, c::C) = colon(convert(promote_type(A,C),a), b, convert(promote_type(A,C),c))

"""
    colon(start, [step], stop)

Called by `:` syntax for constructing ranges.
"""
colon{T<:Real}(start::T, step, stop::T) = StepRange(start, step, stop)

"""
    :(start, [step], stop)

Range operator. `a:b` constructs a range from `a` to `b` with a step size of 1, and `a:s:b`
is similar but uses a step size of `s`. These syntaxes call the function `colon`. The colon
is also used in indexing to select whole dimensions.
"""
colon{T<:Real}(start::T, step::T, stop::T) = StepRange(start, step, stop)
colon{T<:Real}(start::T, step::Real, stop::T) = StepRange(promote(start, step, stop)...)

colon{T}(start::T, step, stop::T) = StepRange(start, step, stop)

"""
    range(start, [step], length)

Construct a range by length, given a starting value and optional step (defaults to 1).
"""
range{T,S}(a::T, step::S, len::Integer) = StepRange{T,S}(a, step, convert(T, a+step*(len-1)))

## floating point ranges

immutable FloatRange{T<:AbstractFloat} <: Range{T}
    start::T
    step::T
    len::T
    divisor::T
end
FloatRange(a::AbstractFloat, s::AbstractFloat, l::Real, d::AbstractFloat) =
    FloatRange{promote_type(typeof(a),typeof(s),typeof(d))}(a,s,l,d)

# float rationalization helper
function rat(x)
    y = x
    a = d = 1
    b = c = 0
    m = maxintfloat(Float32)
    while abs(y) <= m
        f = trunc(Int,y)
        y -= f
        a, c = f*a + c, a
        b, d = f*b + d, b
        max(abs(a),abs(b)) <= convert(Int,m) || return c, d
        oftype(x,a)/oftype(x,b) == x && break
        y = inv(y)
    end
    return a, b
end

function colon{T<:AbstractFloat}(start::T, step::T, stop::T)
    step == 0 && throw(ArgumentError("range step cannot be zero"))
    start == stop && return FloatRange{T}(start,step,1,1)
    (0 < step) != (start < stop) && return FloatRange{T}(start,step,0,1)

    # float range "lifting"
    r = (stop-start)/step
    n = round(r)
    lo = prevfloat((prevfloat(stop)-nextfloat(start))/n)
    hi = nextfloat((nextfloat(stop)-prevfloat(start))/n)
    if lo <= step <= hi
        a0, b = rat(start)
        a = convert(T,a0)
        if a/convert(T,b) == start
            c0, d = rat(step)
            c = convert(T,c0)
            if c/convert(T,d) == step
                e = lcm(b,d)
                a *= div(e,b)
                c *= div(e,d)
                eT = convert(T,e)
                if (a+n*c)/eT == stop
                    return FloatRange{T}(a, c, n+1, eT)
                end
            end
        end
    end
    FloatRange{T}(start, step, floor(r)+1, one(step))
end

colon{T<:AbstractFloat}(a::T, b::T) = colon(a, one(a), b)

colon{T<:Real}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...)
colon{T<:AbstractFloat}(a::T, b::AbstractFloat, c::T) = colon(promote(a,b,c)...)
colon{T<:AbstractFloat}(a::T, b::Real, c::T) = colon(promote(a,b,c)...)

range(a::AbstractFloat, len::Integer) = FloatRange(a,one(a),len,one(a))
range(a::AbstractFloat, st::AbstractFloat, len::Integer) = FloatRange(a,st,len,one(a))
range(a::Real, st::AbstractFloat, len::Integer) = FloatRange(float(a), st, len, one(st))
range(a::AbstractFloat, st::Real, len::Integer) = FloatRange(a, float(st), len, one(a))

## linspace and logspace

immutable LinSpace{T<:AbstractFloat} <: Range{T}
    start::T
    stop::T
    len::T
    divisor::T
end

function linspace{T<:AbstractFloat}(start::T, stop::T, len::T)
    len == round(len) || throw(InexactError())
    0 <= len || error("linspace($start, $stop, $len): negative length")
    if len == 0
        n = convert(T, 2)
        if isinf(n*start) || isinf(n*stop)
            start /= n; stop /= n; n = one(T)
        end
        return LinSpace(-start, -stop, -one(T), n)
    end
    if len == 1
        start == stop || error("linspace($start, $stop, $len): endpoints differ")
        return LinSpace(-start, -start, zero(T), one(T))
    end
    n = convert(T, len - 1)
    len - n == 1 || error("linspace($start, $stop, $len): too long for $T")
    a0, b = rat(start)
    a = convert(T,a0)
    if a/convert(T,b) == start
        c0, d = rat(stop)
        c = convert(T,c0)
        if c/convert(T,d) == stop
            e = lcm(b,d)
            a *= div(e,b)
            c *= div(e,d)
            s = convert(T,n*e)
            if isinf(a*n) || isinf(c*n)
                s, p = frexp(s)
                p2 = oftype(s,2)^p
                a /= p2; c /= p2
            end
            if a*n/s == start && c*n/s == stop
                return LinSpace(a, c, len, s)
            end
        end
    end
    a, c, s = start, stop, n
    if isinf(a*n) || isinf(c*n)
        s, p = frexp(s)
        p2 = oftype(s,2)^p
        a /= p2; c /= p2
    end
    if a*n/s == start && c*n/s == stop
        return LinSpace(a, c, len, s)
    end
    return LinSpace(start, stop, len, n)
end
function linspace{T<:AbstractFloat}(start::T, stop::T, len::Real)
    T_len = convert(T, len)
    T_len == len || throw(InexactError())
    linspace(start, stop, T_len)
end

"""
    linspace(start::Real, stop::Real, n::Real=50)

Construct a range of `n` linearly spaced elements from `start` to `stop`.

```jldoctest
julia> linspace(1.3,2.9,9)
9-element LinSpace{Float64}:
 1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9
```
"""
linspace(start::Real, stop::Real, len::Real=50) =
    linspace(promote(AbstractFloat(start), AbstractFloat(stop))..., len)

function show(io::IO, r::LinSpace)
    print(io, "linspace(")
    show(io, first(r))
    print(io, ',')
    show(io, last(r))
    print(io, ',')
    show(io, length(r))
    print(io, ')')
end

"""
`print_range(io, r)` prints out a nice looking range r in terms of its elements
as if it were `collect(r)`, dependent on the size of the
terminal, and taking into account whether compact numbers should be shown.
It figures out the width in characters of each element, and if they
end up too wide, it shows the first and last elements separated by a
horizontal elipsis. Typical output will look like `1.0,2.0,3.0,…,4.0,5.0,6.0`.

`print_range(io, r, pre, sep, post, hdots)` uses optional
parameters `pre` and `post` characters for each printed row,
`sep` separator string between printed elements,
`hdots` string for the horizontal ellipsis.
"""
function print_range(io::IO, r::Range,
                     pre::AbstractString = " ",
                     sep::AbstractString = ",",
                     post::AbstractString = "",
                     hdots::AbstractString = ",\u2026,") # horiz ellipsis
    # This function borrows from print_matrix() in show.jl
    # and should be called by show and display
    limit = get(io, :limit, false)
    sz = displaysize(io)
    if !haskey(io, :compact)
        io = IOContext(io, compact=true)
    end
    screenheight, screenwidth = sz[1] - 4, sz[2]
    screenwidth -= length(pre) + length(post)
    postsp = ""
    sepsize = length(sep)
    m = 1 # treat the range as a one-row matrix
    n = length(r)
    # Figure out spacing alignments for r, but only need to examine the
    # left and right edge columns, as many as could conceivably fit on the
    # screen, with the middle columns summarized by horz, vert, or diag ellipsis
    maxpossiblecols = div(screenwidth, 1+sepsize) # assume each element is at least 1 char + 1 separator
    colsr = n <= maxpossiblecols ? (1:n) : [1:div(maxpossiblecols,2)+1; (n-div(maxpossiblecols,2)):n]
    rowmatrix = r[colsr]' # treat the range as a one-row matrix for print_matrix_row
    A = alignment(io, rowmatrix, 1:m, 1:length(rowmatrix), screenwidth, screenwidth, sepsize) # how much space range takes
    if n <= length(A) # cols fit screen, so print out all elements
        print(io, pre) # put in pre chars
        print_matrix_row(io,rowmatrix,A,1,1:n,sep) # the entire range
        print(io, post) # add the post characters
    else # cols don't fit so put horiz ellipsis in the middle
        # how many chars left after dividing width of screen in half
        # and accounting for the horiz ellipsis
        c = div(screenwidth-length(hdots)+1,2)+1 # chars remaining for each side of rowmatrix
        alignR = reverse(alignment(io, rowmatrix, 1:m, length(rowmatrix):-1:1, c, c, sepsize)) # which cols of rowmatrix to put on the right
        c = screenwidth - sum(map(sum,alignR)) - (length(alignR)-1)*sepsize - length(hdots)
        alignL = alignment(io, rowmatrix, 1:m, 1:length(rowmatrix), c, c, sepsize) # which cols of rowmatrix to put on the left
        print(io, pre)   # put in pre chars
        print_matrix_row(io, rowmatrix,alignL,1,1:length(alignL),sep) # left part of range
        print(io, hdots) # horizontal ellipsis
        print_matrix_row(io, rowmatrix,alignR,1,length(rowmatrix)-length(alignR)+1:length(rowmatrix),sep) # right part of range
        print(io, post)  # post chars
    end
end

"""
    logspace(start::Real, stop::Real, n::Integer=50)

Construct a vector of `n` logarithmically spaced numbers from `10^start` to `10^stop`.

```jldoctest
julia> logspace(1.,10.,5)
5-element Array{Float64,1}:
   10.0
 1778.28
    3.16228e5
    5.62341e7
    1.0e10
```
"""
logspace(start::Real, stop::Real, n::Integer=50) = 10.^linspace(start, stop, n)

## interface implementations

size(r::Range) = (length(r),)

isempty(r::StepRange) =
    (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start))
isempty(r::AbstractUnitRange) = first(r) > last(r)
isempty(r::FloatRange) = length(r) == 0
isempty(r::LinSpace) = length(r) == 0

"""
    step(r)

Get the step size of a [`Range`](:obj:`Range`) object.
```jldoctest
julia> step(1:10)
1

julia> step(1:2:10)
2

julia> step(2.5:0.3:10.9)
0.3

julia> step(linspace(2.5,10.9,85))
0.1
```
"""
step(r::StepRange) = r.step
step(r::AbstractUnitRange) = 1
step(r::FloatRange) = r.step/r.divisor
step{T}(r::LinSpace{T}) = ifelse(r.len <= 0, convert(T,NaN), (r.stop-r.start)/r.divisor)

unsafe_length(r::Range) = length(r)  # generic fallback

function unsafe_length(r::StepRange)
    n = Integer(div(r.stop+r.step - r.start, r.step))
    isempty(r) ? zero(n) : n
end
length(r::StepRange) = unsafe_length(r)
unsafe_length(r::AbstractUnitRange) = Integer(last(r) - first(r) + 1)
unsafe_length(r::OneTo) = r.stop
length(r::AbstractUnitRange) = unsafe_length(r)
length(r::OneTo) = unsafe_length(r)
length(r::FloatRange) = Integer(r.len)
length(r::LinSpace) = Integer(r.len + signbit(r.len - 1))

function length{T<:Union{Int,UInt,Int64,UInt64}}(r::StepRange{T})
    isempty(r) && return zero(T)
    if r.step > 1
        return checked_add(convert(T, div(unsigned(r.stop - r.start), r.step)), one(T))
    elseif r.step < -1
        return checked_add(convert(T, div(unsigned(r.start - r.stop), -r.step)), one(T))
    else
        checked_add(div(checked_sub(r.stop, r.start), r.step), one(T))
    end
end

function length{T<:Union{Int,Int64}}(r::AbstractUnitRange{T})
    @_inline_meta
    checked_add(checked_sub(last(r), first(r)), one(T))
end
length{T<:Union{Int,Int64}}(r::OneTo{T}) = T(r.stop)

length{T<:Union{UInt,UInt64}}(r::AbstractUnitRange{T}) =
    r.stop < r.start ? zero(T) : checked_add(last(r) - first(r), one(T))

# some special cases to favor default Int type
let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    global length

    function length{T <: smallint}(r::StepRange{T})
        isempty(r) && return Int(0)
        div(Int(r.stop)+Int(r.step) - Int(r.start), Int(r.step))
    end

    length{T <: smallint}(r::AbstractUnitRange{T}) = Int(last(r)) - Int(first(r)) + 1
    length{T <: smallint}(r::OneTo{T}) = Int(r.stop)
end

first{T}(r::OrdinalRange{T}) = convert(T, r.start)
first{T}(r::OneTo{T}) = one(T)
first{T}(r::FloatRange{T}) = convert(T, r.start/r.divisor)
first{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.start/r.divisor)

last{T}(r::OrdinalRange{T}) = convert(T, r.stop)
last{T}(r::FloatRange{T}) = convert(T, (r.start + (r.len-1)*r.step)/r.divisor)
last{T}(r::LinSpace{T}) = convert(T, (r.len-1)*r.stop/r.divisor)

minimum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r)
maximum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r)
minimum(r::Range)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : min(first(r), last(r))
maximum(r::Range)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : max(first(r), last(r))

ctranspose(r::Range) = [x for _=1:1, x=r]
transpose(r::Range) = r'

# Ranges are immutable
copy(r::Range) = r

## iteration

start(r::FloatRange) = 0
done(r::FloatRange, i::Int) = length(r) <= i
next{T}(r::FloatRange{T}, i::Int) =
    (convert(T, (r.start + i*r.step)/r.divisor), i+1)

start(r::LinSpace) = 1
done(r::LinSpace, i::Int) = length(r) < i
next{T}(r::LinSpace{T}, i::Int) =
    (convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor), i+1)

start(r::StepRange) = oftype(r.start + r.step, r.start)
next{T}(r::StepRange{T}, i) = (convert(T,i), i+r.step)
done{T,S}(r::StepRange{T,S}, i) = isempty(r) | (i < min(r.start, r.stop)) | (i > max(r.start, r.stop))
done{T,S}(r::StepRange{T,S}, i::Integer) =
    isempty(r) | (i == oftype(i, r.stop) + r.step)

start{T}(r::UnitRange{T}) = oftype(r.start + one(T), r.start)
next{T}(r::AbstractUnitRange{T}, i) = (convert(T, i), i + one(T))
done{T}(r::AbstractUnitRange{T}, i) = i == oftype(i, r.stop) + one(T)

start{T}(r::OneTo{T}) = one(T)

# some special cases to favor default Int type to avoid overflow
let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    global start
    global next
    start{T<:smallint}(r::StepRange{T}) = convert(Int, r.start)
    next{T<:smallint}(r::StepRange{T}, i) = (i % T, i + r.step)
    start{T<:smallint}(r::UnitRange{T}) = convert(Int, r.start)
    next{T<:smallint}(r::AbstractUnitRange{T}, i) = (i % T, i + 1)
    start{T<:smallint}(r::OneTo{T}) = 1
end

## indexing

function getindex{T}(v::UnitRange{T}, i::Integer)
    @_inline_meta
    ret = convert(T, first(v) + i - 1)
    @boundscheck ((i > 0) & (ret <= v.stop) & (ret >= v.start)) || throw_boundserror(v, i)
    ret
end

function getindex{T}(v::OneTo{T}, i::Integer)
    @_inline_meta
    @boundscheck ((i > 0) & (i <= v.stop)) || throw_boundserror(v, i)
    convert(T, i)
end

function getindex{T}(v::Range{T}, i::Integer)
    @_inline_meta
    ret = convert(T, first(v) + (i - 1)*step(v))
    ok = ifelse(step(v) > zero(step(v)),
                (ret <= v.stop) & (ret >= v.start),
                (ret <= v.start) & (ret >= v.stop))
    @boundscheck ((i > 0) & ok) || throw_boundserror(v, i)
    ret
end

function getindex{T}(r::FloatRange{T}, i::Integer)
    @_inline_meta
    @boundscheck checkbounds(r, i)
    convert(T, (r.start + (i-1)*r.step)/r.divisor)
end

function getindex{T}(r::LinSpace{T}, i::Integer)
    @_inline_meta
    @boundscheck checkbounds(r, i)
    convert(T, ((r.len-i)*r.start + (i-1)*r.stop)/r.divisor)
end

getindex(r::Range, ::Colon) = copy(r)

function getindex{T<:Integer}(r::AbstractUnitRange, s::AbstractUnitRange{T})
    @_inline_meta
    @boundscheck checkbounds(r, s)
    f = first(r)
    st = oftype(f, f + first(s)-1)
    range(st, length(s))
end

function getindex{T}(r::OneTo{T}, s::OneTo)
    @_inline_meta
    @boundscheck checkbounds(r, s)
    OneTo(T(s.stop))
end

function getindex{T<:Integer}(r::AbstractUnitRange, s::StepRange{T})
    @_inline_meta
    @boundscheck checkbounds(r, s)
    st = oftype(first(r), first(r) + s.start-1)
    range(st, step(s), length(s))
end

function getindex{T<:Integer}(r::StepRange, s::Range{T})
    @_inline_meta
    @boundscheck checkbounds(r, s)
    st = oftype(r.start, r.start + (first(s)-1)*step(r))
    range(st, step(r)*step(s), length(s))
end

function getindex(r::FloatRange, s::OrdinalRange)
    @_inline_meta
    @boundscheck checkbounds(r, s)
    FloatRange(r.start + (first(s)-1)*r.step, step(s)*r.step, length(s), r.divisor)
end

function getindex{T}(r::LinSpace{T}, s::OrdinalRange)
    @_inline_meta
    @boundscheck checkbounds(r, s)
    sl::T = length(s)
    ifirst = first(s)
    ilast = last(s)
    vfirst::T = ((r.len - ifirst) * r.start + (ifirst - 1) * r.stop) / r.divisor
    vlast::T = ((r.len - ilast) * r.start + (ilast - 1) * r.stop) / r.divisor
    return linspace(vfirst, vlast, sl)
end

show(io::IO, r::Range) = print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r)))
show(io::IO, r::UnitRange) = print(io, repr(first(r)), ':', repr(last(r)))
show(io::IO, r::OneTo) = print(io, "Base.OneTo(", r.stop, ")")

=={T<:Range}(r::T, s::T) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
==(r::OrdinalRange, s::OrdinalRange) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
=={T<:LinSpace}(r::T, s::T) = (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s))

function ==(r::Range, s::Range)
    lr = length(r)
    if lr != length(s)
        return false
    end
    u, v = start(r), start(s)
    while !done(r, u)
        x, u = next(r, u)
        y, v = next(s, v)
        if x != y
            return false
        end
    end
    return true
end

intersect(r::OneTo, s::OneTo) = OneTo(min(r.stop,s.stop))

intersect{T1<:Integer, T2<:Integer}(r::AbstractUnitRange{T1}, s::AbstractUnitRange{T2}) = max(first(r),first(s)):min(last(r),last(s))

intersect{T<:Integer}(i::Integer, r::AbstractUnitRange{T}) =
    i < first(r) ? (first(r):i) :
    i > last(r)  ? (i:last(r))  : (i:i)

intersect{T<:Integer}(r::AbstractUnitRange{T}, i::Integer) = intersect(i, r)

function intersect{T1<:Integer, T2<:Integer}(r::AbstractUnitRange{T1}, s::StepRange{T2})
    if isempty(s)
        range(first(r), 0)
    elseif step(s) == 0
        intersect(first(s), r)
    elseif step(s) < 0
        intersect(r, reverse(s))
    else
        sta = first(s)
        ste = step(s)
        sto = last(s)
        lo = first(r)
        hi = last(r)
        i0 = max(sta, lo + mod(sta - lo, ste))
        i1 = min(sto, hi - mod(hi - sta, ste))
        i0:ste:i1
    end
end

function intersect{T1<:Integer, T2<:Integer}(r::StepRange{T1}, s::AbstractUnitRange{T2})
    if step(r) < 0
        reverse(intersect(s, reverse(r)))
    else
        intersect(s, r)
    end
end

function intersect(r::StepRange, s::StepRange)
    if isempty(r) || isempty(s)
        return range(first(r), step(r), 0)
    elseif step(s) < 0
        return intersect(r, reverse(s))
    elseif step(r) < 0
        return reverse(intersect(reverse(r), s))
    end

    start1 = first(r)
    step1 = step(r)
    stop1 = last(r)
    start2 = first(s)
    step2 = step(s)
    stop2 = last(s)
    a = lcm(step1, step2)

    # if a == 0
    #     # One or both ranges have step 0.
    #     if step1 == 0 && step2 == 0
    #         return start1 == start2 ? r : Range(start1, 0, 0)
    #     elseif step1 == 0
    #         return start2 <= start1 <= stop2 && rem(start1 - start2, step2) == 0 ? r : Range(start1, 0, 0)
    #     else
    #         return start1 <= start2 <= stop1 && rem(start2 - start1, step1) == 0 ? (start2:step1:start2) : Range(start1, step1, 0)
    #     end
    # end

    g, x, y = gcdx(step1, step2)

    if rem(start1 - start2, g) != 0
        # Unaligned, no overlap possible.
        return range(start1, a, 0)
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    n = min(stop1 - mod(stop1 - b, a), stop2 - mod(stop2 - b, a))
    m:a:n
end

function intersect(r1::Range, r2::Range, r3::Range, r::Range...)
    i = intersect(intersect(r1, r2), r3)
    for t in r
        i = intersect(i, t)
    end
    i
end

# findin (the index of intersection)
function _findin{T1<:Integer, T2<:Integer}(r::Range{T1}, span::AbstractUnitRange{T2})
    local ifirst
    local ilast
    fspan = first(span)
    lspan = last(span)
    fr = first(r)
    lr = last(r)
    sr = step(r)
    if sr > 0
        ifirst = fr >= fspan ? 1 : ceil(Integer,(fspan-fr)/sr)+1
        ilast = lr <= lspan ? length(r) : length(r) - ceil(Integer,(lr-lspan)/sr)
    elseif sr < 0
        ifirst = fr <= lspan ? 1 : ceil(Integer,(lspan-fr)/sr)+1
        ilast = lr >= fspan ? length(r) : length(r) - ceil(Integer,(lr-fspan)/sr)
    else
        ifirst = fr >= fspan ? 1 : length(r)+1
        ilast = fr <= lspan ? length(r) : 0
    end
    ifirst, ilast
end

function findin{T1<:Integer, T2<:Integer}(r::AbstractUnitRange{T1}, span::AbstractUnitRange{T2})
    ifirst, ilast = _findin(r, span)
    ifirst:ilast
end

function findin{T1<:Integer, T2<:Integer}(r::Range{T1}, span::AbstractUnitRange{T2})
    ifirst, ilast = _findin(r, span)
    ifirst:1:ilast
end

## linear operations on ranges ##

-(r::OrdinalRange) = range(-first(r), -step(r), length(r))
-(r::FloatRange)   = FloatRange(-r.start, -r.step, r.len, r.divisor)
-(r::LinSpace)     = LinSpace(-r.start, -r.stop, r.len, r.divisor)

.+(x::Real, r::AbstractUnitRange) = range(x + first(r), length(r))
.+(x::Real, r::Range) = (x+first(r)):step(r):(x+last(r))
#.+(x::Real, r::StepRange)  = range(x + r.start, r.step, length(r))
.+(x::Real, r::FloatRange) = FloatRange(r.divisor*x + r.start, r.step, r.len, r.divisor)
function .+{T}(x::Real, r::LinSpace{T})
    x2 = x * r.divisor / (r.len - 1)
    LinSpace(x2 + r.start, x2 + r.stop, r.len, r.divisor)
end
.+(r::Range, x::Real)      = x + r
#.+(r::FloatRange, x::Real) = x + r

.-(x::Real, r::Range)      = (x-first(r)):-step(r):(x-last(r))
.-(x::Real, r::FloatRange) = FloatRange(r.divisor*x - r.start, -r.step, r.len, r.divisor)
function .-(x::Real, r::LinSpace)
    x2 = x * r.divisor / (r.len - 1)
    LinSpace(x2 - r.start, x2 - r.stop, r.len, r.divisor)
end
.-(r::AbstractUnitRange, x::Real) = range(first(r)-x, length(r))
.-(r::StepRange , x::Real) = range(r.start-x, r.step, length(r))
.-(r::FloatRange, x::Real) = FloatRange(r.start - r.divisor*x, r.step, r.len, r.divisor)
function .-(r::LinSpace, x::Real)
    x2 = x * r.divisor / (r.len - 1)
    LinSpace(r.start - x2, r.stop - x2, r.len, r.divisor)
end

.*(x::Real, r::OrdinalRange) = range(x*first(r), x*step(r), length(r))
.*(x::Real, r::FloatRange)   = FloatRange(x*r.start, x*r.step, r.len, r.divisor)
.*(x::Real, r::LinSpace)     = LinSpace(x * r.start, x * r.stop, r.len, r.divisor)
.*(r::Range, x::Real)        = x .* r
.*(r::FloatRange, x::Real)   = x .* r
.*(r::LinSpace, x::Real)     = x .* r

./(r::OrdinalRange, x::Real) = range(first(r)/x, step(r)/x, length(r))
./(r::FloatRange, x::Real)   = FloatRange(r.start/x, r.step/x, r.len, r.divisor)
./(r::LinSpace, x::Real)     = LinSpace(r.start / x, r.stop / x, r.len, r.divisor)

promote_rule{T1,T2}(::Type{UnitRange{T1}},::Type{UnitRange{T2}}) =
    UnitRange{promote_type(T1,T2)}
convert{T<:Real}(::Type{UnitRange{T}}, r::UnitRange{T}) = r
convert{T<:Real}(::Type{UnitRange{T}}, r::UnitRange) = UnitRange{T}(r.start, r.stop)

promote_rule{T1,T2}(::Type{OneTo{T1}},::Type{OneTo{T2}}) =
    OneTo{promote_type(T1,T2)}
convert{T<:Real}(::Type{OneTo{T}}, r::OneTo{T}) = r
convert{T<:Real}(::Type{OneTo{T}}, r::OneTo) = OneTo{T}(r.stop)

promote_rule{T1,UR<:AbstractUnitRange}(::Type{UnitRange{T1}}, ::Type{UR}) =
    UnitRange{promote_type(T1,eltype(UR))}
convert{T<:Real}(::Type{UnitRange{T}}, r::AbstractUnitRange) = UnitRange{T}(first(r), last(r))
convert(::Type{UnitRange}, r::AbstractUnitRange) = UnitRange(first(r), last(r))

promote_rule{T1a,T1b,T2a,T2b}(::Type{StepRange{T1a,T1b}},::Type{StepRange{T2a,T2b}}) =
    StepRange{promote_type(T1a,T2a),promote_type(T1b,T2b)}
convert{T1,T2}(::Type{StepRange{T1,T2}}, r::StepRange{T1,T2}) = r

promote_rule{T1a,T1b,UR<:AbstractUnitRange}(::Type{StepRange{T1a,T1b}},::Type{UR}) =
    StepRange{promote_type(T1a,eltype(UR)),promote_type(T1b,eltype(UR))}
convert{T1,T2}(::Type{StepRange{T1,T2}}, r::Range) =
    StepRange{T1,T2}(convert(T1, first(r)), convert(T2, step(r)), convert(T1, last(r)))
convert{T}(::Type{StepRange}, r::AbstractUnitRange{T}) =
    StepRange{T,T}(first(r), step(r), last(r))

promote_rule{T1,T2}(::Type{FloatRange{T1}},::Type{FloatRange{T2}}) =
    FloatRange{promote_type(T1,T2)}
convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::FloatRange{T}) = r
convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::FloatRange) =
    FloatRange{T}(r.start,r.step,r.len,r.divisor)

promote_rule{F,OR<:OrdinalRange}(::Type{FloatRange{F}}, ::Type{OR}) =
    FloatRange{promote_type(F,eltype(OR))}
convert{T<:AbstractFloat}(::Type{FloatRange{T}}, r::OrdinalRange) =
    FloatRange{T}(first(r), step(r), length(r), one(T))
convert{T}(::Type{FloatRange}, r::OrdinalRange{T}) =
    FloatRange{typeof(float(first(r)))}(first(r), step(r), length(r), one(T))

promote_rule{T1,T2}(::Type{LinSpace{T1}},::Type{LinSpace{T2}}) =
    LinSpace{promote_type(T1,T2)}
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace{T}) = r
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::LinSpace) =
    LinSpace{T}(r.start, r.stop, r.len, r.divisor)

promote_rule{F,OR<:OrdinalRange}(::Type{LinSpace{F}}, ::Type{OR}) =
    LinSpace{promote_type(F,eltype(OR))}
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::OrdinalRange) =
    linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r)))
convert{T}(::Type{LinSpace}, r::OrdinalRange{T}) =
    convert(LinSpace{typeof(float(first(r)))}, r)

# Promote FloatRange to LinSpace
promote_rule{F,OR<:FloatRange}(::Type{LinSpace{F}}, ::Type{OR}) =
    LinSpace{promote_type(F,eltype(OR))}
convert{T<:AbstractFloat}(::Type{LinSpace{T}}, r::FloatRange) =
    linspace(convert(T, first(r)), convert(T, last(r)), convert(T, length(r)))
convert{T<:AbstractFloat}(::Type{LinSpace}, r::FloatRange{T}) =
    convert(LinSpace{T}, r)

# +/- of ranges is defined in operators.jl (to be able to use @eval etc.)

## non-linear operations on ranges and fallbacks for non-real numbers ##

.+(x::Number, r::Range) = [ x+y for y=r ]
.+(r::Range, y::Number) = [ x+y for x=r ]

.-(x::Number, r::Range) = [ x-y for y=r ]
.-(r::Range, y::Number) = [ x-y for x=r ]

.*(x::Number, r::Range) = [ x*y for y=r ]
.*(r::Range, y::Number) = [ x*y for x=r ]

./(x::Number, r::Range) = [ x/y for y=r ]
./(r::Range, y::Number) = [ x/y for x=r ]

.^(x::Number, r::Range) = [ x^y for y=r ]
.^(r::Range, y::Number) = [ x^y for x=r ]

## concatenation ##

function vcat{T}(rs::Range{T}...)
    n::Int = 0
    for ra in rs
        n += length(ra)
    end
    a = Array{T}(n)
    i = 1
    for ra in rs, x in ra
        @inbounds a[i] = x
        i += 1
    end
    return a
end

convert{T}(::Type{Array{T,1}}, r::Range{T}) = vcat(r)
collect(r::Range) = vcat(r)

reverse(r::OrdinalRange) = colon(last(r), -step(r), first(r))
reverse(r::FloatRange)   = FloatRange(r.start + (r.len-1)*r.step, -r.step, r.len, r.divisor)
reverse(r::LinSpace)     = LinSpace(r.stop, r.start, r.len, r.divisor)

## sorting ##

issorted(r::AbstractUnitRange) = true
issorted(r::Range) = length(r) <= 1 || step(r) >= zero(step(r))

sort(r::AbstractUnitRange) = r
sort!(r::AbstractUnitRange) = r

sort(r::Range) = issorted(r) ? r : reverse(r)

sortperm(r::AbstractUnitRange) = 1:length(r)
sortperm(r::Range) = issorted(r) ? (1:1:length(r)) : (length(r):-1:1)

function sum{T<:Real}(r::Range{T})
    l = length(r)
    # note that a little care is required to avoid overflow in l*(l-1)/2
    return l * first(r) + (iseven(l) ? (step(r) * (l-1)) * (l>>1)
                                     : (step(r) * l) * ((l-1)>>1))
end

function sum(r::FloatRange)
    l = length(r)
    if iseven(l)
        s = r.step * (l-1) * (l>>1)
    else
        s = (r.step * l) * ((l-1)>>1)
    end
    return (l * r.start + s)/r.divisor
end

function mean{T<:Real}(r::Range{T})
    isempty(r) && throw(ArgumentError("mean of an empty range is undefined"))
    (first(r) + last(r)) / 2
end

median{T<:Real}(r::Range{T}) = mean(r)

function in(x, r::Range)
    n = step(r) == 0 ? 1 : round(Integer,(x-first(r))/step(r))+1
    n >= 1 && n <= length(r) && r[n] == x
end

in{T<:Integer}(x::Integer, r::AbstractUnitRange{T}) = (first(r) <= x) & (x <= last(r))
in{T<:Integer}(x, r::Range{T}) = isinteger(x) && !isempty(r) && x>=minimum(r) && x<=maximum(r) && (mod(convert(T,x),step(r))-mod(first(r),step(r)) == 0)
in(x::Char, r::Range{Char}) = !isempty(r) && x >= minimum(r) && x <= maximum(r) && (mod(Int(x) - Int(first(r)), step(r)) == 0)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
#----------------------------------------------------------------------------# Weight
abstract Weight
abstract BatchWeight <: Weight  # only BatchWeight types work for fitting by batch
abstract StochasticWeight <: BatchWeight

#-----------------------------------------------------------------------# EqualWeight
"""
One of the `Weight` types.  Observations are weighted equally.  For analytical
updates, the online algorithm will give results equal to the offline version.

- `EqualWeight()`
"""
type EqualWeight <: BatchWeight
    nobs::Int
    EqualWeight() = new(0)
    EqualWeight(n::Int) = new(n)
end

#-----------------------------------------------------------------# ExponentialWeight
"""
One of the `Weight` types.  Updates are performed with a constant weight
`λ = 2 / (1 + lookback)`.

- `ExponentialWeight(λ::Float64)`
- `ExponentialWeight(lookback::Int)`
"""
type ExponentialWeight <: Weight
    nobs::Int
    λ::Float64
    function ExponentialWeight(λ::Real = 1.0)
            @assert 0 <= λ <= 1
            new(0, λ)
    end
    ExponentialWeight(lookback::Integer) = ExponentialWeight(2.0 / (lookback + 1))
end

#----------------------------------------------------------------# BoundedEqualWeight
"""
One of the `Weight` types.  Uses `EqualWeight` until reaching `λ = 2 / (1 + lookback)`,
then weights are held constant.

- `BoundedEqualWeight(λ::Float64)`
- `BoundedEqualWeight(lookback::Int)`
"""
type BoundedEqualWeight <: Weight
    nobs::Int
    λ::Float64
    function BoundedEqualWeight(λ::Real = 1.0)
            @assert 0 <= λ <= 1
            new(0, λ)
    end
    BoundedEqualWeight(lookback::Integer) = BoundedEqualWeight(2.0 / (lookback + 1))
end

#----------------------------------------------------------------------# LearningRate
"""
One of the `Weight` types.  It's primary use is for the OnlineStats that use stochastic
approximation (`StatLearn`, `QuantReg`, `QuantileMM`, `QuantileSGD`, `NormalMix`, and
`KMeans`).  The weight at update `t` is `1 / t ^ r`.  When weights reach `λ`, they are
held consant.  Compare to `LearningRate2`.

- `LearningRate(r = 0.5, λ = 0.0)`
"""
type LearningRate <: StochasticWeight
    nobs::Int
    nups::Int
    r::Float64
    λ::Float64
    function LearningRate(r::Real = 0.5, λ::Real = 0.0)
        @assert 0 < r <= 1
        @assert λ >= 0
        new(0, 0, r, λ)
    end
end

#---------------------------------------------------------------------# LearningRate2
"""
One of the `Weight` types.  It's primary use is for the OnlineStats that use stochastic
approximation (`StatLearn`, `QuantReg`, `QuantileMM`, `QuantileSGD`, `NormalMix`, and
`KMeans`).  The weight at update `t` is `1 / (1 + c * (t - 1))`.  When weights reach
`λ`, they are held consant.  Compare to `LearningRate`.

- `LearningRate2(c = 0.5, λ = 0.0)`
"""
type LearningRate2 <: StochasticWeight
    nobs::Int
    nups::Int
    c::Float64
    λ::Float64
    function LearningRate2(c::Real = 0.5, λ = 0.0)
        @assert λ >= 0
        @assert c > 0
        new(0, 0, c, λ)
    end
end

#---------------------------------------------------------------------------# methods
StatsBase.nobs(w::Weight) = w.nobs
nups(w::StochasticWeight) = w.nups

# increase nobs by the number of new observations
updatecounter!(w::Weight, n2::Int = 1)              = (w.nobs += n2)
updatecounter!(w::StochasticWeight, n2::Int = 1)    = (w.nobs += n2; w.nups += 1)
updatecounter!(o::OnlineStat, n2::Int = 1) = updatecounter!(o.weight, n2)

# Get weight for update of size n2, typically immediately after updatecounter!
weight(w::EqualWeight, n2::Int = 1)         = n2 / w.nobs
weight(w::BoundedEqualWeight, n2::Int = 1)  = nobs(w) == 0 ? 1.0 : max(w.λ, n2 / w.nobs)
weight(w::ExponentialWeight, n2::Int = 1)   = w.λ
weight(w::LearningRate, n2::Int = 1)        = max(w.λ, exp(-w.r * log(w.nups)))
weight(w::LearningRate2, n2::Int = 1)       = max(w.λ, 1.0 / (1.0 + w.c * (w.nups - 1)))
weight(o::OnlineStat, n2::Int = 1) = weight(o.weight, n2)

ext: jl


text:
# ---
# notebook: false
# ---
y = 1 + 2


text:
# This wouldn't be possible without this https://github.com/R3l3ntl3ss/Meme_Api.
# Also, thank you to Twitter for approving my developer account request in a day. You rock!

using HTTP, JSON, Base64
include("authorizn.jl")

API = "https://meme-api.herokuapp.com/"

post_url = "https://api.twitter.com/1.1/statuses/update.json"

upload_url = "https://upload.twitter.com/1.1/media/upload.json"

function step_one(name)
    image_data = ""
    image_size = 0
    try
        open(name) do file
            image_data = read(file)
            image_size = stat(file).size
        end
    catch
        println("not found...")
        exit(1)
    end
    s1 = Dict("command" => "INIT","total_bytes" => image_size,"media_type" => "image/jpeg","media_category" => "tweet_image")
    bot_init = tweet_oauth(upload_url, s1)
    json = JSON.parse(String(bot_init.body))
    return step_two(image_data, json["media_id"])
end

function step_two(image_data, media_id)
    s2 = Dict("command" => "APPEND", "media_id" => media_id, "media" => base64encode(image_data),"segment_index" => 0)
    r_append = tweet_oauth(upload_url, s2)
    return step_three(media_id)
end

function step_three(media_id)
    s3 = Dict("command" => "FINALIZE", "media_id" => media_id)
    r_final = tweet_oauth(upload_url, s3)
    return media_id
end

function get_meme()
    subs = ["PewdiepieSubmissions","ProgrammerHumor","mathmemes"]
    randomsub = rand(1 : 3, 1)[1]
    url = API*"gimme/"*subs[randomsub]*"/1"
    req = HTTP.request("GET", url)
    json = JSON.parse(String(req.body))
    imagel_ink = json["memes"][1]["url"]
    type = split(imagel_ink, ".")[end]
    img_name = "image.$type"
    println("status: ", req.status)
    download(imagel_ink, img_name)
    return json, img_name
end

function tweet()
    meme, name = get_meme()
    post_link = meme["memes"][1]["postLink"]
    title = meme["memes"][1]["title"]
    media_id = step_one(name)
    options = Dict("status" => "$title\nSource: $post_link\n\n-By http://bit.ly/PseudoCodeNerd-memeBot", "media_ids" => media_id)
    r_tweet = tweet_oauth(post_url, options)
    println("(Testing) Status: ", r_tweet.status)
end

#Finally!
tweet()

max_stars_repo_name: PseudoCodeNerd/julia-memebot-twitter

max_issues_repo_name: PseudoCodeNerd/julia-memebot-twitter

max_forks_repo_name: PseudoCodeNerd/julia-memebot-twitter


text:
using AxisArrays

windowing(x,dim=timedim(x);kwds...) = windowing(hastimes(x),x,dim;kwds...)
map_windowing(fn,x,dim=timedim(x);kwds...) =
  map_windowing(fn,hastimes(x),x,dim;kwds...)

function map_windowing(fn,::HasTimes,x,dim;step=nothing,flatten=false,kwds...)
  windows = windowing(x,dim;step=step,kwds...)
  xs = map(windows) do ixs
    fn(x[Axis{:time}(ixs)])
  end

  if flatten
    result = cat(xs...,dims = ndims(xs[1])+1)
    axisnames = Symbol.("ax".*string.(1:ndims(result)-1))
    AxisArray(result,
              (Axis{ax}(1:n) for (ax,n) in zip(axisnames,
                                               size(result)[1:end-1]))...,
              AxisArrays.axes(windows,Axis{:time}))
  else
    AxisArray(xs,AxisArrays.axes(windows,Axis{:time}))
  end
end
map_windowing(fn,::HasNoTimes,x,dim;kwds...) =
  map(ixs -> fn(x[Axis{:time}(ixs)]),windowing(HasNoTimes(),x,dim;kwds...))

function windowing(::HasNoTimes,x::AbstractArray,dim;
                   length=nothing,step=nothing,minlength=length)
  (max(1,t-length+1):t for t in Base.axes(x,dim)[min(minlength,end):step:end])
end

function windowing(::HasTimes,data::AbstractArray,dim;
                   length=nothing,step=nothing,minlength=length)
  helper(x::Number) = max(1,floor(Int,x*s / Δt(data)))
  helper(x::Quantity) = max(1,floor(Int,x / Δt(data)))
  length_,step_,minlength_ = helper.((length,step,minlength))

  win = windowing(HasNoTimes(),data,dim,
                  length=length_,step=step_,minlength=minlength_)
  AxisArray(collect(win),
            Axis{:time}((min(minlength_,size(data,dim)):step_:size(data,dim)).*
                        Δt(data) .- (length_/2*Δt(data))))
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE

# # Use of heat exchange boundary condition
#
# Boundary conditions needs fields `heat transfer coefficient` and
# `external temperature` defined. Boundary condition is of type
# ```math
# f = h\left(T-T_{\mathrm{u}}\right)
# ```

using FEMBase, Test, HeatTransfer, SparseArrays

# Set up initial data

X = Dict(
    1 => [0.0,0.0],
    2 => [1.0,0.0])

T = Dict(
    1 => 0.0,
    2 => 0.0)

# Create element and update fields

element = Element(Seg2, [1, 2])
update!(element, "geometry", X)
update!(element, "temperature", 0.0 => T)
update!(element, "heat transfer coefficient", 6.0)
update!(element, "external temperature", 1.0)

# Create problem, add elements to problem and assemble at time ``t=0``:

problem = Problem(PlaneHeat, "test problem", 1)
add_elements!(problem, [element])
assemble!(problem, 0.0)

# Test for stiffness matrix ``\boldsymbol{K}`` and force vector ``\boldsymbol{f}``:

K = Matrix(sparse(problem.assembly.K))
f = Vector(sparse(problem.assembly.f)[:])
@test isapprox(K, [2.0 1.0; 1.0 2.0])
@test isapprox(f, [3.0; 3.0])


text:
module SpecialArrays

using Adapt

using Base: @propagate_inbounds, @pure, @_inline_meta
using Base: require_one_based_indexing, tail, unsafe_length

using Compat
using DocStringExtensions
using Requires: @require
using UnsafeArrays

const NestedArray{V,M,N} = AbstractArray{<:AbstractArray{V,M},N}

include("tuple.jl")
include("util.jl")

export inner_eltype, inner_ndims, inner_axes, inner_size, inner_length
include("functions.jl")

export SlicedArray, slice, align
include("slicedarray.jl")

export FlattenedArray, flatten
include("flattenedarray.jl")

export BatchedVector, batch, similarbatch
include("batchedvector.jl")

function __init__()
    @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" include("zygote.jl")
end

end # module


text:
using DependencyTrees
using DependencyTrees.TransitionParsing
using Documenter

DocMeta.setdocmeta!(DependencyTrees, :DocTestSetup, :(using DependencyTrees); recursive=true)

makedocs(
    sitename = "DependencyTrees.jl",
    format = Documenter.HTML(),
    modules = [DependencyTrees],
    pages = [
        "Home" => "index.md",
        "Trees" => ["trees.md", "treebanks.md"],
        "Transition Parsing" => "transition_parsing.md",
        "Evaluation" => "evaluation.md",
        "Errors" => "errors.md"
    ],
    doctest = true
)

deploydocs(repo = "github.com/dellison/DependencyTrees.jl.git")


text:
export read_dataset

function read_dataset(
  dname::AbstractString,  # name of the dataset to load
  nvectors::Union{Integer, UnitRange},      # number of vectors to read
  V::Bool=false )         # whether to print progress

  if V println("Loading $(dname)... "); end

  # === Deep_babenko ===
  if dname == "Deep1M_babenko"

    fname = "./data/deep_babenko/deep1M_learn.fvecs"
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "Deep1M_babenko_groundtruth"

    fname = "./data/deep_babenko/deep1M_groundtruth.ivecs"
    X     = ivecs_read(nvectors, fname)
    return X

  elseif dname == "Deep1M_babenko_query"

    fname = "./data/deep_babenko/deep1M_queries.fvecs"
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "Deep1M_babenko_base"

    fname = "./data/deep_babenko/deep1M_base.fvecs"
    @show nvectors, fname
    X     = fvecs_read(nvectors, fname)
    return X

  # === GIST1M ===
  elseif dname == "GIST1M"

    fname = "./data/gist/gist_learn.fvecs";
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "GIST1M_query"

    fname = "./data/gist/gist_query.fvecs";
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "GIST1M_groundtruth"

    fname = "./data/gist/gist_groundtruth.ivecs";
    X     = ivecs_read(nvectors, fname)
    return X

  elseif dname == "GIST1M_base"

    fname = "./data/gist/gist_base.fvecs";
    X     = fvecs_read(nvectors, fname)
    return X

  # ====  SIFT1M ===
  elseif dname == "SIFT1M"

    fname = "./data/sift/sift_learn.fvecs";
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "SIFT1M_query"

    fname = "./data/sift/sift_query.fvecs"
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "SIFT1M_groundtruth"

    fname = "./data/sift/sift_groundtruth.ivecs"
    X     = ivecs_read(nvectors, fname)
    return X

  elseif dname == "SIFT1M_base"

    fname = "./data/sift/sift_base.fvecs";
    X     = fvecs_read(nvectors, fname)
    return X

  # === Convnet1M ===
  elseif dname == "Convnet1M_base"

    fname = "./data/feats/feats_m_128.mat"
    X     = h5read(fname, "feats_m_128_base");

  elseif dname == "Convnet1M"

    fname = "./data/feats/feats_m_128.mat"
    X     = h5read(fname, "feats_m_128_train");

  elseif dname == "Convnet1M_query"

    fname = "./data/feats/feats_m_128.mat"
    X     = h5read(fname, "feats_m_128_test");

  elseif dname == "Convnet1M_groundtruth"

    fname = "./data/feats/feats_m_128_gt.mat"
    X     = h5read(fname, "gt");

  # === Deep1M ===
  elseif dname == "Deep1M_base"

    fname = "./data/deep/deep.h5"
    X     = h5read(fname, "base");

  elseif dname == "Deep1M"

    fname = "./data/deep/deep.h5"
    X     = h5read(fname, "train");

  elseif dname == "Deep1M_query"

    fname = "./data/deep/deep.h5"
    X     = h5read(fname, "query");

  elseif dname == "Deep1M_groundtruth"

    fname = "./data/deep/deep.h5"
    X     = h5read(fname, "gt");
    return X

  # === Deep1B ===
  elseif dname == "Deep1B"

    fname = "./data/deep1b/learn_00"
    X     = fvecs_read(nvectors, fname)
    return X

  elseif dname == "Deep1B_base"

    # address on curiosity
    fname = "/scratch/julm/deep1b/base_all"
    #fname = "./data/deep1b/base_all"
    X     = fvecs_read(nvectors, fname)
    return X

  # === SIFT1B
  elseif dname == "SIFT1B_query"

    fname = "/hdd/sift1b/bigann_query.bvecs"
    X     = bvecs_read(nvectors,fname)
    X = convert(Matrix{Float32},X)
    return X

  elseif dname == "SIFT1B_base"

    # fname = "./data/sift1b/bigann_base.bvecs"
    fname = "/hdd/sift1b/bigann_base.bvecs"
    # fname = "/scratch/zakhmi/bigann_base.bvecs" # curiosity
    X     = bvecs_read(nvectors,fname)
    X = convert(Matrix{Float32},X)
    return X

  elseif dname == "SIFT1B" || dname == "SIFT10M"

    fname = "/hdd/sift1b/bigann_learn.bvecs"
    X     = bvecs_read(nvectors,fname)
    X = convert(Matrix{Float32},X)
    return X

  elseif dname == "SIFT1B_groundtruth"

    fname = "/hdd/sift1b/gnd/idx_1000M.ivecs"
    X     = ivecs_read(nvectors, fname)
    return X

  elseif dname == "SIFT1B_groundtruth_10M"

    fname = "/hdd/sift1b/gnd/idx_10M.ivecs"
    X     = ivecs_read(nvectors, fname)
    return X

  elseif dname == "SIFT1B_groundtruth_1M"

    fname = "/hdd/sift1b/gnd/idx_1M.ivecs"
    X     = ivecs_read(nvectors, fname)
    return X

  # === MNIST
  elseif dname == "MNIST_query"

    fname = "./data/mnist/mnist.h5";
    X     = h5read(fname, "test");
    X     = reshape( X, 28 * 28, 10000 );
    X     = convert(Matrix{Float32}, X);

  elseif dname == "MNIST" || dname == "MNIST_base"

    fname = "./data/mnist/mnist.h5";
    X     = h5read(fname, "train");
    X     = reshape( X, 28 * 28, 60000 );
    X     = convert(Matrix{Float32}, X);

  elseif dname == "MNIST_groundtruth"

    fname = "./data/mnist/mnist.h5";
    X     = h5read(fname, "gt");

  # === labelme22K
  elseif dname == "labelme_query"

    fname = "./data/labelme/label.h5";
    X     = h5read(fname, "query");
    X     = convert( Matrix{Float32}, X );

  elseif dname == "labelme" || dname == "labelme_base"

    fname = "./data/labelme/label.h5";
    X     = h5read(fname, "train");
    X     = convert( Matrix{Float32}, X );

  elseif dname == "labelme_groundtruth"

    fname = "./data/labelme/label.h5";
    X     = h5read(fname, "gt");

  else

    error("Dataset $(dname) unknown")

  end

  _, n = size( X );

  if typeof( nvectors ) <: Integer
    if nvectors <= n
      X = X[:, 1:nvectors];
    else
      error("Asked to read $nvectors vectors, but the datasets has only $n vectors.");
    end
  else
    X = X[:, nvectors];
  end

  return X
end


text:
############################################################
# joBlock ##############################################
############################################################

export joBlock, joBlockException

struct joBlockException <: Exception
    msg :: String
end

############################################################
## outer constructors

"""
    julia> op = joBlock(rows,op1[,op2][,...];[weights=...,][name=...])

Block operator composed from different square JOLI operators

# Signature

    joBlock(rows::Vector{RVDT},ops::joAbstractLinearOperator...;
        weights::LocalVector{WDT}=zeros(0),name::String="joBlock")
            where {RVDT<:Integer,WDT<:Number}

# Arguments

- `rows`: # number of operator in a blocks' row
- `op#`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `weights`: vector of waights for each operator
    - `name`: custom name

# Notes

- operators are to be given in row-major order
- all operators in a blocs' row must have the same # of rows (M)
- sum of Ns for operators in each row must be the same
- all given operators must have same domain/range types
- the domain/range types of joBlock are equal to domain/range types of the given operators

# Example

define operators

    a=rand(ComplexF64,4,4);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")
    b=rand(ComplexF64,4,8);
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64,name="B")
    c=rand(ComplexF64,6,6);
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64,name="C")
    d=rand(ComplexF64,6,6);
    D=joMatrix(d;DDT=ComplexF32,RDT=ComplexF64,name="D")

define weights if needed

    w=rand(ComplexF64,4)

basic blocks in function syntax

    S=joBlock([2,2],A,B,C,D)

basic blocks in [] syntax

    S=[A B; C D]

weighted blocks

    S=joBlock([2,2],A,B,C,D;weights=w)

"""
function joBlock(rows::Vector{RVDT},ops::joAbstractLinearOperator...;
           weights::LocalVector{WDT}=zeros(0),name::String="joBlock") where {RVDT<:Integer,WDT<:Number}

    isempty(ops) && throw(joBlockException("empty argument list"))
    l=length(ops)
    sum(rows)==l || throw(joBlockException("sum of operators in the rows does not match # of operators"))
    r=length(rows)
    for i=1:l
        deltype(ops[i])==deltype(ops[1]) || throw(joBlockException("domain type mismatch for $i operator"))
        reltype(ops[i])==reltype(ops[1]) || throw(joBlockException("range type mismatch for $i operator"))
    end
    i=0
    mm=1
    for j=1:r
        for k=1:rows[j]
            i+=1
            ops[i].m==ops[mm].m || throw(joBlockException("size mismatch (M) for $k operator in row $j"))
        end
        mm+=rows[j]
    end
    (length(weights)==l || length(weights)==0) || throw(joBlockException("lenght of weights vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    rm=zeros(Int,r)
    rn=zeros(Int,r)
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    mo=zeros(Int,l)
    no=zeros(Int,l)
    i=0
    mm=1
    ro=0
    for j=1:r
        co=0
        rm[j]=ops[mm].m
        for k=1:rows[j]
            i+=1
            ms[i]=ops[mm].m
            ns[i]=ops[i].n
            mo[i]=ro
            no[i]=co
            rn[j]+=ops[i].n
            co+=ops[i].n
        end
        rn[j]==rn[1] || throw(joBlockException("size mismatch (N) for operators in row $j"))
        ro+=ops[mm].m
        mm+=rows[j]
    end
    m=sum(rm)
    n=rn[1]
    weighted=(length(ws)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*ops[i])
            push!(fops_T,ws[i]*transpose(ops[i]))
            push!(fops_A,conj(ws[i])*adjoint(ops[i]))
            push!(fops_C,conj(ws[i])*conj(ops[i]))
        else
            push!(fops,ops[i])
            push!(fops_T,transpose(ops[i]))
            push!(fops_A,adjoint(ops[i]))
            push!(fops_C,conj(ops[i]))
        end
    end
    return joCoreBlock{deltype(fops[1]),reltype(fops[1])}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end


text:
###########################################################################
#
#  Independent Metropolis-Hastings (IMH)
#
#  Parameters:
#    -logCandidate: candidate density
#    -randCandidate: function for sampling from candidate density
#
###########################################################################

export IMH

###########################################################################
#                  IMH type
###########################################################################

immutable IMH <: MCMCSampler
  logCandidate::Function
  randCandidate!::Function
end

IMH(proposal::ContinuousMultivariateDistribution) =
  IMH(pars::Vector{Float64} -> logpdf(proposal, pars), sample::Vector{Float64} -> rand!(proposal, sample))

###########################################################################
#                  IMH task
###########################################################################

# IMH sampling
function SamplerTask(model::MCMCModel, sampler::IMH, runner::MCMCRunner)
	local pars, proposedPars
	local logTarget, logCandidate, proposedLogTarget, proposedLogCandidate

	# hook inside Task to allow remote resetting
	task_local_storage(:reset, (resetPars::Vector{Float64}) -> (pars = copy(resetPars); logTarget = model.eval(pars)))

	# initialization
	pars, proposedPars = copy(model.init), zeros(model.size)
	logTarget, logCandidate = model.eval(pars), sampler.logCandidate(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"
  @assert isfinite(logCandidate) "Initial values out of candidate density support, try other values"

	# main loop
	while true
		sampler.randCandidate!(proposedPars)
		proposedLogTarget, proposedLogCandidate = model.eval(proposedPars), sampler.logCandidate(proposedPars)

	  ratio = proposedLogTarget-logTarget-proposedLogCandidate+logCandidate
		if ratio > 0 || (ratio > log(rand()))
			ms = MCMCSample(proposedPars, proposedLogTarget, pars, logTarget, {"accept" => true})
	    produce(ms)
	    pars, logTarget, logCandidate = copy(proposedPars), copy(proposedLogTarget), copy(proposedLogCandidate)
	  else
			ms = MCMCSample(pars, logTarget, pars, logTarget, {"accept" => false})
      produce(ms)
    end
	end
end


text:
using Manopt, Documenter, Literate

# generate examples using Literate
tutorialsInputPath = joinpath(@__DIR__, "..", "src/tutorials")
tutorialsRelativePath = "tutorials/"
tutorialsOutputPath = joinpath(@__DIR__, "src/" * tutorialsRelativePath)
tutorials = ["MeanAndMedian", "GradientOfSecondOrderDifference", "JacobiFields"]
menuEntries = ["Getting Started: Optimize!", "Gradient of \$d_2\$", "Jacobi Fields"]
TutorialMenu = Array{Pair{String,String},1}()
for (i, tutorial) in enumerate(tutorials)
    global TutorialMenu
    sourceFile = joinpath(tutorialsInputPath, tutorial * ".jl")
    targetFile = joinpath(tutorialsOutputPath, tutorial * "md")
    Literate.markdown(
        sourceFile,
        tutorialsOutputPath;
        name = tutorial,
        # codefence = "```julia" => "```",
        credit = false,
    )
    push!(TutorialMenu, menuEntries[i] => joinpath(tutorialsRelativePath, tutorial * ".md"))
end
makedocs(
    format = Documenter.HTML(prettyurls = false),
    modules = [Manopt],
    sitename = "Manopt.jl",
    pages = [
        "Home" => "index.md",
        "About" => "about.md",
        "Plans" => "plans/index.md",
        "Solvers" => [
            "Introduction" => "solvers/index.md",
            "Conjugate gradient descent" => "solvers/conjugate_gradient_descent.md",
            "Cyclic Proximal Point" => "solvers/cyclic_proximal_point.md",
            "Douglas–Rachford" => "solvers/DouglasRachford.md",
            "Gradient Descent" => "solvers/gradientDescent.md",
            "Nelder–Mead" => "solvers/NelderMead.md",
            "Subgradient method" => "solvers/subgradient.md",
            "Steihaug-Toint TCG Method" => "solvers/truncated_conjugate_gradient_descent.md",
            "Riemannian Trust-Regions Solver" => "solvers/trust_regions.md",
        ],
        "Functions" => [
            "Introduction" => "functions/index.md",
            "cost functions" => "functions/costFunctions.md",
            "Differentials" => "functions/differentials.md",
            "Adjoint Differentials" => "functions/adjointDifferentials.md",
            "Gradients" => "functions/gradients.md",
            "JacobiFields" => "functions/jacobiFields.md",
            "Proximal Maps" => "functions/proximalMaps.md",
            "Specific manifold functions" => "functions/manifold.md",
        ],
        "Helpers" => [
            "Data" => "helpers/data.md",
            "Error Measures" => "helpers/errorMeasures.md",
            "Exports" => "helpers/exports.md",
        ],
        "Tutorials" => TutorialMenu,
        "Function Index" => "list.md",
    ],
)
deploydocs(
    repo = "github.com/JuliaManifolds/Manopt.jl",
    push_preview = true
)


text:
using Plots, LightGraphs, JLD, LaTeXStrings, MTSG, NGWP, Distances
using Plots.PlotMeasures

G = loadgraph("../datasets/RGC100.lgz"); N = nv(G)
X = load("../datasets/RGC100_xyz.jld", "xyz")[:, 1:2]
X3 = load("../datasets/RGC100_xyz.jld", "xyz")
L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L)
standardize_eigenvectors!(𝚽)

dist_X = pairwise(Euclidean(1e-12), X3; dims = 1)
A = 1.0 .* adjacency_matrix(G)
W = zeros(N, N); W[A .> 0] = 1 ./ dist_X[A .> 0]; W = A .* W

Q = incidence_matrix(G; oriented = true)
∇𝚽 = Q' * 𝚽
edge_length = sqrt.(sum((Q' * X3).^2, dims = 2)[:])

G_Sig = GraphSig(W)
GP = partition_tree_fiedler(G_Sig; swapRegion = false)


text:
"""
`module  JAC.TestFrames`
    ... a submodel of JAC that contains all methods for testing individual methods, modules, ....
"""
module TestFrames

    using Printf, JAC, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..TableStrings

    export testDummy

    function testCompareFiles(fold::String, fnew::String, sa::String, noLines::Int64)
        success = true
        printTest, iostream = Defaults.getDefaults("test flag/stream")

        # Compare the test computations with previous results
        oldLines = readlines(fold)
        newLines = readlines(fnew)
        iold = 0;  for i=1:length(oldLines)   line = oldLines[i];   if  occursin(sa, line)  iold = i;   break   end   end
        inew = 0;  for i=1:length(newLines)   line = newLines[i];   if  occursin(sa, line)  inew = i;   break   end   end
        if  iold == 0   ||   inew == 0    success = false
            if printTest   println(iostream, "Tries to compare two inappropriate files fold = $(fold); fnew =$(fnew) on string $sa  ($iold, $inew)")  end
            return( success )
        end
        #
        ii = inew + 1
        for  i = iold+2:iold+noLines
        ii = ii + 1
        if   length(oldLines[i]) < 5    continue    end
        if   oldLines[i] != newLines[ii]    success = false
            if  printTest  println(iostream, "    *** Old::  " * oldLines[i])
                            println(iostream, "    *** New::  " * newLines[ii])
            end
        end
        end

        return( success )
    end

    function testPrint(sa::String, success::Bool)
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        ok(succ) =  succ ? "[OK]" : "[Fail]"
        sb = sa * TableStrings.hBlank(110);   sb = sb[1:100] * ok(success);    println(iostream, sb)
        return( nothing )
    end

    """
    `Basics.testMethod_integrate_ongrid(; short::Bool=true)`  ... tests the integration on grid.
    """
    function testMethod_integrate_ongrid(; short::Bool=true)
        success = true
        printTest, iostream = Defaults.getDefaults("test flag/stream")

        # Test the integration on the grid for an analytical function
        grid = Radial.Grid("grid: by given parameters", rnt=2.0e-10, h=5.0e-3, NoPoints=9000)
        function f1(r)
            A = 10^5;    gamma = 3;    a = 300
            return( A * r^gamma * exp(-a * r^2) )
        end
        exact1 = 5/9

        integrand = Vector{Float64}(size(grid.r, 1))
        for i = 1:size(grid.r, 1)
            integrand[i] = f1(grid.r[i])
        end
        integrand = integrand .* grid.rp[1:size(integrand, 1)]   # adopt to the form of Grasp92

        integral  = Basic.integrate("function: on radial grid, Newton-Cotes", integrand, grid)
        err = abs(integral - exact1)
        if  abs(err) > 1.0e-12
            success = false
            if printTest   info(iostream, "... Newton-Cotes:  I = $integral,  Err = $err")  end
        end

        integral  = Basic.integrate("function: on radial grid, Simpson rule", integrand, grid)
        err = abs(integral - exact1)
        if  abs(err) > 1.0e-12
            success = false
            if printTest   info(iostream, "... Simpson rule:  I = $integral,  Err = $err")  end
        end

        integral  = Basic.integrate("function: on radial grid, trapez rule", integrand, grid)
        err = abs(integral - exact1)
        if  abs(err) > 1.0e-12
            success = false
            if printTest   info(iostream, "... trapez rule:  I = $integral,  Err = $err")  end
        end

        println(iostream, "Warning(testMethod_integrate_ongrid): test of integration with Grasp orbitals has been set silent.")
        #=
        # Test the integration with orbital functions from Grasp92
        grid = Radial.Grid("grid: exponential")

        orbitals1 = Basics.readOrbitalFileGrasp92("../test/approved/Ne-0+-scf.exp.out", grid)
        orbitals2 = Basics.readOrbitalFileGrasp92("../test/approved/Ne-1+-scf.exp.out", grid)

        for i = 1:size(orbitals1, 1)
            for j = 1:size(orbitals2, 1)
                orb1      = orbitals1[i];   orb2 = orbitals2[j]
                if   orb1.subshell.kappa != orb2.subshell.kappa    break   end

                mtp       = min(size(orb1.P, 1), size(orb2.P, 1))
                integrand = ( orb1.P[1:mtp] .* orb2.P[1:mtp] + orb1.Q[1:mtp] .* orb2.Q[1:mtp] ) .* grid.rp[1:mtp]

                integrala = Basic.integrate("function: on radial grid, Newton-Cotes", integrand, grid)^2
                integralb = Basic.integrate("function: on radial grid, Simpson rule", integrand, grid)^2
                integralc = Basic.integrate("function: on radial grid, trapez rule",  integrand, grid)^2

                info(iostream, "<$(string(orb1.subshell)) | $(string(orb2.subshell))> = $integrala, $integralb, $integralc")
            end
        end  =#

        testPrint("testMethod_integrate_ongrid()::", success)
        return(success)
    end

    """
    `TestFrames.testMethod_Wigner_3j(; short::Bool=true)`  ... tests on Wigner 3j symbols.
    """
    function testMethod_Wigner_3j(; short::Bool=true)
        success = true
        printTest, iostream = Defaults.getDefaults("test flag/stream")

        # Wigner_3j(1,2,1,0,0,0) = 0.36514837167011074230
        a = c = AngularJ64(1);    b = AngularJ64(2);    ma = mb = mc = AngularM64(0)
        wa = AngularMomentum.Wigner_3j(a,b,c,ma,mb,mc)
        if  abs(wa - 0.36514837167011074230) > 1.0e-12
            success = false
            if printTest   info(iostream, "Wigner_3j(1,2,1,0,0,0) = 0.36514837167011074230 ... but obtains value = $wa")   end
        end

        # Wigner_3j(3,6,3,0,0,0) = 0.182482967150452976281
        a = c = AngularJ64(3);    b = AngularJ64(6);    ma = mb = mc = AngularM64(0)
        wa = AngularMomentum.Wigner_3j(a,b,c,ma,mb,mc)
        if  abs(wa - 0.18248296715045297628) > 1.0e-12
            success = false
            if printTest   info(iostream, "Wigner_3j(3,6,3,0,0,0) = 0.182482967150452976281 ... but obtains value = $wa")   end
        end

        # Wigner_3j(1/2,1,1/2,1/2,0,-1/2) = 0.40824829046386301637
        a = c = AngularJ64(1//2);    b = AngularJ64(1);    ma =  AngularM64(1//2);   mb =  AngularM64(0);    mc = AngularM64(-1//2)
        wa = AngularMomentum.Wigner_3j(a,b,c,ma,mb,mc)
        if  abs(wa - 0.40824829046386301637) > 1.0e-12
            success = false
            if printTest   info(iostream, "Wigner_3j(1/2,1,1/2,1/2,0,-1/2) = 0.40824829046386301637 ... but obtains value = $wa")  end
        end

        testPrint("testMethod_Wigner_3j()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_AlphaVariation(; short::Bool=true)`  ... tests on module AlphaVariation.
    """
    function testModule_AlphaVariation(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-AlphaVariation-new.sum")
        printstyled("\n\nTest the module  AlphaVariation  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Model(26.); properties=[JAC.AlphaX],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                alphaSettings=AlphaVariation.Settings(true, true, false, Int64[]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-AlphaVariation-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-AlphaVariation-new.sum"), "Alpha variation parameters:", 1)
        testPrint("testModule_AlphaVariation()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_AutoIonization(; short::Bool=true)`  ... tests on module AutoIonization.
    """
    function testModule_AutoIonization(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-AutoIonization-new.sum")
        printstyled("\n\nTest the module  AutoIonization  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Model(36.); grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 1.5e-2, NoPoints = 600),
                                initialConfigs=[Configuration("1s^2 2s^2 2p"), Configuration("1s 2s^2 2p^2")],
                                finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2p^2")], process = JAC.Auger,
                                processSettings = AutoIonization.Settings(true, true, true, Tuple{Int64,Int64}[(3,1), (4,1), (5,1), (6,1)], 0., 1.0e6, 2, "Coulomb") )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-AutoIonization-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-AutoIonization-new.sum"), "AutoIonization rates and intr", 17)
        testPrint("testModule_AutoIonization()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_CoulombExcitation(; short::Bool=true)`  ... tests on module CoulombExcitation.
    """
    function testModule_CoulombExcitation(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-CoulombExcitation-new.sum")
        printstyled("\n\nTest the module  CoulombExcitation  ... \n", color=:cyan)
        ### Make the tests
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        println(iostream, "Make the comparison with approved data for ... test-CoulombExcitation-new.sum")
        success = true
        ## success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-AutoIonization-approved.sum"),
        ##                             joinpath(@__DIR__, "..", "test", "test-AutoIonization-new.sum"), "AutoIonization rates and intrinsic angular parameters:", 25)
        testPrint("testModule_CoulombExcitation()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_DecayYield(; short::Bool=true)`  ... tests on module DecayYield.
    """
    function testModule_DecayYield(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-DecayYield-new.sum")
        printstyled("\n\nTest the module  DecayYield  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(12.); properties=[JAC.Yields],
                                grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600),
                                configs=[Configuration("1s 2s^2 2p^6")],
                                yieldSettings=DecayYield.Settings("SCA", true, false, Int64[]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-DecayYield-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-DecayYield-new.sum"), "Auger rates (without", 6)
        testPrint("testModule_DecayYield()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_Dielectronic(; short::Bool=true)`  ... tests on module Dielectronic.
    """
    function testModule_Dielectronic(; short::Bool=true)     Defaults.setDefaults("print summary: open", "test-Dielectronic-new.sum")
        printstyled("\n\nTest the module  Dielectronic  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600),
                                initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                                intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                                finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                                process = JAC.Dierec,
                                processSettings=Dielectronic.Settings([E1, M1], [UseCoulomb, UseBabushkin], false,
                                                                    true, Tuple{Int64,Int64,Int64}[(1,1,0)], 0., 0., 0., "Coulomb")  )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-Dielectronic-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-Dielectronic-new.sum"),
                                "Partial (Auger) capture", 10)
        testPrint("testModule_Dielectronic()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_Einstein(; short::Bool=true)`  ... tests on module Einstein.
    """
    function testModule_Einstein(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-Einstein-new.sum")
        printstyled("\n\nTest the module  Einstein  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(36.); properties=[JAC.EinsteinX],
                                configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                                einsteinSettings=Einstein.Settings([E1, M1, E2, M2], true,
                                true, Tuple{Int64,Int64}[(5,0), (7,0), (10,0), (11,0), (12,0), (13,0), (14,0)], 0., 0., 10000. ) )

        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-Einstein-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-Einstein-new.sum"), "Einstein coefficients, t", 100)
        testPrint("testModule_Einstein()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_FormFactor(; short::Bool=true)`  ... tests on module FormFactor.
    """
    function testModule_FormFactor(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-FormFactor-new.sum")
        printstyled("\n\nTest the module  FormFactor  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=[JAC.FormF],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                formSettings=FormFactor.Settings([0.1], true, false, Int64[]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-FormFactor-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-FormFactor-new.sum"), "Standard and modifi", 6)
        testPrint("testModule_FormFactor()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_GreenFunction(; short::Bool=true)`  ... tests on module GreenFunction.
    """
    function testModule_GreenFunction(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-GreenFunction-new.sum")
        printstyled("\n\nTest the module  GreenFunction  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=[JAC.Green],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                greenSettings=GreenFunction.Settings(true, true, false, Int64[]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-GreenFunction-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-GreenFunction-new.sum"), "Green's function summary:", 6)
        testPrint("testModule_GreenFunction()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_Hfs(; short::Bool=true)`  ... tests on module Hfs.
    """
    function testModule_Hfs(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-Hfs-b-new.sum")
        printstyled("\n\nTest the module  Hfs  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0); properties=[JAC.HFS],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                hfsSettings=Hfs.Settings(true, true, true, true, true, true, false, Int64[] ) )

        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-Hfs-b-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-Hfs-b-new.sum"), "Level        J^P           Energy", 20)
        testPrint("testModule_Hfs()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_IsotopeShift(; short::Bool=true)`  ... tests on module IsotopeShift.
    """
    function testModule_IsotopeShift(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-IsotopeShift-new.sum")
        printstyled("\n\nTest the module  IsotopeShift  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=[JAC.Isotope],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                isotopeSettings=IsotopeShift.Settings(true, false, false, true, false, Int64[], "method-1") )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-IsotopeShift-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-IsotopeShift-new.sum"), "IsotopeShift parameters and amplitudes:", 15)
        testPrint("testModule_IsotopeShift()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_LandeZeeman(; short::Bool=true)`  ... tests on module LandeZeeman.
    """
    function testModule_LandeZeeman(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-LandeZeeman-new.sum")
        printstyled("\n\nTest the module  LandeZeeman  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.75, AngularJ64(5//2), 1.0, 2.0);
                                properties=[JAC.LandeJ],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                zeemanSettings=LandeZeeman.Settings(true, true, true, true, 0., true, false, Int64[] ) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-LandeZeeman-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-LandeZeeman-new.sum"), "Lande g_J factors and Zeeman amplitudes:", 30)
        testPrint("testModule_LandeZeeman()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_MultiPhotonDeExcitation(; short::Bool=true)`  ... tests on module MultiPhotonDeExcitation.
    """
    function testModule_MultiPhotonDeExcitation(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-MultiPhotonDeExcitation-new.sum")
        printstyled("\n\nTest the module  MultiPhotonDeExcitation  ... \n", color=:cyan)
        ### Make the tests
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        println(iostream, "Make the comparison with approved data for ... test-MultiPhotonDeExcitation-new.sum")
        success = true
        ## success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-MultiPhotonDeExcitation-approved.sum"),
        ##                             joinpath(@__DIR__, "..", "test", "test-MultiPhotonDeExcitation-new.sum"), "xxx", 100)
        testPrint("testModule_MultiPhotonDeExcitation()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_MultipoleMoment(; short::Bool=true)`  ... tests on module MultipoleMoment.
    """
    function testModule_MultipoleMoment(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-MultipoleMoment-new.sum")
        ### Make the tests
        printstyled("\n\nTest the module  MultipoleMoment  ... \n", color=:cyan)
        grid = Radial.Grid("grid: exponential")
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); grid=grid, properties=Basics.AtomicLevelProperty[],
                                configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")] )

        wxa  = perform(wa; output=true)
        wma  = wxa["multiplet:"]

        flow = 6;    fup = 8;   ilow = 1;   iup = 3
        println("\n\nDipole amplitudes:\n")
        for  finalLevel in wma.levels
            for  initialLevel in wma.levels
                if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
                    MultipoleMoment.dipoleAmplitude(finalLevel, initialLevel, grid; display=true)
                end
            end
        end
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-MultipoleMoment-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-MultipoleMoment-new.sum"), "Dipole amplitude", 2)
        testPrint("testModule_MultipoleMoment()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_MultipolePolarizibility(; short::Bool=true)`  ... tests on module MultipolePolarizibility.
    """
    function testModule_MultipolePolarizibility(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-MultipolePolarizibility-new.sum")
        printstyled("\n\nTest the module  MultipolePolarizibility  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=[JAC.Polarity],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                polaritySettings=MultipolePolarizibility.Settings(EmMultipole[], 0, 0, Float64[], false, false, Int64[]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-MultipolePolarizibility-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-MultipolePolarizibility-new.sum"),
                                "Multipole polarizibilities and amplitudes:", 5)
        testPrint("testModule_MultipolePolarizibility()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_ParityNonConservation(; short::Bool=true)`  ... tests on module ParityNonconservation.
    """
    function testModule_ParityNonConservation(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-ParityNonConservation-new.sum")
        printstyled("\n\nTest the module  ParityNonConservation  ... \n", color=:cyan)
        ### Make the tests
        grid = Defaults.getDefaults("standard grid")
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=Basics.AtomicLevelProperty[],
                                configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")] )

        wxa  = perform(wa; output=true)
        wma  = wxa["multiplet:"]

        nModel = Nuclear.Model(26.);    flow = 6;    fup = 8;   ilow = 1;   iup = 3

        println("\n\nWeak-charge and Schiff-moment amplitudes:\n")
        for  finalLevel in wma.levels
            for  initialLevel in wma.levels
                if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
                    ParityNonConservation.weakChargeAmplitude(finalLevel, initialLevel, nModel, grid; display=true)
                    ParityNonConservation.schiffMomentAmplitude(finalLevel, initialLevel, nModel, grid; display=true)
                end
            end
        end
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-ParityNonConservation-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-ParityNonConservation-new.sum"), "weak-charge amplitude", 12)
        testPrint("testModule_ParityNonConservation()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PlasmaShift(; short::Bool=true)`  ... tests on module PlasmaShift.
    """
    function testModule_PlasmaShift(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PlasmaShift-new.sum")
        printstyled("\n\nTest the module  PlasmaShift  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); properties=[JAC.Plasma],
                                configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")],
                                plasmaSettings=PlasmaShift.Settings() )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PlasmaShift-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PlasmaShift-new.sum"), "Plasma screening", 3)
        testPrint("testModule_PlasmaShift()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoExcitation(; short::Bool=true)`  ... tests on module PhotoExcitation.
    """
    function testModule_PhotoExcitation(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoExcitation-new.sum")
        printstyled("\n\nTest the module  PhotoExcitation  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                                initialConfigs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                                finalConfigs  =[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                                process = JAC.PhotoExc,
                                processSettings=PhotoExcitation.Settings([E1, M1], [UseCoulomb, UseBabushkin], true, true, true, true,
                                                                    false, Tuple{Int64,Int64}[], 0., 0., 1.0e6, Basics.ExpStokes(0., 0., 0.) ) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoExcitation-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoExcitation-new.sum"),
                                "Photo-excitation cross sections for (completely)", 200)
        testPrint("testModule_PhotoExcitation()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoExcitationAutoion(; short::Bool=true)`  ... tests on module PhotoExcitationAutoion.
    """
    function testModule_PhotoExcitationAutoion(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoExcitationAutoion-new.sum")
        printstyled("\n\nTest the module  PhotoExcitationAutoion  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.); grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 1.5e-2, NoPoints = 600),
                                initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                                intermediateConfigs=[Configuration("1s 2s^2"), Configuration("1s 2p^2")],
                                finalConfigs=[Configuration("1s^2")],
                                process = JAC.PhotoExcAuto,
                                processSettings=PhotoExcitationAutoion.Settings([E1, M1], [UseCoulomb, UseBabushkin], true, true, Tuple{Int64,Int64,Int64}[(1,1,1)], 2)  )

    wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        println(iostream, "Make the comparison with approved data for ... test-PhotoExcitationAutoion-new.sum")
        success = true
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoExcitationAutoion-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoExcitationAutoion-new.sum"), "xxx", 10)
        testPrint("testModule_PhotoExcitationAutoion()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoExcitationFluores(; short::Bool=true)`  ... tests on module PhotoExcitationFluores.
    """
    function testModule_PhotoExcitationFluores(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoExcitationFluores-new.sum")
        printstyled("\n\nTest the module  PhotoExcitationFluores  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(26.),
                                initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                                intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                                finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                                process = JAC.PhotoExcFluor,
                                processSettings=PhotoExcitationFluores.Settings([E1, M1], [UseCoulomb, UseBabushkin],
                                                                                true, true, Tuple{Int64,Int64,Int64}[(1,1,1)])  )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        println(iostream, "Make the comparison with approved data for ... test-PhotoExcitationFluores-new.sum")
        success = true
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoExcitationFluores-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoExcitationFluores-new.sum"),  "xxx", 10)
        testPrint("testModule_PhotoExcitationFluores()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoIonization(; short::Bool=true)`  ... tests on module PhotoIonization.
    """
    function testModule_PhotoIonization(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoIonization-new.sum")
        printstyled("\n\nTest the module  PhotoIonization  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(36.); grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 800),
                                initialConfigs=[Configuration("1s^2 2s^2 2p^6")],
                                finalConfigs  =[Configuration("1s^2 2s^2 2p^5"), Configuration("1s^2 2s 2p^6") ],
                                process = JAC.Photo,
                                processSettings=PhotoIonization.Settings([E1, M1], [UseCoulomb, UseBabushkin], [3000., 4000.], false, true, true, true,
                                                true, Tuple{Int64,Int64}[(1,1), (1,2)], ExpStokes(1., 0., 0.)) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoIonization-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoIonization-new.sum"), "Total photoionization c", 8)
        testPrint("testModule_PhotoIonization()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoRecombination(; short::Bool=true)`  ... tests on module PhotoRecombination.
    """
    function testModule_PhotoRecombination(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoRecombination-new.sum")
        printstyled("\n\nTest the module  PhotoRecombination  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(36.); grid=Radial.Grid("grid: by given parameters"; rnt = 2.0e-6, h = 5.0e-2, hp = 2.0e-2, NoPoints = 600),
                                initialConfigs=[Configuration("1s^2 2s^2 2p^5"), Configuration("1s^2 2s 2p^6") ],
                                finalConfigs  =[Configuration("1s^2 2s^2 2p^6")],
                                process = JAC.Rec,
                                processSettings=PhotoRecombination.Settings([E1, M1], [UseCoulomb, UseBabushkin], [10., 20.], [0.],
                                                false, true, true, true, true, Tuple{Int64,Int64}[(1,1)]) )
        wb = perform(wa)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoRecombination-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoRecombination-new.sum"), "Photorecombination cross sections for ", 12)
        testPrint("testModule_PhotoRecombination()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_PhotoEmission(; short::Bool=true)`  ... tests on module PhotoEmission.
    """
    function testModule_PhotoEmission(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-PhotoEmission-new.sum")
        printstyled("\n\nTest the module  PhotoEmission  ... \n", color=:cyan)
        ### Make the tests
        wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                                initialConfigs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                                finalConfigs  =[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                                process = JAC.Radiative,
                                processSettings=PhotoEmission.Settings([E1, M1, E2, M2], [UseCoulomb, UseBabushkin], true, true,
                                true, Tuple{Int64,Int64}[(5,0), (7,0), (10,0), (11,0), (12,0), (13,0), (14,0), (15,0), (16,0)], 0., 0., 10000. ) )
        ##x streamDummy = open(pwd() * "/runtests.dummy", "w")
        ##x redirect_stdout(streamDummy) do
        wb = perform(wa)
        ##x end
        ##x close(streamDummy)
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-PhotoEmission-approved.sum"),
                                    joinpath(@__DIR__, "..", "test", "test-PhotoEmission-new.sum"), "Einstein coefficients, t", 100)
        testPrint("testModule_PhotoEmission()::", success)
        return(success)
    end

    """
    `TestFrames.testModule_RayleighCompton(; short::Bool=true)`  ... tests on module RayleighCompton.
    """
    function testModule_RayleighCompton(; short::Bool=true)
        Defaults.setDefaults("print summary: open", "test-RayleighCompton-new.sum")
        printstyled("\n\nTest the module  RayleighCompton  ... \n", color=:cyan)
        ### Make the tests
        ###
        Defaults.setDefaults("print summary: close", "")
        # Make the comparison with approved data
        printTest, iostream = Defaults.getDefaults("test flag/stream")
        println(iostream, "Make the comparison with approved data for ... test-RayleighCompton-new.sum")
        success = true
        ## success = testCompareFiles( joinpath(@__DIR__, "..", "test", "approved", "test-RayleighCompton-approved.sum"),
        ##                             joinpath(@__DIR__, "..", "test", "test-RayleighCompton-new.sum"), "xxx", 100)
        testPrint("testModule_RayleighCompton()::", success)
        return(success)
    end

end # module


text:
using PharmaceuticalClassification
using Test

import Serialization

include("test_utils/generate-synthetic-data.jl")
include("test_utils/utils.jl")

@testset "PharmaceuticalClassification.jl" begin
    include("unit.jl")
    include("integration.jl")
end