text:
struct StackedView{T<:Number,N,A<:Tuple{Vararg{AbstractArray{T}}}} <: AbstractArray{T,N}
    parents::A

    function StackedView{T,N,A}(parents::A) where {T<:Number,N,A<:Tuple{Vararg{AbstractArray{T}}}}
        inds = axes(parents[1])
        length(inds) == N-1 || throw(DimensionMismatch("component arrays must be of dimension \$(N-1), got \$(length(inds))"))
        for i = 2:length(parents)
            axes(parents[i]) == inds || throw(DimensionMismatch("all arrays must have the same indices, got \$inds and \$(axes(parents[i]))"))
        end
        new(parents)
    end
end

"""
    StackedView(B, C, ...) -> A

Present arrays `B`, `C`, etc, as if they are separate channels along
the first dimension of `A`. In particular,

    B == A[1,:,:...]
    C == A[2,:,:...]

and so on. Combined with `colorview`, this allows one to combine two
or more grayscale images into a single color image.

See also: [`colorview`](@ref).
"""
StackedView

@inline Base.size(V::StackedView) = (length(V.parents), size(V.parents[1])...)
@inline Base.axes(V::StackedView) = (Base.OneTo(length(V.parents)), axes(V.parents[1])...)

@inline function Base.getindex(V::StackedView{T,N}, I::Vararg{Int,N}) where {T,N}
    i1, itail = I[1], tail(I)
    P = V.parents
    @boundscheck (1 <= i1) & (i1 <= length(P)) & checkbounds(Bool, P[1], itail...) ||
        Base.throw_boundserror(V, I)
    _unsafe_getindex(i1, itail, P...)
end

@inline function Base.setindex!(V::StackedView{T,N}, val, I::Vararg{Int,N}) where {T,N}
    i1, itail = I[1], tail(I)
    P = V.parents
    @boundscheck (1 <= i1) & (i1 <= length(P)) & checkbounds(Bool, P[1], itail...) ||
        Base.throw_boundserror(V, I)
    _unsafe_setindex!(i1, itail, convert(T, val), P...)
    val
end

# This does the same thing that V.parents[i1][itail...] would do,
# except in a type-stable way. The cost is the introduction of
# branches.
@inline function _unsafe_getindex(idx, I, A::AbstractArray, As...)
    if idx == 1
        @inbounds ret = A[I...]
        return ret
    end
    _unsafe_getindex(idx-1, I, As...)
end
_unsafe_getindex(idx, I) = error("ran out of arrays; this shouldn't happen")

# For getting all of the channels (e.g., for ColorView)
@inline function _unsafe_getindex_all(I, A, As...)
    @inbounds ret = A[I...]
    (ret, _unsafe_getindex_all(I, As...)...)
end
_unsafe_getindex_all(I) = ()

@inline function _unsafe_setindex!(idx, I, val, A::AbstractArray, As...)
    if idx == 1
        @inbounds A[I...] = val
        return val
    end
    _unsafe_setindex!(idx-1, I, val, As...)
end
_unsafe_setindex!(idx, I, val) = error("ran out of arrays; this shouldn't happen")

# For setting all of the channels (e.g., for ColorView), one `val` per parent
@inline function _unsafe_setindex_all!(I, vals::NTuple{N}, As::NTuple{N,AbstractArray}) where N
    val1, valrest = vals[1], tail(vals)
    A1, Arest = As[1], tail(As)
    @inbounds A1[I...] = val1
    _unsafe_setindex_all!(I, valrest, Arest)
end
_unsafe_setindex_all!(I, ::Tuple{}, ::Tuple{}) = nothing

# Performance optimizations
@inline getchannels(P::StackedView, ::Type{C}, I) where {C<:Color2} = _unsafe_getindex_all(I, P.parents...)
@inline getchannels(P::StackedView, ::Type{C}, I) where {C<:Color3} = _unsafe_getindex_all(I, P.parents...)
@inline getchannels(P::StackedView, ::Type{C}, I) where {C<:Color4} = _unsafe_getindex_all(I, P.parents...)

@inline setchannels!(P::StackedView, val::Color2, I) = _unsafe_setindex_all!(I, (comp1(val),alpha(val)), P.parents)
@inline setchannels!(P::StackedView, val::Color3, I) = _unsafe_setindex_all!(I, (comp1(val),comp2(val),comp3(val)), P.parents)
@inline setchannels!(P::StackedView, val::Color4, I) = _unsafe_setindex_all!(I, (comp1(val),comp2(val),comp3(val),alpha(val)), P.parents)

# When overlaying 2 images, you often might want one color channel to be black
struct ZeroArrayPromise{T} end
const zeroarray = ZeroArrayPromise{Union{}}()

struct ZeroArray{T,N,R<:AbstractUnitRange} <: AbstractArray{T,N}
    inds::NTuple{N,R}
end

ZeroArrayPromise{T}(inds::NTuple{N,R}) where {T,N,R<:AbstractUnitRange} = ZeroArray{T,N,R}(inds)
Base.eltype(::Type{ZeroArrayPromise{T}}) where {T} = T

Base.axes(A::ZeroArray) = A.inds
Base.size(A::ZeroArray) = length.(A.inds)
Base.getindex(A::ZeroArray{T,N}, I::Vararg{Int,N}) where {T,N} = zero(T)

@inline function StackedView(arrays::Union{AbstractArray,ZeroArrayPromise}...)
    T = promote_eleltype_all(arrays...)
    stackedview(T, arrays...)
end
@inline function StackedView{T}(arrays::Union{AbstractArray,ZeroArrayPromise}...) where T<:Number
    stackedview(T, arrays...)
end

# Now, it seems we should be able to do this with uppercase-typed
# calls, but there seems to be some inference bug, and using a
# function with a different name works around it.
@inline function stackedview(::Type{T}, arrays::Union{AbstractArray,ZeroArrayPromise}...) where T<:Number
    inds = firstinds(arrays...)
    arrays_take = take_zeros(T, inds, arrays...)
    arrays_T = map(A->of_eltype(zero(T), A), arrays_take)
    # To compute N = length(inds)+1 in a type-stable fashion, we have
    # to use tuple tricks (i.e., make a tuple of length(inds)+1)
    _stackedview(T, (length(arrays), inds...), arrays_T)
end
_stackedview(::Type{T}, ::Tuple{Vararg{Any,N}}, arrays) where {T,N} = StackedView{T,N,typeof(arrays)}(arrays)

@inline firstinds(A::AbstractArray, Bs...) = axes(A)
@inline firstinds(::ZeroArrayPromise, Bs...) = firstinds(Bs...)
firstinds() = error("not all arrays can be zeroarray")

@inline take_zeros(T, inds, ::ZeroArrayPromise, Bs...) = (ZeroArrayPromise{T}(inds), take_zeros(T, inds, Bs...)...)
@inline take_zeros(T, inds, A::AbstractArray, Bs...) = (A, take_zeros(T, inds, Bs...)...)
take_zeros(T, inds) = ()

# Extensions of PaddedViews
function PaddedViews.paddedviews(fillvalue, As::Union{AbstractArray,ZeroArrayPromise}...)
    inds = PaddedViews.outerinds(As...)
    map(A->PaddedView(fillvalue, A, inds), As)
end

PaddedViews.PaddedView(fillvalue, zap::ZeroArrayPromise, indices) = zap

@inline PaddedViews.outerinds(A::ZeroArrayPromise, Bs...) = PaddedViews.outerinds(Bs...)
@inline PaddedViews._outerinds(inds, A::ZeroArrayPromise, Bs...) =
    PaddedViews._outerinds(inds, Bs...)


text:
######################################################################
# An example of a Combined chart in XlsxWriter.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_combined.xlsx")
	ws = add_worksheet!(wb)

	# Add a format for the headings.
	bold = add_format!(wb, Dict("bold"  => 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Number" "Batch 1" "Batch 2"]
	data = [
		[2 3 4 5 6 7];
		[10 40 50 20 10 50];
		[30 60 70 50 40 30];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])
	write_column!(ws, "C2", data[3,:])

	#
	# In the first example we will create a combined column and line chart.
	# They will share the same X and Y axes.
	#

	# Create a new column  This will use this as the primary chart.
	column_chart1 = add_chart!(wb, Dict("type" => "column"))

	# Configure the data series for the primary chart.
	add_series!(column_chart1, Dict(
		"name" => "=Sheet1!\$B\$1",
		"categories" => "=Sheet1!\$A\$2:\$A\$7",
		"values" => "=Sheet1!\$B\$2:\$B\$7",
	))

	# Create a new column  This will use this as the secondary chart.
	line_chart1 = add_chart!(wb, Dict("type" => "line"))

	# Configure the data series for the secondary chart.
	add_series!(line_chart1, Dict(
		"name" => "=Sheet1!\$C\$1",
		"categories" => "=Sheet1!\$A\$2:\$A\$7",
		"values" => "=Sheet1!\$C\$2:\$C\$7",
	))

	# Combine the
	combine!(column_chart1, line_chart1)

	# Add a chart title and some axis  Note, this is done via the
	# primary chart.
	set_title!(column_chart1, Dict( "name" => "Combined chart - same Y axis"))
	set_x_axis!(column_chart1, Dict("name" => "Test number"))
	set_y_axis!(column_chart1, Dict("name" => "Sample length (mm)"))

	# Insert the chart into the worksheet
	insert_chart!(ws, "E2", column_chart1)

	#
	# In the second example we will create a similar combined column and line
	# chart except that the secondary chart will have a secondary Y axis.
	#

	# Create a new column  This will use this as the primary chart.
	column_chart2 = add_chart!(wb, Dict("type" => "column"))

	# Configure the data series for the primary chart.
	add_series!(column_chart2, Dict(
		"name" => "=Sheet1!\$B\$1",
		"categories" => "=Sheet1!\$A\$2:\$A\$7",
		"values" => "=Sheet1!\$B\$2:\$B\$7",
	))

	# Create a new column  This will use this as the secondary chart.
	line_chart2 = add_chart!(wb, Dict("type" => "line"))

	# Configure the data series for the secondary  We also set a
	# secondary Y axis via (y2_axis). This is the only difference between
	# this and the first example, apart from the axis label below.
	add_series!(line_chart2, Dict(
		"name" => "=Sheet1!\$C\$1",
		"categories" => "=Sheet1!\$A\$2:\$A\$7",
		"values" => "=Sheet1!\$C\$2:\$C\$7",
		"y2_axis" => true,
	))

	# Combine the charts.
	combine!(column_chart2, line_chart2)

	# Add a chart title and some axis labels.
	set_title!(column_chart2, Dict(  "name" => "Combine chart - secondary Y axis"))
	set_x_axis!(column_chart2, Dict( "name" => "Test number"))
	set_y_axis!(column_chart2, Dict( "name" => "Sample length (mm)"))

	# Note: the y2 properties are on the secondary chart.
	set_y2_axis!(line_chart2, Dict("name" => "Target length (mm)"))

	# Insert the chart into the worksheet
	insert_chart!(ws, "E18", column_chart2)

	close(wb)
	isfile("chart_combined.xlsx")
end

test()

ext: jl


text:
function discretizeperiphery!(system::CVSystem)
    for i in 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            system.branches.term[i].P = zeros(system.solverparams.numsteps+1,5);
            system.branches.term[i].V = zeros(system.solverparams.numsteps+1,5);
            system.branches.term[i].Q = zeros(system.solverparams.numsteps+1,5);
        end
    end

    append!(system.svc.P,zeros(system.solverparams.numsteps+1))
    append!(system.svc.V,zeros(system.solverparams.numsteps+1))
    append!(system.svc.Q,zeros(system.solverparams.numsteps+1))
    append!(system.ivc.P,zeros(system.solverparams.numsteps+1))
    append!(system.ivc.V,zeros(system.solverparams.numsteps+1))
    append!(system.ivc.Q,zeros(system.solverparams.numsteps+1))
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# varinteran.jl - operaciones básicas para computar variaciones interanuales de
# ìndices de precios

"""
    varinteran(idx::AbstractVector, base_index::Real = 100) -> Vector{<:AbstractFloat}
    varinteran(cpimat::AbstractMatrix, base_index::Real = 100) -> Matrix{<:AbstractFloat}
    varinteran(cpimat::AbstractMatrix, base_index::AbstractVector) -> Matrix{<:AbstractFloat}

Obtiene variaciones interanuales del vector `idx` o de la matriz `cpimat`
utilizando como índice base el número o vector `base_index`.

- Si `base_index` es un vector, se obtienen las variaciones interanuales
  utilizando diferentes índices base para cada columna de `cpimat`. El vector
  `base_index` debe tener la misma cantidad de columnas que `cpimat`.
"""
function varinteran end

"""
    varinteran!(v::AbstractVector, idx::AbstractVector, base_index::Real = 100) -> Vector{<:AbstractFloat}

Computa las variaciones interanuales del vector `idx` utilizando como índice
base `base_index`. Si se provee `v`, los resultados son guardados en este
vector, en vez del mismo `idx`.

- El vector `v` tiene 11 observaciones menos que `idx`.
"""
function varinteran!(v::AbstractVector, idx::AbstractVector, base_index::Real = 100)
    l = length(v)
    for i in l:-1:2
        @inbounds v[i] = 100 * (idx[i+11] / idx[i-12+11] - 1)
    end
    v[1] = 100 * (idx[12] / base_index - 1)
end

# Implementación sobre AbstractVector
function varinteran(idx::AbstractVector, base_index::Real = 100)
    r = length(idx)
    v = zeros(eltype(idx), r-11)
    varinteran!(v, idx, base_index)
    v
end

# Implementación sobre AbstractMatrix
function varinteran(cpimat::AbstractMatrix, base_index::Real = 100)
    r, c = size(cpimat)
    vmat = zeros(eltype(cpimat), r-11, c)
    for j in 1:c
        vcol = @view vmat[:, j]
        idxcol = @view cpimat[:, j]
        varinteran!(vcol, idxcol, base_index)
    end
    vmat
end

# Implementación que utiliza diferentes índices base
function varinteran(cpimat::AbstractMatrix, base_index::AbstractVector)
    r, c = size(cpimat)
    vmat = zeros(eltype(cpimat), r-11, c)
    for j in 1:c
        vcol = @view vmat[:, j]
        idxcol = @view cpimat[:, j]
        varinteran!(vcol, idxcol, base_index[j])
    end
    vmat
end


text:
module test_lupin_enums

using Test
using Lupin

# @enum DepthTest Never Less Equal LessEqual Greater NotEqual GreaterEqual Always
#
# @test DepthTest.Never == Never

end # module test_lupin_enums


text:
module RecursiveSparseBlocks

export SparseMatrixRSB

include("../deps/deps.jl")
include("constants.jl")

const RSBCooIdx = Cint
const RSBNnzIdx = Cint
const RSBType   = Cchar
const RSBFlags  = Cint
const RSBErr    = Cint
const RSBChar   = Cchar
const RSBBlkIdx = Cint
const RSBTrans  = RSBFlags

function __init__()
    # init library with default options
    err = ccall((:rsb_lib_init, librsb), RSBErr,
                (Ptr{Void},), C_NULL)

    if err != 0
        error("rsb_lib_init failed with code $(err)")
    end
end

type SparseMatrixRSB{T} <: AbstractSparseMatrix{T, Cint}
    ptr::Ptr{Void}
end

function rsb_mtx_free(A::SparseMatrixRSB)
    ccall((:rsb_mtx_free, librsb), Ptr{Void}, (Ptr{Void},), A.ptr)
end

function check_rsb_error(err::RSBErr)
    if err != RSB_ERR_NO_ERROR
        buflen = 5000
        buf = Array{UInt8}(buflen)
        ccall((:rsb_strerror_r, librsb), RSBErr, (RSBErr, Ptr{RSBChar}, Csize_t),
              err, buf, buflen)
        i = findfirst(0x0, buf)
        if i == 0
            error("SparseMatrixRSB operation failed with code $(err)")
        else
            error(String(buf[1:i-1]))
        end
    end
end

@inline function rsb_type_code{T}(::Type{T})
    typ = T == Float32    ? 'S' :
          T == Float64    ? 'D' :
          T == Complex64  ? 'C' :
          T == Complex128 ? 'Z' :
          error("librsb does not support matrices of type $T")
    return typ
end

function SparseMatrixRSB{T}(A::SparseMatrixCSC{T})
    typ = rsb_type_code(T)
    m, n = size(A)
    nnz = length(A.nzval)

    IA = Vector{RSBCooIdx}(A.rowval)
    CP = Vector{RSBCooIdx}(A.colptr)

    # make zero based
    IA .-= 1
    CP .-= 1

    err = Ref{RSBErr}()
    ptr = ccall((:rsb_mtx_alloc_from_csc_const, librsb), Ptr{Void},
                (Ptr{Void}, Ptr{RSBCooIdx}, Ptr{RSBCooIdx}, RSBNnzIdx, RSBType,
                 RSBCooIdx, RSBCooIdx, RSBBlkIdx, RSBBlkIdx, RSBFlags,
                 Ptr{RSBErr}),
                A.nzval, IA, CP, nnz, typ, m, n,
                RSB_DEFAULT_ROW_BLOCKING, RSB_DEFAULT_COL_BLOCKING,
                RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS, err)
    check_rsb_error(err.x)

    ret = SparseMatrixRSB{T}(ptr)
    finalizer(ret, rsb_mtx_free)
    return ret
end

function SparseMatrixRSB{T}(I_, J_, V::Vector{T}, m::Integer, n::Integer)
    typ = rsb_type_code(T)
    nnz = length(V)
    I = Vector{RSBCooIdx}(I_)
    J = Vector{RSBCooIdx}(J_)

    # make zero based
    I .-= 1
    J .-= 1

    err = Ref{RSBErr}()
    ptr = ccall((:rsb_mtx_alloc_from_coo_const, librsb), Ptr{Void},
                (Ptr{Void}, Ptr{RSBCooIdx}, Ptr{RSBCooIdx}, RSBNnzIdx, RSBType,
                 RSBCooIdx, RSBCooIdx, RSBBlkIdx, RSBBlkIdx, RSBFlags,
                 Ptr{RSBErr}),
                V, I, J, nnz, typ, m, n,
                RSB_DEFAULT_ROW_BLOCKING, RSB_DEFAULT_COL_BLOCKING,
                RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS, err)
    check_rsb_error(err.x)

    ret = SparseMatrixRSB{T}(ptr)
    finalizer(ret, rsb_mtx_free)
    return ret
end

function Base.convert{T}(::Type{SparseMatrixCSC}, A::SparseMatrixRSB{T})
    I, J, V = findnz(A)
    m, n = size(A)
    return sparse(I, J, V, m, n)
end

function Base.findnz{T}(A::SparseMatrixRSB{T})
    N = nnz(A)
    V = Array{T}(N)
    I = Array{RSBCooIdx}(N)
    J = Array{RSBCooIdx}(N)

    err = ccall((:rsb_mtx_get_coo, librsb), RSBErr,
                (Ptr{Void}, Ptr{Void}, Ptr{RSBCooIdx}, Ptr{RSBCooIdx},
                 RSBFlags),
                A.ptr, V, I, J, RSB_FLAG_FORTRAN_INDICES_INTERFACE)
    check_rsb_error(err)

    return (I, J, V)
end

function Base.copy{T}(A::SparseMatrixRSB{T})
    typ = rsb_type_code(T)
    Bptr = Ref{Ptr{Void}}()
    err = ccall((:rsb_mtx_clone, librsb), RSBErr,
                (Ptr{Ptr{Void}}, RSBType, RSBTrans, Ptr{Void}, Ptr{Void}, RSBFlags),
                BPtr, typ, RSB_TRANSPOSITION_N, C_NULL, A.ptr,
                RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS)
    check_rsb_error(err)

    ret = SparseMatrixRSB{T}(Bptr.x)
    finalizer(ret, rsb_mtx_free)
    return ret
end

function Base.eltype{T}(A::SparseMatrixRSB{T})
    return T
end

function Base.show(io::IO, ::MIME"text/plain", A::SparseMatrixRSB)
    m, n = size(A)
    print(io, m, "×", n, " sparse RSB matrix with ", nnz(A), " ", eltype(A),
          " stored entries")
    if nnz(A) != 0
        print(io, ":")
        show(io, A)
    end
end

Base.show(io::IO, S::SparseMatrixRSB) = Base.show(convert(IOContext, io),
                                                  S::SparseMatrixRSB)

function Base.show(io::IOContext, S::SparseMatrixRSB)
    # Taken mostly from SparseMatrixCSC show function

    if nnz(S) == 0
        return show(io, MIME("text/plain"), S)
    end

    limit::Bool = get(io, :limit, false)
    if limit
        rows = displaysize(io)[1]
        half_screen_rows = div(rows - 8, 2)
    else
        half_screen_rows = typemax(Int)
    end
    pad = ndigits(max(size(S)...))
    k = 0
    sep = "\n  "
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end

    I, J, V = findnz(S)
    k = 0
    for (i, j, v) in zip(I, J, V)
        if k < half_screen_rows || k > nnz(S)-half_screen_rows
            print(io, sep, '[',
                  rpad(i, pad), ", ",
                  lpad(j, pad), "]  =  ")
            Base.show(io, v)
        elseif k == half_screen_rows
            print(io, sep, '\u22ee')
        end
        k += 1
    end
end

function Base.nnz(A::SparseMatrixRSB)
    out = Ref{RSBNnzIdx}()
    err = ccall((:rsb_mtx_get_info, librsb), RSBErr,
                (Ptr{Void}, Cint, Ptr{Void}),
                A.ptr, RSB_MIF_MATRIX_NNZ__TO__RSB_NNZ_INDEX_T, out)
    check_rsb_error(err)
    return Int(out.x)
end

function rsb_get_nrows(A::SparseMatrixRSB)
    out = Ref{RSBCooIdx}()
    err = ccall((:rsb_mtx_get_info, librsb), RSBErr,
                (Ptr{Void}, Cint, Ptr{Void}),
                A.ptr, RSB_MIF_MATRIX_ROWS__TO__RSB_COO_INDEX_T, out)
    check_rsb_error(err)
    return out.x
end

function rsb_get_ncols(A::SparseMatrixRSB)
    out = Ref{RSBCooIdx}()
    err = ccall((:rsb_mtx_get_info, librsb), RSBErr,
                (Ptr{Void}, Cint, Ptr{Void}),
                A.ptr, RSB_MIF_MATRIX_COLS__TO__RSB_COO_INDEX_T, out)
    check_rsb_error(err)
    return out.x
end

function Base.size(A::SparseMatrixRSB)
    return (Int(rsb_get_nrows(A)),
            Int(rsb_get_ncols(A)))
end

function Base.size(A::SparseMatrixRSB, k::Integer)
    if k == 1
        return Int(rsb_get_nrows(A))
    elseif k == 2
        return Int(rsb_get_ncols(A))
    else
        return 1
    end
end

function Base.getindex{T}(A::SparseMatrixRSB{T}, i::Integer, j::Integer)
    yptr = Ref{T}()
    iptr = Ref{RSBCooIdx}(i)
    jptr = Ref{RSBCooIdx}(j)
    err = ccall((:rsb_mtx_get_vals, librsb), RSBErr,
                (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, RSBNnzIdx,
                 RSBFlags),
                A.ptr, yptr, iptr, jptr, 1, RSB_FLAG_FORTRAN_INDICES_INTERFACE)
    if err == -1
        return zero(T)
    end
    check_rsb_error(err)
    return yptr.x
end

function Base.setindex!{T}(A::SparseMatrixRSB{T}, value, i::Integer, j::Integer)
    yptr = Ref{T}(value)
    iptr = Ref{RSBCooIdx}(i)
    jptr = Ref{RSBCooIdx}(j)
    err = ccall((:rsb_mtx_set_vals, librsb), RSBErr,
                (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, RSBNnzIdx,
                 RSBFlags),
                A.ptr, yptr, iptr, jptr, 1, RSB_FLAG_FORTRAN_INDICES_INTERFACE)
    check_rsb_error(err)
    return yptr.x
end

function Base.diag{T}(A::SparseMatrixRSB{T})
    v = Vector{T}(min(size(A)...))
    err = ccall((:rsb_mtx_get_vec, librsb), RSBErr,
                (Ptr{Void}, Ptr{Void}, Cint),
                A.ptr, v, RSB_EXTF_DIAG)
    check_rsb_error(err)
    return v
end

function Base.full{T}(A::SparseMatrixRSB{T})
    B = zeros(T, size(A))
    I, J, V = findnz(A)
    for (i, j, v) in zip(I, J, V)
        B[i, j] = v
    end
    return B
end

# Sparse matrix by vector multiplication
for (f, transA) in ((:A_mul_B!,  RSB_TRANSPOSITION_N),
                    (:At_mul_B!, RSB_TRANSPOSITION_T),
                    (:Ac_mul_B!, RSB_TRANSPOSITION_C))
    @eval begin
        function Base.$f{T}(y::Vector{T}, A::SparseMatrixRSB{T}, x::Vector{T})
            m, n = size(A)
            if $transA == RSB_TRANSPOSITION_N
                @assert length(x) == n
                @assert length(y) == m
            else
                @assert length(x) == m
                @assert length(y) == n
            end

            alpha = Ref{T}(one(T))
            beta  = Ref{T}(zero(T))
            err = ccall((:rsb_spmv, librsb), RSBErr,
                        (RSBTrans, Ptr{Void}, Ptr{Void}, Ptr{Void}, RSBCooIdx,
                         Ptr{Void}, Ptr{Void}, RSBCooIdx),
                        $transA, alpha, A.ptr, x, 1, beta, y, 1)
            check_rsb_error(err)
            return y
        end
    end
end

# Sparse matrix by dense matrix multiplication
for (f, transA) in ((:A_mul_B!,  RSB_TRANSPOSITION_N),
                    (:At_mul_B!, RSB_TRANSPOSITION_T),
                    (:Ac_mul_B!, RSB_TRANSPOSITION_C))
    @eval begin
        function Base.$f{T}(C::Matrix{T}, A::SparseMatrixRSB{T}, B::Matrix{T})
            ma, na = size(A)
            mb, nb = size(B)
            mc, nc = size(C)

            if $transA == RSB_TRANSPOSITION_N
                @assert na == mb
                @assert (mc, nc) == (ma, nb)
            else
                @assert ma == mb
                @assert (mc, nc) == (na, nb)
            end

            alpha = Ref{T}(one(T))
            beta  = Ref{T}(zero(T))
            err = ccall((:rsb_spmm, librsb), RSBErr,
                        (RSBTrans, Ptr{Void}, Ptr{Void}, RSBCooIdx, RSBFlags,
                         Ptr{Void}, RSBNnzIdx, Ptr{Void}, Ptr{Void}, RSBNnzIdx),
                        $transA, alpha, A.ptr, nb,
                        RSB_FLAG_WANT_COLUMN_MAJOR_ORDER, B, mb, beta, C, na)
            check_rsb_error(err)
            return C
        end
    end
end

# TODO:
# get column and row sums
# convert{S, T}(::Type{SparseMatrixRSB{S}}, ::SparseMatrixRSB{T})
# transpose / transpose!
# conj / conj!
# A_ldiv_B! / trisolve (BLAS_?ussv, BLAS_?ussm)
# + (I think we just insert elements with the right flag set if we want to add)

end


text:
#=
Description:  Convert string to upper case
Author: Tushar Gautam
Date: 9th Jan 2022
=#

print("Enter the string: ")

str = readline()

# Strings are immutable, hence converting it to
# array of characters
str = Vector{Char}(str)
# This would also work: str = collect(str)

# Converting to uppercase
for (i, char) in enumerate(str)
    str[i] = uppercase(char)
end

# Converting array of chars to string
str = join(str)

println(str)


text:
@testset "Ambiguities" begin
    @test isempty(Base.Test.detect_ambiguities(TransposedVectors, Base, Core))
end

@testset "Core" begin
    v = [1,2,3]
    z = [1+im,2,3]

    @test TransposedVector(v) == [1 2 3]
    @test TransposedVector{Int}(v) == [1 2 3]
    @test size(TransposedVector{Int}(3)) === (1,3)
    @test size(TransposedVector{Int}(1,3)) === (1,3)
    @test size(TransposedVector{Int}((3,))) === (1,3)
    @test size(TransposedVector{Int}((1,3))) === (1,3)
    @test_throws Exception TransposedVector{Float64, Vector{Int}}(v)

    @test (v.')::TransposedVector == [1 2 3]
    @test (v')::TransposedVector == [1 2 3]
    @test (z.')::TransposedVector == [1+im 2 3]
    @test (z')::TransposedVector == [1-im 2 3]

    tv = v.'
    tz = z.'

    @test (tv.')::Vector == [1, 2, 3]
    @test (tv')::Vector == [1, 2, 3]
    @test (tz.')::Vector == [1+im, 2, 3]
    @test (tz')::Vector == [1-im, 2, 3]

    @test conj(tv) === tv
    @test conj(tz) == [1-im 2 3]

    @test isa(similar(tv), TransposedVector)
    @test isa(similar(tv, Float64), TransposedVector)
    @test isa(copy(tv), TransposedVector)

    @test tv[2] === v[2]
    @test tv[1,2] === v[2]

    @test (tv2 = copy(tv); tv2[2] = 6; tv2[2] === 6)
    @test (tv2 = copy(tv); tv2[1,2] = 6; tv2[2] === 6)

    @test length(tv) === 3
    @test size(tv) === (1,3)
    @test size(tv,1) === 1
    @test size(tv,2) === 3

    @test map(-, tv)::TransposedVector == [-1 -2 -3]
    @test (-).(tv)::TransposedVector == [-1 -2 -3]
    @test (-).(tv,1)::TransposedVector == [0  1  2]

    @test AbstractVector(tv)::Vector{Int} == [1, 2, 3]
    @test Vector{Float64}(tv)::Vector{Float64} == [1, 2, 3]
    @test TransposedVector{Float64,Vector{Float64}}(tv)::TransposedVector{Float64} == [1.0  2.0  3.0]

    y = rand(Complex{Float64},3)
    @test sumabs2(imag.(diag(y .+ y'))) < 1e-20
end


text:
# An (incredibly simplistic) example of how you might wish to handle logging in a TUI environment.

using Ansillary
using Logging

# Warning! Not threadsafe!
struct TUILogger <: AbstractLogger
    queue::Vector
    displayed::Vector
end

TUILogger() = TUILogger([], [])

Logging.shouldlog(::TUILogger, _, _, _, _) = true

Logging.min_enabled_level(::TUILogger) = Logging.Info

Logging.handle_message(
    logger::TUILogger,
    level,
    message,
    module_,
    group,
    id,
    file,
    line;
    kwargs...,
) = push!(logger.queue, (level, (id, message, module_, group, file, line), kwargs))

function render_logs!(logger)
    MAX_DISPLAYED = 5
    TIMEOUT = 5

    filter!(l -> time() < l[1], logger.displayed)

    n = length(logger.displayed)

    for _ in 1:(MAX_DISPLAYED - n)
        if length(logger.queue) == 0
            break
        end
        push!(logger.displayed, (time() + TIMEOUT, popfirst!(logger.queue)))
    end

    (rows, _) = displaysize(stdout)
    Cursor.save() do
        Cursor.move!(Cursor.Coordinate(rows - MAX_DISPLAYED, 1))
        Screen.clear!(Screen.FromCursorDown())
        for (_, data) in logger.displayed
            println(data)
        end
    end
end

const LOGGER = TUILogger()
global_logger(LOGGER)

Screen.alternative() do
    for event in Inputs.EventLoop(Inputs.Second(1))
        render_logs!(LOGGER)

        if event == Inputs.CTRL_C
            break
        elseif event isa Inputs.Modified && Inputs.Ctrl() in event.modifiers
            @error "D:" event
        elseif event isa Inputs.Character
            @warn ":|" event
        elseif event isa Inputs.Key
            Cursor.move!(Cursor.Coordinate(1, 1))
            Screen.clear!(Screen.CurrentLine())
            println(event)
        end
    end
end


text:
using Documenter, Literate
using PhaseSlopeIndex

EXAMPLE_DIR = joinpath(@__DIR__, "literate")
OUT_DIR = joinpath(@__DIR__, "src/generated")

## Use Literate.jl to generate docs and notebooks of examples
for example in readdir(EXAMPLE_DIR)
    EXAMPLE = joinpath(EXAMPLE_DIR, example)

    Literate.markdown(EXAMPLE, OUT_DIR; documenter=true) # markdown for Documenter.jl
    Literate.notebook(EXAMPLE, OUT_DIR) # .ipynb notebook
    Literate.script(EXAMPLE, OUT_DIR) # .jl script
end

## Build docs
makedocs(;
    sitename="PhaseSlopeIndex.jl",
    format=Documenter.HTML(),
    modules=[PhaseSlopeIndex],
    pages=["Home" => "index.md", "Examples" => "generated/examples.md"],
)

deploydocs(;
    repo="github.com/ssnio/PhaseSlopeIndex.jl.git", devbranch="master", branch="gh-pages"
)


text:
"""
    CurrentStep(num_steps, fmt, lc, header)

Trace column display the current step and total number of iteration,
e.g. `[ 4/10]`, with the figures colored according to progress made.
"""
struct CurrentStep <: TraceColumn
    num_steps::Integer
    fmt::FormatExpr
    lc::LinearColorant{Int}
    header::String
end

"""
    CurrentStep(num_step[; lc, header])

Create a current step trace column for `num_step` iterations. The
default color scale goes from red to green, and the default header is
blank.
"""
function CurrentStep(num_steps::Integer;
                     lc::LinearColorant=LinearColorant(1, num_steps, red_green_scale()),
                     header::String=repeat(" ", 2length(digits(num_steps))+3))
    nd = length(digits(num_steps))
    spacer = repeat(" ", max(0,length(header)-2nd-3))
    CurrentStep(num_steps,
                FormatExpr("[{1:s}{2:$(nd)d}/{3:$(nd)d}$(crayon"reset")]$spacer"),
                lc, rpad(header, 2nd+3))
end
(c::CurrentStep)(i::Integer) = (c.lc(i), i, c.num_steps)

export CurrentStep


text:
using Documenter
using TianJun

include("common.jl")
recursive_replace(@__DIR__, r"\.en(?<ext>(\.md)?)$" => s"\g<ext>")

const BUILD_DIR = "build.en"

makedocs(
    modules = [TianJun],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-132847825-3",
        footer = "Powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl) and the [Julia Programming Language](https://julialang.org/). All contents published at this site follows [CC-BY-4.0](https://creativecommons.org/licenses/by/4.0/) by default. For the Chinese version, please visit [tianjun.me](https://tianjun.me).",
        lang = "en-US",
        assets = [
            "assets/favicon.ico",
            "assets/custom.css",
        ]
    ),
    sitename = "Jun Tian",
    linkcheck = haskey(ENV, "CI"),
    build = BUILD_DIR,
    pages = [
        "👋 About" => "index.md",
        "💻 Programming" => [
            "A Deep Dive into Distributed.jl" => "programming/A_Deep_Dive_into_Distributed.jl/index.md"
        ],
        "📖 Reading" => [
            "Notes on Distributional Reinforcement Learning" => "reading/Notes_on_Distributional_Reinforcement_Learning/index.md"
        ],
        "🙋 Ask Me Anything" => "AMA.md",
        "🔗 Blogroll" => "blogroll.md",
    ]
)

cp(joinpath(@__DIR__, "CNAME"), joinpath(@__DIR__, BUILD_DIR, "CNAME"); force = true)


text: nldata = Dict(:y=>[2.0, 1.0, -0.5, 3.0, 0.25])


text:
module PiecewiseFunctions
import MultivariatePolynomials: AbstractVariable, variable
export PiecewiseFunction, from_above, from_below, pieces, breaks

struct PiecewiseFunction{F, B<:Real, FV<:AbstractVector{F}, BV<:AbstractVector{B}} <: Function
    pieces::FV
    breaks::BV

    function PiecewiseFunction{F, B, FV, BV}(pieces::FV, breaks::BV) where {F, B, FV<:AbstractVector{F}, BV<:AbstractVector{B}}
        @assert issorted(breaks)
        @assert length(breaks) == length(pieces) + 1
        new{F, B, FV, BV}(pieces, breaks)
    end
end

PiecewiseFunction(pieces::FV, breaks::BV) where {F, B, FV<:AbstractVector{F}, BV<:AbstractVector{B}} =
    PiecewiseFunction{F, B, FV, BV}(pieces, breaks)

breaks(p::PiecewiseFunction) = p.breaks
pieces(p::PiecewiseFunction) = p.pieces
(pf::PiecewiseFunction)(t) = from_above(pf, t)

function from_above(pf::PiecewiseFunction, t)
    i = searchsortedlast(pf.breaks, t)
    if i <= 0 || i >= length(pf.breaks)
        error("Input value $t is out of the allowable range [$(pf.breaks[1]), $(pf.breaks[end]))")
    end
    pf.pieces[i](t - pf.breaks[i])
end

function from_below(pf::PiecewiseFunction, t)
    i = searchsortedfirst(pf.breaks, t)
    if i <= 1 || i >= length(pf.breaks) + 1
        error("Input value $t is out of the allowable range ($(pf.breaks[1]), $(pf.breaks[end])]")
    end
    pf.pieces[i - 1](t - pf.breaks[i - 1])
end

Base.broadcast(f, pf::PiecewiseFunction) = PiecewiseFunction(f.(pieces(pf)), breaks(pf))
end


text:
using HTTP
using EzXML
using DataFrames
using Dates

function extract(date)
    result = DataFrame(
        date = String[],
        pitcherid = String[],
        pitcher_teamid = String[],
        pitcher_firstname = String[],
        pitcher_lastname = String[],
        pitcher_teamname = String[],
        pitcherthrow = String[],
        batterid = String[],
        batter_teamid = String[],
        batter_firstname = String[],
        batter_lastname = String[],
        batter_teamname = String[],
        batterstand = String[],
        eventdesc = String[],
        pitchresult = String[],
        x = String[],
        y = String[],
        px = String[],
        pz = String[],
        pfxx = String[],
        pfxz = String[],
        zone = String[],
        sztop = String[],
        szbottom = String[],
        pitchtype = String[],
        startspeed = String[],
        endspeed = String[],
        spindir = String[],
        spinrate = String[],
        nasty = String[]
    )
    year = Dates.year(date)
    month = lpad(Dates.month(date), 2, "0")
    day = lpad(Dates.day(date), 2, "0")
    base = "https://gd2.mlb.com/components/game/mlb/year_$year/month_$month/day_$day"
    gidlist = gethtml(base)
    for li in findall("//li/a/text()",gidlist)
        litext = strip(li.content)
        if occursin(r"gid_*", litext)
            @info "Process dataset: " date "/" litext ": Start"
            playerdf = getplayers(base, litext, date)
            pitchdf = getpitches(base, litext, date)
            processdf = join(pitchdf, playerdf, on=(:pitcherid, :id), kind=:inner, makeunique=true)
            df = join(processdf, playerdf, on=(:batterid, :id), kind=:inner, makeunique=true)
            df = rename(df, [
                    :firstname => :pitcher_firstname,
                    :lastname => :pitcher_lastname,
                    :teamname => :pitcher_teamname,
                    :teamid => :pitcher_teamid,
                    :firstname_1 => :batter_firstname,
                    :lastname_1 => :batter_lastname,
                    :teamname_1 => :batter_teamname,
                    :teamid_1 => :batter_teamid,
                ]
            )
            result = vcat(result, df[
                [
                    :date,
                    :pitcherid,
                    :pitcher_teamid,
                    :pitcher_firstname,
                    :pitcher_lastname,
                    :pitcher_teamname,
                    :pitcherthrow,
                    :batterid,
                    :batter_teamid,
                    :batter_firstname,
                    :batter_lastname,
                    :batter_teamname,
                    :batterstand,
                    :eventdesc,
                    :pitchresult,
                    :x,
                    :y,
                    :px,
                    :pz,
                    :pfxx,
                    :pfxz,
                    :zone,
                    :sztop,
                    :szbottom,
                    :pitchtype,
                    :startspeed,
                    :endspeed,
                    :spindir,
                    :spinrate,
                    :nasty
                ]
            ])
            @info "Process dataset: " date "/" litext ": Finish!"
        end
    end
    return result
end

function gethtml(url)
    r = HTTP.get(url)
    if 200 <= r.status < 300
        return root(parsehtml(String(r.body)))
    else
        error("Page is not accessable.")
    end
end

function getxml(url)
    r = HTTP.get(url)
    if 200 <= r.status < 300
        return root(parsexml(String(r.body)))
    else
        error("Page is not accessable.")
    end
end

function getplayers(base, gid, date)
    playerdf = DataFrame(
        date = String[],
        id = String[],
        firstname = String[],
        lastname = String[],
        teamid = String[],
        teamname = String[]
    )
    url = base * "/" * gid * "players.xml"
    players = getxml(url)
    for player in findall("//player", players)
        try
            id = player["id"]
            firstname = player["first"]
            lastname = player["last"]
            teamid = player["team_id"]
            teamname = player["team_abbrev"]
            push!(playerdf, [Dates.format(date, "yyyy-mm-dd"), id, firstname, lastname, teamid, teamname])
        catch
            @warn "player tag " player " is not extractable."
        end
    end
    return playerdf
end

function getpitches(base, gid, date)
    pitchdf = DataFrame(
        date = String[],
        pitcherid = String[],
        batterid = String[],
        pitcherthrow = String[],
        batterstand = String[],
        eventdesc = String[],
        pitchresult = String[],
        x = String[],
        y = String[],
        px = String[],
        pz = String[],
        pfxx = String[],
        pfxz = String[],
        zone = String[],
        sztop = String[],
        szbottom = String[],
        pitchtype = String[],
        startspeed = String[],
        endspeed = String[],
        spindir = String[],
        spinrate = String[],
        nasty = String[]
    )
    url = base * "/" * gid * "inning/inning_all.xml"
    pitches = getxml(url)
    for atbat in findall("//atbat", pitches)
        pitcherid = atbat["pitcher"]
        batterid = atbat["batter"]
        pitcherthrow = atbat["p_throws"]
        batterstand = atbat["stand"]
        eventdesc = atbat["des"]
        for pitch in findall("./pitch", atbat)
            try
                pitchresult = pitch["des"]
                x = pitch["x"]
                y = pitch["y"]
                px = pitch["px"]
                pz = pitch["pz"]
                pfxx = pitch["pfx_x"]
                pfxz = pitch["pfx_z"]
                zone = pitch["zone"]
                sztop = pitch["sz_top"]
                szbottom = pitch["sz_bot"]
                pitchtype = pitch["pitch_type"]
                startspeed = pitch["start_speed"]
                endspeed = pitch["end_speed"]
                spindir = pitch["spin_dir"]
                spinrate = pitch["spin_rate"]
                nasty = pitch["nasty"]
                push!(pitchdf, [
                    Dates.format(date, "yyyy-mm-dd"),
                    pitcherid,
                    batterid,
                    pitcherthrow,
                    batterstand,
                    eventdesc,
                    pitchresult,
                    x,
                    y,
                    px,
                    pz,
                    pfxx,
                    pfxz,
                    zone,
                    sztop,
                    szbottom,
                    pitchtype,
                    startspeed,
                    endspeed,
                    spindir,
                    spinrate,
                    nasty
                    ]
                )
            catch
               @warn "pitch tag " pitch " is not extractable."
            end
        end
    end
    return pitchdf
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

pos_training_path = joinpath(main_image_path, "trainset", "faces")
neg_training_path = joinpath(main_image_path, "trainset", "non-faces")
pos_testing_path = joinpath(main_image_path, "testset", "faces")#joinpath(homedir(), "Desktop", "faces")#"$main_image_path/testset/faces/"
neg_testing_path = joinpath(main_image_path, "testset", "non-faces")

# pos_training_path = joinpath(data_path, "lfw-all")
# neg_training_path = joinpath(data_path, "all-non-faces")
# pos_testing_path = joinpath(data_path, "lizzie-testset", "faces")
# neg_testing_path = joinpath(data_path, "lizzie-testset", "nonfaces")


text:
using Counters
using Test

# write your own tests here
@test 1 == 1
c = Counter{String}()
c["alpha"] = 1
incr!(c, "bravo")
c["charlie"] += 1

@test c["alpha"] == 1
@test c["bravo"] == 1
@test c["charlie"] == 1

c = c + c

@test c["alpha"] == 2
@test nnz(c) == 3
@test sum(c) == 6

c["delta"] = 2
c["delta"] += -10
@test c["delta"] == 0

d = deepcopy(c)
clean!(d)
@test c == d

@test nnz(c) == 3
@test sum(c) == 6
@test length(c) == 4
@test length(d) == 3

X = Set(1:3)
c = counter(X)
Y = BitSet(1:3)
d = counter(Y)
@test c == d

c = counter([1, 1, 2, 3])
@test mean(c) == 1.75

x = [mod(t, 10) for t in rand(Int, 100)]
c = counter(x)
xx = collect_by_counts(c)
top = xx[1]
@test c[top] == maximum(values(c))


text:
@testset "semi_safe_load" begin

@test AMDGPU.PAGESIZE == 0x1000

struct MediumObject end
Base.sizeof(::MediumObject) = 0x1000
Base.sizeof(::Type{MediumObject}) = 0x1000
struct BigObject end
Base.sizeof(::BigObject) = 0x2000
Base.sizeof(::Type{BigObject}) = 0x2000

for (ptr, bases) in (
    # Less-than-page object
    (Ptr{Int}(UInt64(0x0)), UInt64[0x0]),
    (Ptr{Int}(UInt64(0x1)), UInt64[0x0]),
    (Ptr{Int}(UInt64(0xffc)), UInt64[0x0, 0x1000]),
    (Ptr{Int}(UInt64(0x1000)), UInt64[0x1000]),
    (Ptr{Int}(UInt64(0x1001)), UInt64[0x1000]),

    # Full-page object
    (Ptr{MediumObject}(UInt64(0x0)), UInt64[0x0]),
    (Ptr{MediumObject}(UInt64(0x1)), UInt64[0x0, 0x1000]),
    (Ptr{MediumObject}(UInt64(0xffc)), UInt64[0x0, 0x1000]),
    (Ptr{MediumObject}(UInt64(0x1000)), UInt64[0x1000]),
    (Ptr{MediumObject}(UInt64(0x1001)), UInt64[0x1000, 0x2000]),

    # Multi-page object
    (Ptr{BigObject}(UInt64(0x0)), UInt64[0x0, 0x1000]),
    (Ptr{BigObject}(UInt64(0x1)), UInt64[0x0, 0x1000, 0x2000]),
    (Ptr{BigObject}(UInt64(0xffc)), UInt64[0x0, 0x1000, 0x2000]),
    (Ptr{BigObject}(UInt64(0x1000)), UInt64[0x1000, 0x2000]),
    (Ptr{BigObject}(UInt64(0x1001)), UInt64[0x1000, 0x2000, 0x3000]),
)
    @test AMDGPU._pages_spanned(ptr) == length(bases)
    AMDGPU._per_page((x,y)->nothing, ptr) do base
        @test base in bases
    end
end

end


text:
using CLIMA.PlanetParameters
export NoRadiation, StevensRadiation

abstract type RadiationModel end

vars_state(::RadiationModel, FT) = @vars()
vars_aux(::RadiationModel, FT) = @vars()
vars_integrals(::RadiationModel, FT) = @vars()

function atmos_nodal_update_aux!(::RadiationModel, ::AtmosModel, state::Vars,
                                 aux::Vars, t::Real)
end
function preodefun!(::RadiationModel, aux::Vars, state::Vars, t::Real)
end
function integrate_aux!(::RadiationModel, integ::Vars, state::Vars, aux::Vars)
end
function flux_radiation!(::RadiationModel, flux::Grad, state::Vars,
                         aux::Vars, t::Real)
end

struct NoRadiation <: RadiationModel
end

"""
  StevensRadiation <: RadiationModel

Stevens et. al (2005) version of the δ-four stream model used to represent radiative transfer.
Analytical description as a function of the liquid water path and inversion height zᵢ
"""
struct StevensRadiation{FT} <: RadiationModel
  "κ [m^2/s] "
  κ::FT
  "α_z Troposphere cooling parameter [m^(-4/3)]"
  α_z::FT
  "z_i Inversion height [m]"
  z_i::FT
  "ρ_i Density"
  ρ_i::FT
  "D_subsidence Large scale divergence [s^(-1)]"
  D_subsidence::FT
  "F₀ Radiative flux parameter [W/m^2]"
  F_0::FT
  "F₁ Radiative flux parameter [W/m^2]"
  F_1::FT
end
vars_integrals(m::StevensRadiation, FT) = @vars(∂κLWP::FT)
function integrate_aux!(m::StevensRadiation, integrand::Vars, state::Vars, aux::Vars)
  FT = eltype(state)
  integrand.radiation.∂κLWP = state.ρ * m.κ * aux.moisture.q_liq
end
function flux_radiation!(m::StevensRadiation, flux::Grad, state::Vars,
                         aux::Vars, t::Real)
  FT = eltype(flux)
  z = aux.orientation.Φ/grav
  Δz_i = max(z - m.z_i, -zero(FT))
  # Constants
  cloud_top_cooling = m.F_0 * exp(-aux.∫dnz.radiation.∂κLWP)
  cloud_base_warming = m.F_1 * exp(-aux.∫dz.radiation.∂κLWP)
  free_troposphere_cooling = m.ρ_i * FT(cp_d) * m.D_subsidence * m.α_z * ((cbrt(Δz_i))^4 / 4 + m.z_i * cbrt(Δz_i))
  F_rad = cloud_base_warming + cloud_base_warming + free_troposphere_cooling
  flux.ρe += SVector(FT(0),
                     FT(0),
                     F_rad)
end
function preodefun!(m::StevensRadiation, aux::Vars, state::Vars, t::Real)
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using Distributions
using Random
using LinearAlgebra: diagm
using SMC

struct LinearGaussianInfer{T}
    α::T
    K::Int
    Y::Matrix{Float64}
    end_t::Int
end

function LinearGaussianInfer(K, Y, end_t)
    sigma_dist = LogNormal(log(9.0), 0.1)
    tau_dist = LogNormal(log(1.0), 0.1)
    return LinearGaussianInfer((σ = sigma_dist, τ = tau_dist), K, Y, end_t)
end

function rprior(m::LinearGaussianInfer, K)
    [map(rand, m.α) for _ = 1:K]
end

function dprior(m::LinearGaussianInfer, p)
    sum(map(logpdf, m.α, p))
end

function dprior(m::LinearGaussianInfer, ps::Array)
    [dprior(m, p) for p in ps]
end

function loglik_fn(m::LinearGaussianInfer, p)
    F0(t) = [0.0, 0.0, 0.0]
    G0(t) = [0.0, 0.0, 0.0]
    F(t) = diagm([1.0, 1.0, 1.0])
    G(t) = diagm(([1.0, 1.0, 1.0]))
    Σ(t) = p.σ * diagm(([1.0, 1.0, 1.0]))
    Tau(t) = diagm([p.τ, p.τ, p.τ])
    μ0 = vcat(0.0, 0.0, 0.0)
    Tau0 = 100.0 * diagm([1.0, 1.0, 1.0])

    smc_fns = lgc_smc_model(m.K, F0, F, G0, G, Tau, Σ, μ0, Tau0, m.Y)

    smc_filter = SMC.SMCModel(
        record_history = true,
        rinit = smc_fns.rinit,
        rproposal = smc_fns.rp,
        dinit = smc_fns.dinit,
        dprior = smc_fns.dpr,
        dproposal = smc_fns.dp,
        dtransition = smc_fns.dt,
        dmeasure = smc_fns.dm,
        dpre = smc_fns.dpre,
        threshold = 1.0,
    )

    return smc_filter(Random.GLOBAL_RNG, m.end_t).loglik
end

function (m::LinearGaussianInfer)(p)
    loglik_fn(m, p)
end


text:
using Test
using NFFT
using Random
using LinearAlgebra
using FFTW
using NFFTTools

Random.seed!(123)

include("accuracy.jl")
include("constructors.jl")
include("performance.jl")
include("testToeplitz.jl")
include("samplingDensity.jl")
include("cuda.jl")
# Need to run after the other tests since the overload plan_*
include("wrappers.jl")


text:
abstract type Agent end
abstract type Parameters end
abstract type SimulationOptions end
abstract type SimulationResults end
abstract type SimulationState end

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
@testset "constructor <from connection>"  begin
    io = IOBuffer(httptest1)
    c = Connection(io)
    Spirit.processrequest!(c)

    req = Spirit.Request(c)

    @test req.method == "POST"
    @test req.uri == "/"
    @test req.httpversion == "1.1"
    @test length(req.headers) == 2
    @test get(req.headers, "Host") == "www.somewhere.com"
    @test get(req.headers, "Content-Length") == "9"
    @test req.resource == "/"
    @test req.resource_base == ""

    close(io)

    io = IOBuffer(httptest2)
    c = Connection(io)
    Spirit.processrequest!(c)

    req = Spirit.Request(c)

    @test req.method == "GET"
    @test req.uri == "/docs/index.html"
    @test req.httpversion == "1.1"
    @test length(req.headers) == 5
    @test get(req.headers, "Host") == "www.nowhere123.com"
    @test get(req.headers, "Accept") == "image/gif, image/jpeg, */*"
    @test get(req.headers, "Accept-Language") == "en-us"
    @test get(req.headers, "Accept-Encoding") == "gzip, deflate"
    @test get(req.headers, "User-Agent") == "Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1)"
    @test req.resource == "/docs/index.html"
    @test req.resource_base == ""

    close(io)
end

@testset "getheader" begin
    io = IOBuffer(httptest1)
    c = Connection(io)
    Spirit.processrequest!(c)

    req = Spirit.Request(c)

    @test Spirit.getheader(req, "Host") == "www.somewhere.com"
    @test Spirit.getheader(req, "Content-Length") == "9"
    @test Spirit.getheader(req, "Does Not Exist") == nothing
    @test Spirit.getheader(req, "Does Not Exist", "sasquatch") == "sasquatch"

    close(io)
end

@testset "setheader" begin
    io = IOBuffer(httptest1)
    c = Connection(io)
    Spirit.processrequest!(c)

    req = Spirit.Request(c)

    req2 = Spirit.setheader(req, "X-Special-Thing", "some value")
    @test Spirit.getheader(req, "X-Special-Thing") == nothing
    @test Spirit.getheader(req2, "X-Special-Thing") == "some value"

    req3 = Spirit.setheader(req2, "X-Special-Thing", nothing)
    @test Spirit.getheader(req2, "X-Special-Thing") == "some value"
    @test Spirit.getheader(req3, "X-Special-Thing") == nothing
    @test Spirit.getheader(req3, "X-Special-Thing", "other value") == "other value"

    close(io)
end

@testset "appendheader" begin
    io = IOBuffer(httptest1)
    c = Connection(io)
    Spirit.processrequest!(c)

    req = Spirit.Request(c)

    req2 = Spirit.appendheader(req, "X-Special-Thing", "some value")
    @test Spirit.getheader(req, "X-Special-Thing") == nothing
    @test Spirit.getheader(req2, "X-Special-Thing") == "some value"

    req3 = Spirit.appendheader(req2, "X-Special-Thing", "some other value")
    @test Spirit.getheader(req2, "X-Special-Thing") == "some value"
    @test Spirit.getheader(req3, "X-Special-Thing") == "some value, some other value"

    close(io)
end


text:
#md # !!! note
#
#     This tutorial is under construction, but the code below is already functional.
#
# Driver that computes the lid-driven cavity benchmark at low Reynolds numbers
# when using a mixed FE Q(k)/Pdisc(k-1).

# Load Gridap library
using Gridap

# Discrete model
n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)

# Define Dirichlet boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])

# Define test FESpaces (Q2/P1(disc) pair)
V = TestFESpace(
  reffe=:QLagrangian, conformity=:H1, valuetype=VectorValue{2,Float64},
  model=model, labels=labels, order=2, dirichlet_tags=["diri0","diri1"])
Q = TestFESpace(
  reffe=:PLagrangian, conformity=:L2, valuetype=Float64,
  model=model, order=1, constraint=:zeromean)
Y = MultiFieldFESpace([V,Q])

# Define trial FESpaces from Dirichlet values
u0 = VectorValue(0,0)
u1 = VectorValue(1,0)
U = TrialFESpace(V,[u0,u1])
P = TrialFESpace(Q)
X = MultiFieldFESpace([U,P])

# Define integration mesh and quadrature
trian = get_triangulation(model); degree = 2
quad = CellQuadrature(trian,degree)

# Define and solve the FE problem
function a(x,y)
  v,q = y
  u,p = x
  ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)
end
t_Ω = LinearFETerm(a,trian,quad)
op = AffineFEOperator(X,Y,t_Ω)
uh, ph = solve(op)

# Export results to vtk
writevtk(trian,"results",order=2,cellfields=["uh"=>uh,"ph"=>ph])

max_stars_repo_name: andre-massing/Tutorials

max_issues_repo_name: andre-massing/Tutorials

max_forks_repo_name: andre-massing/Tutorials


text:
module Sentinel

include("./constants.jl")
using JuMP
using .Constants

export init_sentinel

function init_sentinel(m)

    fix(m[:flbx_front_left_bot_x][length(boxes)], 0, force=true)
    fix(m[:flby_front_left_bot_y][length(boxes)], 0, force=true)
    fix(m[:flbz_front_left_bot_z][length(boxes)], 0, force=true)

    fix(m[:lpx_length_of_box_parallel_to_x][length(boxes)], 1, force=true)
    fix(m[:wpy_width_of_box_parallel_to_y][length(boxes)], 1, force=true)
    return m
end

end


text:
module BezierRefFEsTests

using Gridap

using Gridap.ReferenceFEs
using Gridap.Polynomials
using Gridap.Fields
using Gridap.Arrays
using FillArrays

using Test

using Gridap.ReferenceFEs: _berstein_matrix
using Gridap.ReferenceFEs: berstein_basis
using Gridap.ReferenceFEs: rationalize_bernstein_basis

## Test Bernstein basis

p_filter(e,o) = sum(e) ≤ o

# 1D
p = 2
prebasis_seg = MonomialBasis{1}(Float64,p,p_filter)
C = _berstein_matrix(prebasis_seg,SEGMENT)
C12 =
[
 1  -2   1
 0   2  -2
 0   0   1
]

@test transpose(C) == C12

ϕ = berstein_basis(prebasis_seg,SEGMENT)
X = [ Point(0.0,0.0), Point(0.5,0.5), Point(1.0,0.0) ]
ψ = linear_combination(X,ϕ)
ξ = [Point(0.0),Point(0.5),Point(1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == [ Point(0.0,0.0), Point(0.5,0.25), Point(1.0,0.0) ]

p = 3
prebasis_seg = MonomialBasis{1}(Float64,p,p_filter)
C = _berstein_matrix(prebasis_seg,SEGMENT)
C13 =
[
 1  -3   3  -1
 0   3  -6   3
 0   0   3  -3
 0   0   0   1
]

@test transpose(C) == C13

ϕ = berstein_basis(prebasis_seg,SEGMENT)
X = [ Point(0.0,0.0), Point(0.25,0.5), Point(0.75,0.5), Point(1.0,0.0) ]
ψ = linear_combination(X,ϕ)
ξ = [Point(0.0),Point(0.5),Point(1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == [ Point(0.0,0.0), Point(0.5,0.375), Point(1.0,0.0) ]

# 2D

p = 2
prebasis_tri = MonomialBasis{2}(Float64,p,p_filter)
C = _berstein_matrix(prebasis_tri,TRI)
C22 =
[
 1  -2   1  -2   2   1
 0   2  -2   0  -2   0
 0   0   1   0   0   0
 0   0   0   2  -2  -2
 0   0   0   0   2   0
 0   0   0   0   0   1
]

@test transpose(C) == C22

ϕ = berstein_basis(prebasis_tri,TRI)
X = [ Point(0.0,0.0), Point(0.5,0.0), Point(1.0,0.0), Point(0.0,0.5), Point(0.5,0.5), Point(0.0,1.0) ]
ψ = linear_combination(X,ϕ)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(0.3,0.3),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ

p = 3
prebasis_tri = MonomialBasis{2}(Float64,p,p_filter)
C = _berstein_matrix(prebasis_tri,TRI)
C23 =
[
 1  -3   3  -1  -3   6  -3   3  -3  -1
 0   3  -6   3   0  -6   6   0   3   0
 0   0   3  -3   0   0  -3   0   0   0
 0   0   0   1   0   0   0   0   0   0
 0   0   0   0   3  -6   3  -6   6   3
 0   0   0   0   0   6  -6   0  -6   0
 0   0   0   0   0   0   3   0   0   0
 0   0   0   0   0   0   0   3  -3  -3
 0   0   0   0   0   0   0   0   3   0
 0   0   0   0   0   0   0   0   0   1
]

@test transpose(C) == C23

X = [
  Point(0.0,0.0), Point(0.25,0.0), Point(0.75,0.0), Point(1.0,0.0),
  Point(0.0,0.25), Point(0.25,0.25), Point(0.75,0.25),
  Point(0.0,0.75), Point(0.25,0.75),
  Point(0.0,1.0) ]
ϕ = berstein_basis(prebasis_tri,TRI)
ψ = linear_combination(X,ϕ)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(0.5,0.5),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ

## BezierRefFE

tri = BezierRefFE(Float64,TRI,(3,3))
nodes = get_node_coordinates(tri) * 5
ϕ = get_shapefuns(tri)
ψ = linear_combination(nodes,ϕ)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(0.5,0.5),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ.*5

f = ConstantField(1)
q = Quadrature(TRI,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) ≈ (5*5)/2

tri = BezierRefFE(Float64,TRI,(3,3))
nodes = get_node_coordinates(tri)
nodes[10] = Point(0.1,0.1)
ϕ = get_shapefuns(tri)
ψ = linear_combination(nodes,ϕ)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(0.5,0.5),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ

f = ConstantField(1)
q = Quadrature(TRI,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) ≈ 1/2

tri = BezierRefFE(Float64,TRI,(3,3))
nodes = get_node_coordinates(tri)
nodes[4] = Point(1/3,-0.1)
nodes[5] = Point(2/3,-0.1)
ϕ = get_shapefuns(tri)
ψ = linear_combination(nodes,ϕ)
f = ConstantField(1)
q = Quadrature(TRI,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) > 1/2

tet = BezierRefFE(Float64,TET,(3,3,3))
nodes = get_node_coordinates(tet) .* 5
ϕ = get_shapefuns(tet)
ψ = linear_combination(nodes,ϕ)
ξ = [Point(0.0,0.0,0.0),Point(0.0,0.0,0.5),Point(0.5,0.5,0.5),Point(1.0,0.0,0.0),Point(0.0,1.0,0.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ .* 5

f = ConstantField(1)
q = Quadrature(TET,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) ≈ (5*5*5)/6

quad = BezierRefFE(Float64,QUAD,(3,3))
nodes = get_node_coordinates(quad) .* 5
ϕ = get_shapefuns(quad)
ψ = linear_combination(nodes,ϕ)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(1.0,0.5),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ .* 5

f = ConstantField(1)
q = Quadrature(QUAD,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) ≈ (5*5)

hex = BezierRefFE(Float64,HEX,(3,3,3))
nodes = get_node_coordinates(hex) .* 5
ϕ = get_shapefuns(hex)
ψ = linear_combination(nodes,ϕ)
ξ = [Point(0.0,0.0,0.0),Point(0.0,0.0,0.5),Point(0.5,1.0,0.5),Point(1.0,0.0,0.0),Point(0.0,1.0,0.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ .* 5

f = ConstantField(1)
q = Quadrature(HEX,3*2)
p = get_coordinates(q)
w = get_weights(q)
J = ∇(ψ)

@test integrate(f,p,w,J) ≈ (5*5*5)

tri = ReferenceFE(TRI,bezier,Float64,(3,3))
_tri = BezierRefFE(Float64,TRI,(3,3))

@test tri == _tri

tri = BezierRefFE(Float64,TRI,(3,3))
ϕ = get_shapefuns(tri)
w = ones(length(ϕ))
ϕr = rationalize_bernstein_basis(ϕ,w)
nodes = get_node_coordinates(tri)
ψ = linear_combination(nodes,ϕr)
ξ = [Point(0.0,0.0),Point(0.0,0.5),Point(0.5,0.5),Point(1.0,0.0),Point(0.0,1.0)]
Ψ = Fill(ψ,length(ξ))
Xi = lazy_map( evaluate, Ψ, ξ )

@test Xi == ξ

reffe = BezierRefFE(Float64,SEGMENT,(3,))
test_lagrangian_reference_fe(reffe)

reffe = BezierRefFE(Float64,TRI,(3,3))
test_lagrangian_reference_fe(reffe)

reffe = BezierRefFE(Float64,TET,(3,3,3))
test_lagrangian_reference_fe(reffe)

end # module


text:
using DifferentialEquations

f(u,t) = 0.98*u
u0 = 1.0
tspan = (0.0,5.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

using Plots; gr()
plot(sol)

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!")

plot!(sol.t,t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")

###############
# Solving system of Equaitons

function lorenz(du,u,p,t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3] - u[2])
    du[3] = u[1] * u[2] - (8.0/3.0) * u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,1.0)

prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob)
plot(sol , vars=[1,2,3])
# animate(sol,lw=3,every=4)
function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
gr()
sol = solve(prob,vars=[1,2,3])

using StaticArrays

#######################
using DiffEqOperators

order = 2
deriv = 2

Δx =  0.1
N = 9
A = DerivativeOperator{Float64}(order,deriv,Δx,N,:Dirichlet0,:Dirichlet0)

full(DerivativeOperator{Float64}(order,deriv,Δx,N,:Dirichlet0,:Dirichlet0))

max_stars_repo_path: Missing around and learning DifferentialEquations.jl/Juliacon_2017_code.jl


text:
function main(r)
    x=y=0
    koaf=1
    while (!isborder(r,Sud))
        move!(r,Sud)
        y+=1
        koaf=(koaf+1)%2
    end
    while (!isborder(r,West))
        move!(r,West)
        x+=1
        koaf=(koaf+1)%2
    end
    moves(r,Ost,koaf)
    while (!isborder(r,Nord))
        for i in (West,Ost)
            if (!isborder(r,Nord))
                move!(r,Nord)
                moves(r,i,koaf)
            end
        end
    end
    while (!isborder(r,Sud))
        move!(r,Sud)
    end
    while (!isborder(r,West))
        move!(r,West)
    end
    go_back(r,x,y)
end

function moves(r, side,koaf::Int)
    if (koaf==1)
        putmarker!(r)
    end
    while (!isborder(r,side))
        move!(r,side)
        koaf=(koaf+1)%2
        if (koaf==1)
            putmarker!(r)
        end
    end
end

function go_back(r, x::Int , y::Int)
    for i in 1:x
        move!(r, Ost)
    end
    for i in 1:y
        move!(r, Nord)
    end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
using FiniteElement
using Test

@testset "FiniteElement.jl" begin
    # Write your tests here.
end


text:
# ---
# title: 1032. Stream of Characters
# id: problem1032
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Trie
# link: <https://leetcode.com/problems/stream-of-characters/description/>
# hidden: true
# ---
#
# Implement the `StreamChecker` class as follows:
#
#   * `StreamChecker(words)`: Constructor, init the data structure with the given words.
#   * `query(letter)`: returns true if and only if for some `k >= 1`, the last `k` characters queried (in order from oldest to newest, including this letter just queried) spell one of the words in the given list.
#
#
#
# **Example:**
#
#
#
#     StreamChecker streamChecker = new StreamChecker(["cd","f","kl"]); // init the dictionary.
#     streamChecker.query('a');          // return false
#     streamChecker.query('b');          // return false
#     streamChecker.query('c');          // return false
#     streamChecker.query('d');          // return true, because 'cd' is in the wordlist
#     streamChecker.query('e');          // return false
#     streamChecker.query('f');          // return true, because 'f' is in the wordlist
#     streamChecker.query('g');          // return false
#     streamChecker.query('h');          // return false
#     streamChecker.query('i');          // return false
#     streamChecker.query('j');          // return false
#     streamChecker.query('k');          // return false
#     streamChecker.query('l');          // return true, because 'kl' is in the wordlist
#
#
#
#
# **Note:**
#
#   * `1 <= words.length <= 2000`
#   * `1 <= words[i].length <= 2000`
#   * Words will only consist of lowercase English letters.
#   * Queries will only consist of lowercase English letters.
#   * The number of queries is at most 40000.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module MultiObjectiveProblems

import JuMP
using OSQP
import ForwardDiff
using Parameters: @with_kw

import NLopt;    # for ZDT4 local critical set

export MOP, SamplingFunction, FixedPointSet, Constraints, Box
export num_objectives, num_vars;
export get_objectives, get_vector_objective, get_gradients, get_omega_function, constraints,
    get_pareto_set, get_pareto_front, get_points, get_scatter_arrays, get_scatter_points,
    get_random_point, get_ideal_point, get_critical_set, get_critical_front,
    get_starting_point;
export sampling_methods;

abstract type Constraints end;
struct Box <:Constraints
    lb :: Vector{Float64}
    ub :: Vector{Float64}
end

@doc """
        abstract type MOP end;

    A supertype for all problem types.
    Each problem must implement the following methods:
    * `num_vars`
    * `num_objectives`
    * `get_objectives`
"""
abstract type MOP end;
Broadcast.broadcastable(m::M where M<:MOP) = Ref(m)

# required
@doc "Return on array of multivariate objective functions."
function get_objectives( :: MOP ) end;
@doc "Return the number of decision variables."
function num_vars( :: MOP) end;
@doc "Return the number of objective functions."
function num_objectives( :: MOP ) end;

# non-mandatory but generically implemented

@doc "Return a constraint object for `mop` or nothing."
function constraints( mop :: M where M <: MOP ) end;
function get_vector_objective( mop :: M where M<:MOP )
    func_array = get_objectives( mop );
    return (x::Vector{R} where R<:Real) -> [ f(x) for f in func_array ];
end

function get_gradients( mop :: M where M <: MOP )
    return [ (x :: Vector{R} where R<:Real ) -> ForwardDiff.gradient( f, x )
        for f in get_objectives(mop) ]
end

function get_omega_function( mop :: M where M <: MOP)
    omega_function = function( x :: Vector{R} where {R<:Real} )
        try
            ∇F = transpose( hcat( ( ∇f(x) for ∇f in get_gradients( mop ) )... ) );

            prob = JuMP.Model( OSQP.Optimizer );

            JuMP.set_silent(prob);
            JuMP.set_optimizer_attribute(prob,"eps_rel",1e-5);
            JuMP.set_optimizer_attribute(prob,"eps_abs", 1e-5);
            JuMP.set_optimizer_attribute(prob,"polish",true);

            JuMP.@variable(prob, α );     # negative of marginal problem value
            JuMP.@variable(prob, d[1:num_vars(mop)] );   # direction vector

            JuMP.@objective(prob, Min, α);

            JuMP.@constraint(prob, descent_contraints, ∇F*d .<= α);
            JuMP.@constraint(prob, norm_constraints, -1.0 .<= d .<= 1.0);

            global_constraints = constraints(mop);
            if isa( global_constraints, Box )
                JuMP.@constraint(prob, box_constraints,  global_constraints.lb .<= x .+ d .<= global_constraints.ub );
            end

            JuMP.optimize!(prob)
            ω = -JuMP.value(α);
            return ω
        catch e
            return -Inf
        end
    end
    return omega_function
end

function get_random_point( mop:: M where M <: MOP )
    if isa( constraints(mop), Box )
        lb = constraints(mop).lb;
        w = constraints(mop).ub .- lb;
        return lb .+ w .* rand( length(lb) );
    end
end

@doc """
Get a starting point for the problem.
If the problem `mop` does not specify its own function, then
a random point is returned for box constrained problems and `nothing` else.
"""
get_starting_point( mop :: M where M<: MOP, args...; kwargs...) = get_random_point(mop);

# non-mandatory and not implemented by default
@doc "Return the ideal (or an utopia) point of the problem."
function get_ideal_point( ::MOP ) end;

@doc """
    get_pareto_set( mop )

Get an object representing the Pareto Set of `mop`.
Possible return types are `FixedPointSet` and `SamplingFunction`.

See also: [`get_points`](@ref), [`get_scatter_points`](@ref), [`SamplingFunction`](@ref),
[`FixedPointSet`](@ref)
"""
function get_pareto_set( :: MOP ) end;

@doc """
    get_pareto_front( mop )

Get an object representing the Pareto Front of `mop`.
Possible return types are `FixedPointSet` and `SamplingFunction`.

See also: [`get_points`](@ref), [`get_scatter_points`](@ref), [`SamplingFunction`](@ref),
[`FixedPointSet`](@ref)
"""
function get_pareto_front( :: MOP ) end;

function get_critical_set( :: MOP ) end;
function get_critical_front( :: MOP ) end;

@doc "Abstract Super Type for Pareto Set and Pareto Frontier."
abstract type CompareSet end;
function get_points( ::CompareSet, ::Int; kwargs...) end;

@doc "A wrapper around a `list_of_points`."
struct FixedPointSet <: CompareSet
    list_of_points :: Vector{Vector{R}} where R<:Real;
end;     # for calculated point clouds

@doc """
    A `SamplingFunction` enables you to sample points from the Pareto Set
    or the Pareto Front.
    The basic usage for `samp_func :: SamplingFunction` is
    ```
    num_points = 10;
    list_of_points = get_points( samp_func, num_points );
    ```

    Check `sampling_methods(samp_func)` to see whether multiple sampling techniques
    are available.
    You can pass a sampling method via `method` keyword argument:
    ```
    list_of_points = get_points( samp_func, num_points; method = :regular );
    ```

    See also: [`sampling_methods`](@ref), [`get_points`](@ref), [`get_scatter_points`](@ref),
"""
@with_kw struct SamplingFunction <: CompareSet
    mop :: M where M <: MOP;
    type :: Symbol;
    methods :: Vector{Symbol};

    @assert type ∈ [:ParetoSet, :ParetoFront, :ParetoCriticalSet, :ParetoCriticalFront]
end

@doc "Return a list of Symbols that can be pased to `get_points` for sampling from
a `SamplingFunction`."
sampling_methods( samp_func :: SamplingFunction ) = samp_func.methods;

@doc """
    get_points( samp_func, N; kwargs... )

Sample `N` points using `samp_func::SamplingFunction`.
Additional `kwargs` are passed down to the problem-specfic implementations.
"""
function get_points( samp_func :: SamplingFunction, N :: Int; kwargs... )
    return get_points( samp_func, samp_func.mop, Val(samp_func.type), N; kwargs...);
end

@doc "
    get_points( point_set, N; kwargs... )

Get all points using `point_set::FixedPointSet`."
function get_points( point_set :: FixedPointSet, N :: Int; kwargs... )
    println("Retrieving *all* points from a `FixedPointSet`, all other arguments are ignored.")
    return point_set.list_of_points;
end

@doc "Like `get_points` but returns data suited for plotting.\n
Specify variable axes by kwarg `dims` which default to `[1,2]`.

See also: [`get_points`](@ref)"
function get_scatter_points( samp_func :: SamplingFunction, N :: Int; dims = [1,2], kwargs... )
    get_scatter_arrays( get_points( samp_func, N; kwargs... ); dims = dims )
end

# not called directly
function get_points( samp_func :: SamplingFunction, mop :: M where M<:MOP,
        ::Val{:ParetoFront}, N; kwargs... )
    pset = get_points( samp_func, samp_func.mop, Val(:ParetoSet), N; kwargs...);
    F = get_vector_objective( mop );
    return F.(pset)
end

# not called directly
function get_points( samp_func :: SamplingFunction, mop :: M where M<:MOP,
    ::Val{:ParetoCriticalFront}, N; kwargs... )
pset = get_points( samp_func, samp_func.mop, Val(:ParetoCriticalSet), N; kwargs...);
F = get_vector_objective( mop );
return F.(pset)
end

# also not called directly
@doc "Return x and y data suited for scatter plots when `arr` is a list of vectors."
function get_scatter_arrays( arr :: Vector{Vector{R}} where R <: Real;
        dims :: Vector{Int} = [1, 2] )
    X = [ point[dims[1]] for point in arr ];
    Y = [ point[dims[2]] for point in arr ];
    if length( dims ) == 2
        return X, Y
    else
        Z = [ point[dims[3]] for point in arr ];
        return X,Y,Z
    end
end

include("two_parabolas.jl");
include("ZDT_problems.jl");
include("DTLZ_problems.jl");
include("lovison_problems.jl")
include("MHT_problems.jl")

end


text:
# This file was generated, do not modify it. # hide
train_rows = vcat(1:60, 91:150); # some row indices (observations are rows not columns)
fit!(mach, rows=train_rows)
fitted_params(mach)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export DiscountedReward, discounted

"Discounts reward of underlying specification."
struct DiscountedReward{S <: Specification} <: Specification
    spec::S # Underlying specification
    discount::Float64 # Discount factor
end

DiscountedReward(spec) = DiscountedReward(spec, 1.0, 1.0)
DiscountedReward(spec::DiscountedReward, discount) =
    DiscountedReward(spec.spec, discount * spec.discount)

Base.hash(spec::DiscountedReward, h::UInt) =
    hash(spec.discount, hash(spec.spec, h))
Base.isequal(s1::DiscountedReward, s2::DiscountedReward) =
    s1.discount == s2.discount && s1.spec == s2.spec

is_goal(spec::DiscountedReward, domain::Domain, state::State) =
    is_goal(spec.spec, domain, state)
is_violated(spec::DiscountedReward, domain::Domain, state::State) =
    is_violated(spec.spec, domain, state)
get_cost(spec::DiscountedReward, domain::Domain, s1::State, a::Term, s2::State) =
    get_cost(spec.spec, domain, s1, a, s2)
get_reward(spec::DiscountedReward, domain::Domain, s1::State, a::Term, s2::State) =
    get_reward(spec.spec, domain, s1, a, s2)
get_goal_terms(spec::DiscountedReward) =
    get_goal_terms(spec.spec)

"Discount the rewards or costs associated with `spec`."
discounted(spec::Specification, discount::Float64) =
    DiscountedReward(spec, discount)


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

fname = scm_path("..", "examples", "Shipley", "fig2.6", "fig2.6.a.dot")
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

# Read `=>` as `~` in regression models, or `<-` in causal models.
d = OrderedDict(
  :u => [:x, :v, :s1],
  :w => [:v, :y, :s2]
);

dag = DAG("fig2_6_a", d);
show(dag)

f = [:x]; s = [:v]
g1 = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) == $g1")
g1 = d_separation(dag, f, s; c=[:u])
println("d_separation($(dag.name), $f, $s; c=[:u) == $g1")
g1 = d_separation(dag, f, s; c=[:s1])
println("d_separation($(dag.name), $f, $s; c=[:s1]) == $g1\n")

f = [:x]; s = [:w]
g2 = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) == $g2")
g2 = d_separation(dag, f, s; c=[:u])
println("d_separation($(dag.name), $f, $s; c=[:u]) == $g2")
g2 = d_separation(dag, f, s; c=[:v])
println("d_separation($(dag.name), $f, $s; c=[:v]) == $g2")
g2 = d_separation(dag, f, s; c=[:u, :v])
println("d_separation($(dag.name), $f, $s; c=[:u, :v]) == $g2\n")

f = [:x]; s = [:y]
g3 = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) == $g3")
g3 = d_separation(dag, f, s; c=[:u, :w])
println("d_separation($(dag.name), $f, $s; c=[:u, :w]) == $g3")
g3 = d_separation(dag, f, s; c=[:s1, :s2])
println("d_separation($(dag.name), $f, $s; c=[:s1, :s2]) == $g3")
g3 = d_separation(dag, f, s; c=[:u, :v, :w])
println("d_separation($(dag.name), $f, $s; c=[:u, :v, :w]) == $g3")
g3 = d_separation(dag, f, s; c=[:s1, :s2, :u])
println("d_separation($(dag.name), $f, $s; c=[:s1, :s2, :u]) == $g3")
g3 = d_separation(dag, f, s; c=[:s1])
println("d_separation($(dag.name), $f, $s; c=[:s1]) == $g3\n")

f = [:u]; s = [:w]
g4 = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) == $g4"), [:s1]
g4 = d_separation(dag, f, s; c=[:v])
println("d_separation($(dag.name), $f, $s; c=[:v]) == $g4\n")

bs = basis_set(dag)
bs |> display

adjustment_sets(dag, :u, :w) |> display

#end


text:
#!/usr/bin/julia
#
# Structures to contain wave functions and collections of wave functions.
#
import FFTW
export WaveObject, WaveFunction, WaveMatrix, WaveBasis, WaveBasisMatrix, transform;

"""
    WaveObject{S<:Space, N, M, T<:Number, G<:Grid{S,M}} <: AbstractArray{T,N}

Abstract type from which implementations of discrete space representations of
objects related to quantum-mechanical wave functions are derived.

Grids should be defined on the same space and have the same element type as the
`WaveObject`.

`N` is the dimension of the `WaveObject` while `M` is the dimension of the grid.
"""
abstract type WaveObject{T<:Number, N, S, M, G<:Grid{S,M}} <: AbstractArray{T,N} end;
# Use the data array in the composite type implementation to act like a array.
Base.size(W::WaveObject) = size(W.data);
Base.getindex(W::WaveObject{T,N}, I::Vararg{Int,N}) where {T,N} =
  get(W.data, I, zero(T));
Base.setindex!(W::WaveObject{T,N}, v, I::Vararg{Int,N}) where {T,N} =
  W.data[I...] = v;
# Use linear indexing for 1D `WaveObjects`.
Base.IndexStyle(::Type{<:WaveObject{T,1}}) where T = IndexLinear();
Base.getindex(W::WaveObject{T,1}, I::Int) where T = W.data[I];
Base.setindex!(W::WaveObject{T,1}, v, I::Int) where T = (W.data[I] = v);

"Type alias for matrix-like collections of 1D `WaveFunctions`."
WaveMatrix{T,S,G} = WaveObject{T,2,S,1,G}

"""
`WaveFunction` objects are same dimension as underlying `Grid` object they are
defined upon.
"""
struct WaveFunction{T,N,S,G} <: WaveObject{T,N,S,N,G}
  data::Array{T,N};

  "Spacing between the elements in the discrete grid."
  Δ::Array;

  """
      WaveFunction{T,N,S,G}(data::AbstractArray{K,N}, Δ::AbstractArray)

  Construct `WaveFunction` object with the provided data.
  """
  WaveFunction{T,N,S,G}(data::AbstractArray{K,N}, Δ::AbstractArray) where {T,K,N,S,G} =
    new(data, Δ);

  """
      WaveFunction{T,1,S,G}(data::AbstractVector, Δ::Number) where G<:UniformGrid

  Construct one dimensional `WaveFunction` on a `UniformGrid` which is a special
  case where there is only a single value for `Δ`.
  """
  WaveFunction{T,1,S,G}(data::AbstractVector, Δ::Number) where {T,S,G<:UniformGrid} =
    new(data, [Δ]);

end

"""
    WaveFunction{T}(grid::G) where {T,G<:Grid{S,N}} =

Construct `WaveFunction` with zero values everywhere on the grid.
"""
WaveFunction{T}(grid::G) where {T,N,S,G<:Grid{S,N}} =
  WaveFunction{T,N,S,G}(zeros(T,size(grid)), grid.Δ);

"""
    WaveFunction(grid::Grid)

Construct `WaveFunction` with zero values of type `Complex{Float64}` everywhere.
"""
WaveFunction(grid::Grid) = WaveFunction{Complex{Float64}}(grid);

"""
    WaveFunction{T}(grid::G, ψ::AbstractArray{K,N}) where {T,K,N,S,G<:Grid{S,N}}

Construct wave function on the grid with values at the grid points given by ψ.
"""
function WaveFunction{T}(grid::G, ψ::AbstractArray{K,N}) where {T,K,N,S,G<:Grid{S,N}}
  check_same_size(grid, ψ);
  WaveFunction{T,N,S,G}(ψ, grid.Δ);
end
function WaveFunction{T}(grid::G, ψ::AbstractArray{K,N}) where {A,B,U,T<:Unitful.Quantity{A,B,U},K,N,S,G<:Grid{S,N}}
  check_same_size(grid, ψ);
  WaveFunction{T,N,S,G}(ustrip.(ψ)*unit(T), grid.Δ);
end
"""
    WaveFunction(grid::G, ψ::AbstractArray{T,N}) where {T,N,S,G<:Grid{S,N}}

Construct `WaveFunction` of the same type as the elements of `ψ`.
"""
WaveFunction(grid::G, ψ::AbstractArray{T,N}) where {T,N,S,G<:Grid{S,N}} =
  WaveFunction{T}(grid, ψ);

"""
    WaveFunction(ψ::WaveObject{T,N,S,M,G}) where {T,N,S,M,G}

Construct `WaveFunction` from a compatible `WaveObject`.

Works by dropping any empty dimensions until the `WaveObject` dimensions match
those of the `Grid` it is defined upon.
"""
function WaveFunction{TN}(ψ::WaveObject{T,N,S,M,G}) where {TN,T,N,S,M,G}
  Δ = ψ.Δ;
  for i ∈ N:-1:M+1
    ψ = dropdims(convert(Array{TN}, ψ), dims=i);
  end
  WaveFunction{TN,M,S,G}(ψ, Δ);
end
WaveFunction(ψ::WaveObject{T}) where T = WaveFunction{T}(ψ);

"""
    getproperty(ψ::WaveFunction{S,1,T,G}, s::Symbol) where {S,N,T,G<:UniformGrid}

Override `getproperty` so just one spacing is returned when `WaveFunction` is
defined on a one-dimensional `UniformGrid`.
"""
function Base.getproperty(ψ::WaveFunction{T,1,S,G}, s::Symbol) where
          {T,N,S,G<:UniformGrid}
  if s === :Δ return Base.getfield(ψ, s)[1];
  else return Base.getfield(ψ, s); end
end

"Generic basis of `WaveFunction` objects."
struct WaveBasis{T,N,S,M,G<:Grid{S,M}} <: WaveObject{T,N,S,M,G}
  data::Array{T,N};

  "Spacing between the elements in the discrete grid."
  Δ::Array;

  """
      WaveBasis{T,N,S,M,G}(data::AbstractArray{T,N}, Δ::AbstractArray)

  Construct `WaveBasis` ensuring that dimensionality is consistent with storing
  `WaveFunction` objects with `M` dimensions.
  """
  function WaveBasis{T,N,S,M,G}(data::AbstractArray{T,N}, Δ::AbstractArray) where
      {T,N,S,M,G}
    if N !== M+1
      throw(ArgumentError("WaveBasis must have exactly one more dimension than grid."));
    end
    new(data, Δ);
  end

  """
      WaveBasis{T,2,S,1,G}(data::AbstractMatrix{T}, Δ::Number) where G<:UniformGrid

  Construct `WaveBasis` which will contain one dimensional `WaveFunction`
  objects defined on a `UniformGrid`.

  For such objects the spacing, `Δ` is a single value.
  """
  WaveBasis{T,2,S,1,G}(data::AbstractMatrix{T}, Δ::Number) where {T,S,G<:UniformGrid} =
    new(data, [Δ]);
end

"Alias for `WaveBasis` of 1D `WaveFunction` objects"
WaveBasisMatrix{T,S,G} = WaveBasis{T,2,S,1,G}

"""
    WaveBasis{T}(grid::G, size::Int) where {T,S,M::Number,G<:Grid{S,M}}

Construct `WaveBasis` object with space for `size` basis elements with entries
of type `T`, on `grid`.
"""
WaveBasis{T}(grid::G, size::Int) where {T,S,M,G<:Grid{S,M}} =
  WaveBasis{T,M+1,S,M,G}(zeros(T,length(grid),size), grid.Δ);

"""
    WaveBasis(grid::Grid, size::Int)

Construct `WaveBasis` with space for `size` basis elements of type
`Complex{Float64}` on `grid`.
"""
WaveBasis(grid::Grid, size::Int) = WaveBasis{Complex{Float64}}(grid, size);

"""
    getproperty(Φ::WaveBasis{T,2,S,1,G}, s::Symbol) where {T,S,G<:UniformGrid}

Override `getproperty` so just one spacing is returned when `WaveBasis` stores
one-dimensional `WaveFunction` objects on a `UniformGrid`.
"""
function Base.getproperty(Φ::WaveBasis{T,2,S,1,G}, s::Symbol) where
         {T,S,G<:UniformGrid}
  if s === :Δ return Base.getfield(Φ, s)[1];
  else return Base.getfield(Φ, s); end
end

"""
    check_same_size(grid::Grid, v::AbstractArray)

Check the size of the grid matches that of the passed array.
"""
function check_same_size(grid::Grid, a::AbstractArray)
  if size(a) != size(grid)
    throw(DimensionMismatch("array must be same size as the grid."));
  end
end
function check_same_size(grid::Grid{S,1}, v::Vector) where S
  if length(grid) != length(v)
    throw(DimensionMismatch("vector must be same length as the grid."));
  end
end
function check_same_size(v::AbstractArray, grid::Grid)
  check_same_size(grid, v)
end
function check_same_size(v::AbstractVector, w::AbstractVector)
  if length(v) != length(w)
    throw(DimensionMismatch("abstract vectors must be same length."));
  end
end

"""
    transform(ψ::WaveFunction{T,N,S}, g::Grid{M}) where
        {T,N,S<:ConfigurationSpace,M<:MomentumSpace}

Construct WaveFunction from `ψ` transformed onto momentum space.
"""
function transform(ψ::WaveFunction{T,N,S}, g::Grid{M}) where
            {T,N,S<:ConfigurationSpace,M<:MomentumSpace}
  return WaveFunction(g, FFTW.fft(ψ));
end
function transform(ψ::WaveFunction{T,N,S}, g::Grid{M,R,Q}) where
    {T<:Unitful.Quantity,N,S<:ConfigurationSpace,M<:MomentumSpace,
     R,A,D,U,Q<:Unitful.Quantity{A,D,U}}
    # Determine quantity of transformed WaveFunction.
    newquant = typeof((1.0+1.0im)*unit(Q)^(-1/2));
    return WaveFunction{newquant}(g, FFTW.fft(ψ));
end
function transform(ψ::WaveFunction{T,N,S}, g::Grid{S}) where
            {T,N,S<:ConfigurationSpace}
  return transform(ψ, Grid(g));
end
"""
    transform(ψ::WaveFunction{T,N,S}, g::Grid{S}) where {T,N,S<:MomentumSpace}

Construct WaveFunction from `ψ` transformed onto configuration space.
"""
function transform(ψ::WaveFunction{T,N,S}, g::Grid{C}) where
            {T,N,S<:MomentumSpace,C<:ConfigurationSpace}
  return WaveFunction(g, FFTW.ifft(ψ));
end
function transform(ψ::WaveFunction{T,N,S}, g::Grid{C,R,Q}) where
    {T<:Unitful.Quantity,N,S<:MomentumSpace,C<:ConfigurationSpace,
     R,A,D,U,Q<:Unitful.Quantity{A,D,U}}
    # Determine quantity of transformed WaveFunction.
    newquant = typeof((1.0+1.0im)*unit(Q)^(-1/2));
    return WaveFunction{newquant}(g, FFTW.ifft(ψ));
end
"""
    FFTW.fft(ψ::WaveFunction{Complex{Float64}})

Convert complex valued `WaveFunction` to `Vector` objects prior to calling
`FFTW.fft`.

Workaround for a dispatch bug in FFTW.
"""
FFTW.fft(ψ::WaveFunction) = FFTW.fft(convert(Vector{Complex{Float64}}, ustrip.(ψ)));
"""
    FFTW.ifft(ψ::WaveFunction{Complex{Float64}})

Convert complex valued `WaveFunction` to `Vector` objects prior to calling
`FFTW.ifft`.

Workaround for a dispatch bug in FFTW.
"""
FFTW.ifft(ψ::WaveFunction) = FFTW.ifft(convert(Vector{Complex{Float64}}, ustrip.(ψ)));

#Base.:*(Φ::WaveMatrix{T,S,G}, ψ::WaveFunction{K,1,S,G}) where {T,K,S,G} =
#  WaveFunction{T,1,S,G}(Base.:*(Φ,convert(Vector,ψ)), ψ);

# Take the real part of a wavefunction.
Base.real(ψ::WaveFunction{Complex{T},N,S,G}) where {T,N,S,G} =
  WaveFunction{Complex{T},N,S,G}(real.(ψ), ψ.Δ);