text:
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal
"""
a "safe" version of `floatmin`, such that `1/sfmin` does not overflow.
"""
function safemin(T)
    sfmin = floatmin(T)
    small = one(T) / floatmax(T)
    if small >= sfmin
        sfmin = small * (one(T) + eps(T))
    end
    sfmin
end

function _scale!(A::AbstractArray{T}) where {T}
    smlnum = sqrt(safemin(real(T))) / eps(real(T))
    bignum = 1 / smlnum
    anrm = norm(A,Inf)
    scaleA = false
    cscale = one(real(T))
    if anrm > 0 && anrm < smlnum
        scaleA = true
        cscale = smlnum
    elseif anrm > bignum
        scaleA = true
        cscale = bignum
    end
    scaleA && safescale!(A, anrm, cscale)
    scaleA, cscale, anrm
end

abs1(z) = abs(z)
abs1(z::T) where {T <: Complex} = abs(real(z))+abs(imag(z))

# translated from xLASCL (LAPACK)
"""
multiply an array `A` by a real number `cto/cfrom` expressed as a ratio,
computing without over/underflow where possible.

`cfrom` must not be zero.
"""
function safescale!(A::AbstractArray{T}, cfrom::Real, cto::Real) where T
    # CHECKME: what are correct type constraints on cfrom,cto?
    isnan(cto) && throw(ArgumentError("cto must be a valid number"))
    (isnan(cfrom) || (cfrom == 0)) && throw(ArgumentError("cfrom must be a valid nonzero number"))
    smlnum = safemin(real(T))
    bignum = 1/smlnum
    cfromc = cfrom
    ctoc = cto
    done = false
    while !done
        cfrom1 = cfromc * smlnum
        if cfrom1 == cfromc
            # cfromc is ±∞; multiplier is signed 0 or NaN
            mul = ctoc / cfromc
            done = true
            cto1 = ctoc
        else
            cto1 = ctoc / bignum
            if cto1 == ctoc
                # ctoc is 0 or ±∞, so is correct multiplier
                mul = cto
                done = true
                cfromc = one(real(T))
            elseif abs(cfrom1) > abs(ctoc) && (ctoc != 0)
                mul = smlnum
                done = false
                cfromc = one(real(T))
            elseif abs(cto1) > abs(cfromc)
                mul = bignum
                done = false
                cfromc = cfrom1
            else
                mul = ctoc / cfromc
                done = true
            end
        end
        A .= A .* mul
    end
end

# Rank one update

## General
### BLAS
rankUpdate!(α::T, x::StridedVector{T}, y::StridedVector{T}, A::StridedMatrix{T}) where {T<:BlasReal} = BLAS.ger!(α, x, y, A)

### Generic
function rankUpdate!(α::Number, x::StridedVector, y::StridedVector, A::StridedMatrix)
    m, n = size(A, 1), size(A, 2)
    m == length(x) || throw(DimensionMismatch("x vector has wrong length"))
    n == length(y) || throw(DimensionMismatch("y vector has wrong length"))
    for j = 1:n
        yjc = y[j]'
        for i = 1:m
            A[i,j] += x[i]*α*yjc
        end
    end
end

# index-unsafe upper triangular solve
# with scaling to prevent overflow
# Actually solves A[1:n,1:n] * x = s * b and returns s, a scaling factor
# such that components of x are manageable.
# NB: works with leading n×n block of A and leading n terms of b, x
# Hence all the loops must be explicit.
# This is a translation of (part of) LAPACK::zlatrs, q.v. for details.
# Provide `b` in argument `x`; note that `cnorm` may be updated.
function _usolve!(A::StridedMatrix{T}, n, x, cnorm) where {T}

    cabs1half(z) = abs(real(z)/2) + abs(imag(z)/2)

    RT = real(T)
    half = one(RT) / 2
    smallnum = safemin(RT) / eps(RT)
    bignum = one(RT)/ smallnum
    xscale = one(RT)

    tmax = abs(cnorm[1])
    @inbounds for j=2:n
        tmax = max(tmax, abs(cnorm[j]))
    end

    if tmax <= bignum * half
        tscale = one(RT)
    else
        tscale = half / (smallnum * tmax)
        cnorm[1:n] .*= tscale
    end
    # bound on solution vector:
    xmax = cabs1half(x[1])
    @inbounds for j=2:n
        xmax = max(xmax,cabs1half(x[j]))
    end
    xbound = xmax
    if tscale != one(RT)
        grow = zero(RT)
    else
        # compute grow
        grow = half / max(xbound, smallnum)
        xbound = grow
        toosmall = false
        for j=n:-1:1
            # give up if too small
            toosmall = (grow <= smallnum)
            toosmall && break
            tjjs = A[j,j]
            tjj = abs1(tjjs)
            if tjj >= smallnum
                xbound = min(xbound, min(one(RT), tjj)*grow)
            else
                # M[j] could overflow
                xbound = zero(RT)
            end
            if tjj + cnorm[j] >= smallnum
                grow *= (tjj / (tjj + cnorm[j]))
            else
                # grow could overflow, clamp
                grow = zero(RT)
            end
        end
        if !toosmall
            grow = xbound
        end
    end

    if grow * tscale > smallnum
        # bound is ok; use standard arithmetic
        @inbounds for j in n:-1:1
            xj = x[j] = A[j,j] \ x[j]
            for i in j-1:-1:1
                x[i] -= A[i,j] * xj
            end
        end
    else
        if xmax > bignum * half
            # scale so all(abs.(x) .<= bignum)
            xscale = (bignum * half) / xmax
            @inbounds for i=1:n; x[i] *= xscale; end
            xmax = bignum
        else
            xmax *= 2
        end
        for j in n:-1:1
            # compute x[j] = b[j] / A[j,j], scaling if necessary
            xj = abs1(x[j])
            tjjs = A[j,j] * tscale
            tjj = abs1(tjjs)
            if tjj > smallnum
                # abs(A[j,j] > smallnum):
                if tjj < one(RT)
                    if xj > tjj*bignum
                        # scale x by 1/b[j]
                        rec = one(RT) / xj
                        @inbounds for i=1:n; x[i] *= rec; end
                        xscale *= rec
                        xmax *= rec
                    end
                end
                x[j] = x[j] / tjjs
                xj = abs1(x[j])
            elseif tjj > 0
                # 0 < abs(A[j,j]) <= smallnum
                if xj > tjj*bignum
                    # scale by (1/abs(x[j]))*abs(A[j,j]))*bignum
                    # to avoid overflow when dividing by A[j,j]
                    rec = (tjj * bignum) / xj
                    if cnorm[j] > one(RT)
                        # scale by 1/cnorm[j] to avoid overflow
                        # when multiplying x[j] by column j
                        rec /= cnorm[j]
                    end
                    @inbounds for i=1:n; x[i] *= rec; end
                    xscale *= rec
                    xmax *= rec
                end
                x[j] = x[j] / tjjs # FIXME: use zladiv scheme
                xj = abs1(x[j])
            else
                # A[j,j] = 0; set x to eⱼ, xscale to 0
                # and compute a null vector.
                @inbounds for i=1:n; x[i] = zero(T); end
                x[j] = one(T)
                xscale = zero(RT)
                xmax = zero(RT)
            end
            # scale x if necc. to avoid overflow when adding multiple of A[:,j]
            if xj > one(RT)
                rec = one(RT) / xj
                if cnorm[j] > (bignum - xmax) * rec
                    # scale by 1/(2abs(x[j]))
                    rec *= half
                    @inbounds for i=1:n; x[i] *= rec; end
                    xscale *= rec
                end
            elseif (xj * cnorm[j] > bignum - xmax)
                # scale by 1/2
                @inbounds for i=1:n; x[i] *= half; end
                xscale *= half
            end
            if j > 1
                # compute update
                xjt = x[j] * tscale
                @inbounds for i=1:j-1
                    x[i] -= xjt * A[i,j]
                end
                xmax = abs1(x[1])
                @inbounds for i=2:j-1
                    xmax = max(xmax, abs1(x[i]))
                end
            end
        end
        xscale /= tscale
    end
    if tscale != one(RT)
        @inbounds for i=1:n; cnorm[i] *= one(RT)/tscale; end
    end
    return xscale
end

# conjugate transpose version
function _cusolve!(A::StridedMatrix{T}, n, x, cnorm) where {T}

    cabs1half(z) = abs(real(z)/2) + abs(imag(z)/2)

    RT = real(T)
    half = one(RT) / 2
#    smallnum = safemin(RT) / (2*eps(RT)) # WARNING: assumes base 2
    smallnum = safemin(RT) / eps(RT)
    bignum = one(RT)/ smallnum
    xscale = one(RT)

    tmax = abs(cnorm[1])
    @inbounds for j=2:n
        tmax = max(tmax, abs(cnorm[j]))
    end

    if tmax <= bignum * half
        tscale = one(RT)
    else
        tscale = half / (smallnum * tmax)
        cnorm[1:n] .*= tscale
    end
    # bound on solution vector:
    xmax = cabs1half(x[1])
    @inbounds for j=2:n
        xmax = max(xmax,cabs1half(x[j]))
    end
    xbound = xmax
    if tscale != one(RT)
        grow = zero(RT)
    else
        # compute grow
        grow = half / max(xbound, smallnum)
        xbound = grow
        toosmall = false
        for j=1:n
            # give up if too small
            toosmall = (grow <= smallnum)
            toosmall && break
            # G[j] = max(G[j-1], M[j-1]*(1+cnorm[j]))
            xj = one(RT) + cnorm[j]
            grow = min(grow, xbound / xj)
            tjjs = A[j,j]
            tjj = abs1(tjjs)
            if tjj >= smallnum
                if xj > tjj
                    xbound *= (tjj / xj)
                end
            else
                # M[j] could overflow
                xbound = zero(RT)
            end
        end
        if !toosmall
            grow = min(grow, xbound)
        end
    end

    if grow * tscale > smallnum
        # bound is ok; use standard arithmetic
        @inbounds for j in 1:n
            z = x[j]
            for i in 1:j-1
                z -= conj(A[i,j]) * x[i]
            end
            x[j] = conj(A[j,j]) \ z
        end
    else
        if xmax > bignum * half
            # scale so all(abs.(x) .<= bignum)
            xscale = (bignum * half) / xmax
            @inbounds for i=1:n; x[i] *= xscale; end
            xmax = bignum
        else
            xmax *= 2
        end
        for j=1:n
            # compute x[j] = b[j] - Σ_(k≠j) A[k,j] x[k]
            xj = abs1(x[j])
            uscale = complex(tscale)
            rec = one(RT) / max(xmax, one(RT))
            if cnorm[j] > (bignum - xj) * rec
                # if x[j] could overflow scale x by 1/(2 xmax)
                rec *= half
                tjjs = conj(A[j,j]) * tscale
                tjj = abs1(tjjs)
                if tjj > one(RT)
                    # divide by A[j,j] when scaling if A[j,j] > 1
                    rec = min(one(RT), rec*tjj)
                    uscale /= tjjs
                end
                if rec < one(RT)
                    @inbounds for i=1:n; x[i] *= rec; end
                    xscale *= rec
                    xmax *= rec
                end
            end
            csumj = zero(T)
            @inbounds for i=1:j-1
                csumj += (conj(A[i,j]) * uscale) * x[i]
            end
            if uscale == complex(tscale)
                # if diagonal wasn't used to scale
                # compute x[j] = (x[j] - csumj) / A[j,j]
                x[j] -= csumj
                xj = abs1(x[j])
                tjjs = conj(A[j,j]) * tscale

                # compute x[j] /= A[j,j], scaling if necessary
                tjj = abs1(tjjs)
                if tjj > smallnum
                    if tjj < one(RT)
                        if xj > tjj * bignum
                            rec = one(RT) / xj
                            @inbounds for i=1:n; x[i] *= rec; end
                            xscale *= rec
                            xmax *= rec
                        end
                    end
                    x[j] = x[j] / tjjs
                elseif tjj > zero(RT)
                    # tiny diag
                    if xj > tjj * bignum
                        rec = (tjj * bignum) / xj
                        @inbounds for i=1:n; x[i] *= rec; end
                        xscale *= rec
                        xmax *= rec
                    end
                    x[j] = x[j] / tjjs
                else
                    # zero diag: compute a null vector of Aᴴ
                    x[1:n] .= zero(T)
                    x[j] = one(T)
                    xscale = zero(RT)
                    xmax = zero(RT)
                end
            else
                # if dot product was already divided by A[j,j]
                # compute x[j] = x[j] / A[j,j] - csumj
                x[j] = x[j] / tjjs - csumj
            end
            xmax = max(xmax, abs1(x[j]))
        end
        xscale /= tscale
    end
    if tscale != one(RT)
        @inbounds for i=1:n; cnorm[i] *= one(RT)/tscale; end
    end
    xscale
end

# update scale and sumsq s.t. (scale^2 * sumsq) is increased by ∑ |x[j]|^2
# scale is a running quasi-inf-norm.
# based on LAPACK::zlassq
function _ssq(x::AbstractVector{T}, scale, sumsq) where {T}
    n = length(x)
    rone = one(real(T))
    for ix=1:n
        t1 = abs(real(x[ix]))
        if t1 > 0 || isnan(t1)
            if scale < t1
                sumsq = rone + sumsq * (scale / t1)^2
                scale = t1
            else
                sumsq += (t1 / scale)^2
            end
        end
        if T <: Complex
            t1 = abs(imag(x[ix]))
            if t1 > 0 || isnan(t1)
                if scale < t1
                    sumsq = rone + sumsq * (scale / t1)^2
                    scale = t1
                else
                    sumsq += (t1 / scale)^2
                end
            end
        end
    end
    return scale, sumsq
end

# stdlib norm2 uses a poorly implemented generic scheme for short vectors.
function _norm2(x::AbstractVector{T}) where {T<:Real}
    require_one_based_indexing(x)
    n = length(x)
    n < 1 && return zero(T)
    n == 1 && return abs(x[1])
    scale = zero(T)
    ssq = zero(T)
    for xi in x
        if !iszero(xi)
            a = abs(xi)
            if scale < a
                ssq = one(T) + ssq * (scale / a)^2
                scale = a
            else
                ssq += (a / scale)^2
            end
        end
    end
    return scale * sqrt(ssq)
end

function _norm2(x::AbstractVector{T}) where {T<:Complex}
    require_one_based_indexing(x)
    n = length(x)
    RT = real(T)
    n < 1 && return zero(RT)
    n == 1 && return abs(x[1])
    scale = zero(RT)
    ssq = zero(RT)
    for xx in x
        xr,xi = reim(xx)
        if !iszero(xr)
            a = abs(xr)
            if scale < a
                ssq = one(RT) + ssq * (scale / a)^2
                scale = a
            else
                ssq += (a / scale)^2
            end
        end
        if !iszero(xi)
            a = abs(xi)
            if scale < a
                ssq = one(RT) + ssq * (scale / a)^2
                scale = a
            else
                ssq += (a / scale)^2
            end
        end
    end
    return scale * sqrt(ssq)
end

# As of v1.5, Julia hypot() w/ >2 args is unprotected (the documentation lies),
# so we need this.
# translation of dlapy3, assuming NaN propagation
function _hypot3(x::T, y::T, z::T) where {T}
    xa = abs(x)
    ya = abs(y)
    za = abs(z)
    w = max(xa, ya, za)
    rw = one(real(T)) / w
    r::real(T) = w * sqrt((rw * xa)^2 + (rw * ya)^2 + (rw * za)^2)
    return r
end

function _enormalize!(v::AbstractMatrix{T}) where T <: STypes
    n = size(v,1)
    for j=1:n
        s = one(real(T)) / norm(v[:,j],2)
        t = abs2(v[1,j])
        i0 = 1
        for i=2:n
            u = abs2(v[i,j])
            if  u > t
                i0 = i
                t = u
            end
        end
        t = s * conj(v[i0,j]) / sqrt(t)
        for i=1:n
            v[i,j] *=  t
        end
        v[i0,j] = real(v[i0,j])
    end
end


text:
using Symata, Test

println()
@sym begin
  xi = [1, 3, 6]
  yi = [1, 5, 10]
  L0(x_) := Simplify( Expand( ((x-xi[2])*(x-xi[3]))/((xi[1]-xi[2])*(xi[1]-xi[3])) ) )
  L1(x_) := Simplify( Expand( ((x-xi[1])*(x-xi[3]))/((xi[2]-xi[1])*(xi[2]-xi[3])) ) )
  L2(x_) := Simplify( Expand( ((x-xi[1])*(x-xi[2]))/((xi[3]-xi[1])*(xi[3]-xi[2])) ) )
  Q(x_) := Simplify(yi[1]*L0(x) + yi[2]*L1(x) + yi[3]*L2(x));
  SetJ(q, ToString(Q(x)))
end

@sym Print("Q(x) = ", Q(x))
println()
@sym Println("Q(4.5) = ", Q(4.5))

println()
q = @sym ToString(Q(x));
q |> display
println()

@eval f1(x) = $(Meta.parse(q))

@test round(f1(4.5), digits=5) == 7.65

#=
using GR

xint = 1:0.1:6 # FloatRange
p = plot(xint, f1.(xint),"g", [4.5], [f1(4.5)], "o", title="Interpolated curve")
=#

ext: jl


text:
using CrypticServer, Sockets

CrypticServer.main(port=8080, host=ip"0.0.0.0")


text:
module fiber_reinf_cant_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
# using IterativeSolvers
using Statistics: mean
using LinearAlgebra: Symmetric, cholesky, norm

function fiber_reinf_cant_iso()
    println("""
    Cantilever example.  Isotropic material.
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    t0 = time()
    # # Orthotropic material parameters of the external cylinder
    # E1s = 130.0*phun("GPa")
    # E2s = 5.0*phun("GPa")
    # E3s = E2s
    # nu12s = nu13s = 0.25
    # nu23s = 0.0
    # G12s = 10.0*phun("GPa")
    # G13s = G12s
    # G23s = 5.0*phun("GPa")
    # CTE1 = 3.0e-6
    # CTE2 = 2.0e-5
    # CTE3 = 2.0e-5
    # Isotropic material
    E = 1.0e9*phun("Pa")
    nu = 0.25
    CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref =  -7.516310912734678e-06;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    # Generate mesh
    n = 16
    na = n # number of elements lengthwise
    nb = n # number of elements through the wwith
    nt = n # number of elements through the thickness
    xs = collect(linearspace(0.0, a, na+1))
    ys = collect(linearspace(0.0, b, nb+1))
    ts = collect(linearspace(0.0, t, nt+1))
    fens,fes = H8blockx(xs, ys, ts)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

    MR = DeforModelRed3D
    # externalmaterial = MatDeforElastOrtho(MR,
    #   0.0, E1s, E2s, E3s,
    #   nu12s, nu13s, nu23s,
    #   G12s, G13s, G23s,
    #   CTE1, CTE2, CTE3)
    material = MatDeforElastIso(MR,
    0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = GaussRule(3, 2)

    region = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, gr), material))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

    ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
    ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

    function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(forceout, q0*[0.0; 0.0; 1.0])
    end

    Trac = FDataDict("traction_vector"=>getshr!, "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region],
    "essential_bcs"=>[ex01, ex02, ex03],
    "traction_bcs"=>[Trac],
    "temperature_change"=>FDataDict("temperature"=>dT)
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println(" Normalized deflection: $(utip/uz_ref)")
    println("Solution: $(  time()-t0 )")

    # File =  "NAFEMS-R0031-2-plate.vtk"
    # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20;
    #     scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    # @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso",
    "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>5)
    # modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

    println("Done: $(  time()-t0 )")
    true

end # fiber_reinf_cant_iso

function fiber_reinf_cant_iso_stresses()
    println("""
    Cantilever example.  Isotropic material.
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    # Isotropic material
    E = 1.0e9*phun("Pa")
    nu = 0.25
    CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref =  -7.516310912734678e-06;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    # Generate mesh
    for n = [2 4 8 16 32]

        na = n # number of elements lengthwise
        nb = n # number of elements through the wwith
        nt = n # number of elements through the thickness
        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        ts = collect(linearspace(0.0, t, nt+1))
        fens,fes = H8blockx(xs, ys, ts)
        # fens,fes = H8toH20(fens,fes)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

        MR = DeforModelRed3D
        # externalmaterial = MatDeforElastOrtho(MR,
        #   0.0, E1s, E2s, E3s,
        #   nu12s, nu13s, nu23s,
        #   G12s, G13s, G23s,
        #   CTE1, CTE2, CTE3)
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = GaussRule(3, 2)

        region = FDataDict("femm"=>FEMMDeforLinearMSH8(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

        function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(forceout, q0*[0.0; 0.0; 1.0])
        end

        Trac = FDataDict("traction_vector"=>getshr!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), GaussRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03],
        "traction_bcs"=>[Trac],
        "temperature_change"=>FDataDict("temperature"=>dT)
        )
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println(" Normalized deflection: $(utip/uz_ref)")

        stressfields = NodalField[]
        for c = 1:6
            modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso",
            "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>c)
            modeldata = AlgoDeforLinearModule.exportstress(modeldata)
            push!(stressfields, modeldata["postprocessing"]["exported"][1]["field"])
        end

        # println("@__FILE__ = $(@__FILE__)")
        # jldopen("fiber_reinf_cant_iso_stresses" * "n=$(n)" * ".jld", "w") do file
        #     write(file, "n", n)
        #     write(file, "fens", fens)
        #     write(file, "fes", fes)
        #     write(file, "geom", geom)
        #     write(file, "u", u)
        #     # write(file, "femm", region["femm"])
        #     write(file, "integrationrule", region["femm"].integdomain.integration_rule)
        #     write(file, "stressfields", stressfields)
        #     write(file, "tolerance", tolerance)
        # end

    end
    true

end # fiber_reinf_cant_iso_stresses

function fiber_reinf_cant_iso_stresses_MST10()
    elementtag = "MST10"
    println("""
    Cantilever example.  Isotropic material.
    ############### $(elementtag) ###############
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    # Isotropic material
    E = 1.0e9*phun("Pa")
    nu = 0.25
    CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref =  -7.516310912734678e-06;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    convergencestudy = FDataDict[]
    for n = [1 2 4 8 16]

        na = n # number of elements lengthwise
        nb = n # number of elements through the wwith
        nt = n # number of elements through the thickness
        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        ts = collect(linearspace(0.0, t, nt+1))
        fens,fes = T10blockx(xs, ys, ts)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

        MR = DeforModelRed3D
        # externalmaterial = MatDeforElastOrtho(MR,
        #   0.0, E1s, E2s, E3s,
        #   nu12s, nu13s, nu23s,
        #   G12s, G13s, G23s,
        #   CTE1, CTE2, CTE3)
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

        function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(forceout, q0*[0.0; 0.0; 1.0])
        end

        Trac = FDataDict("traction_vector"=>getshr!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03],
        "traction_bcs"=>[Trac],
        "temperature_change"=>FDataDict("temperature"=>dT)
        )
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println(" Normalized deflection: $(utip/uz_ref)")

        modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / n,
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>region["femm"],
        "integrationrule"=>region["femm"].integdomain.integration_rule,
        "stressfields"=>stressfields,
        "tolerance"=>tolerance)
        )
    end

    # File = "fiber_reinf_cant_iso_stresses_$(elementtag)"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

    true

end # fiber_reinf_cant_iso_stresses_MST10

function fiber_reinf_cant_iso_stresses_T10()
    elementtag = "T10"
    println("""
    Cantilever example.  Isotropic material.
    ############### $(elementtag) ###############
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    # Isotropic material
    E = 1.0e9*phun("Pa")
    nu = 0.25
    CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref =  -7.516310912734678e-06;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    convergencestudy = FDataDict[]
    for n = [1 2 4 8 16]

        na = n # number of elements lengthwise
        nb = n # number of elements through the wwith
        nt = n # number of elements through the thickness
        xs = collect(linearspace(0.0, a, na+1))
        ys = collect(linearspace(0.0, b, nb+1))
        ts = collect(linearspace(0.0, t, nt+1))
        fens,fes = T10blockx(xs, ys, ts)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

        MR = DeforModelRed3D
        # externalmaterial = MatDeforElastOrtho(MR,
        #   0.0, E1s, E2s, E3s,
        #   nu12s, nu13s, nu23s,
        #   G12s, G13s, G23s,
        #   CTE1, CTE2, CTE3)
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

        function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(forceout, q0*[0.0; 0.0; 1.0])
        end

        Trac = FDataDict("traction_vector"=>getshr!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03],
        "traction_bcs"=>[Trac],
        "temperature_change"=>FDataDict("temperature"=>dT)
        )
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println(" Normalized deflection: $(utip/uz_ref)")

        modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_iso_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / n,
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>region["femm"],
        "integrationrule"=>region["femm"].integdomain.integration_rule,
        "stressfields"=>stressfields,
        "tolerance"=>tolerance)
        )
    end

    # File = "fiber_reinf_cant_iso_stresses_$(elementtag)"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

    true

end # fiber_reinf_cant_iso_stresses_T10

function fiber_reinf_cant_yn_strong()
    println("""
    Cantilever example.  Strongly orthotropic material. Orientation "y".
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    t0 = time()
    # # Orthotropic material
    E1s = 100000.0*phun("GPa")
    E2s = 1.0*phun("GPa")
    E3s = E2s
    nu23s = nu12s = nu13s = 0.25
    G12s = 0.2*phun("GPa")
    G23s = G13s = G12s
    CTE1 = 0.0
    CTE2 = 0.0
    CTE3 = 0.0
    # # Isotropic material
    # E = 1.0e9*phun("Pa")
    # nu = 0.25
    # CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref = -1.027498445054843e-05;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    # Generate mesh
    n = 4
    na = 8*n # number of elements lengthwise
    nb = n # number of elements through the wwith
    nt = n # number of elements through the thickness
    xs = collect(linearspace(0.0, a, na+1))
    ys = collect(linearspace(0.0, b, nb+1))
    ts = collect(linearspace(0.0, t, nt+1))
    fens,fes = H8blockx(xs, ys, ts)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    # material = MatDeforElastIso(MR,
    #   0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = zeros(3, 3)
    rotmat3!(csmat, -45.0/180.0*pi*[0,1,0])

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = GaussRule(3, 2)

    region = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, gr), CSys(3, 3, updatecs!), material))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

    ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
    ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

    function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(forceout, q0*[0.0; 0.0; 1.0])
    end

    Trac = FDataDict("traction_vector"=>getshr!, "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region],
    "essential_bcs"=>[ex01, ex02, ex03],
    "traction_bcs"=>[Trac],
    "temperature_change"=>FDataDict("temperature"=>dT)
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection $utip, normalized: $(utip/uz_ref)")
    println("Solution: $(  time()-t0 )")

    # File =  "NAFEMS-R0031-2-plate.vtk"
    # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20;
    #     scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    # @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_yn_strong", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>5)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
    # File = modeldata["postprocessing"]["exported"][1]["file"]
    # @async run(`"paraview.exe" $File`)

    println("Done: $(  time()-t0 )")
    true

end # fiber_reinf_cant_yn_strong

function fiber_reinf_cant_yn_strong_no_algo()
    println("""
    Cantilever example.  Strongly orthotropic material. Orientation "y".
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    t0 = time()
    pu = ustring -> phun(ustring; system_of_units = :SIMM)

    # # Orthotropic material
    E1s = 100000.0*pu("GPa")
    E2s = 1.0*pu("GPa")
    E3s = E2s
    nu23s = nu12s = nu13s = 0.25
    G12s = 0.2*pu("GPa")
    G23s = G13s = G12s
    CTE1 = 0.0
    CTE2 = 0.0
    CTE3 = 0.0
    # # Isotropic material
    # E = 1.0e9*pu("Pa")
    # nu = 0.25
    # CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref = -1.027498445054843e-05*pu("m");

    a = 90.0*pu("mm") # length of the cantilever
    b = 10.0*pu("mm") # width of the cross-section
    t = 20.0*pu("mm") # height of the cross-section
    q0 = -1000.0*pu("Pa") # shear traction
    dT = 0*pu("K") # temperature rise

    tolerance = 0.00001*t

    # Generate mesh
    n = 10
    na = n # number of elements lengthwise
    nb = n # number of elements through the wwith
    nt = n # number of elements through the thickness
    xs = collect(linearspace(0.0, a, na+1))
    ys = collect(linearspace(0.0, b, nb+1))
    ts = collect(linearspace(0.0, t, nt+1))
    println("fens,fes = H8blockx(xs, ys, ts)")
    @time fens,fes = H8blockx(xs, ys, ts)
    println("fens,fes = H8toH20(fens,fes)")
    @time fens,fes = H8toH20(fens,fes)
    println("bfes = meshboundary(fes)")
    @time bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    # material = MatDeforElastIso(MR,
    #   0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = zeros(3, 3)
    rotmat3!(csmat, -45.0/180.0*pi*[0,1,0])

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = GaussRule(3, 2)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, gr), CSys(3, 3, updatecs!), material)

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
    nnodes(geom)

    setebc!(u, lx0, true, 1, zeros(size(lx0)))
    setebc!(u, lx0, true, 2, zeros(size(lx0)))
    setebc!(u, lx0, true, 3, zeros(size(lx0)))
    applyebc!(u)

    S = connectionmatrix(femm, nnodes(geom))

    numberdofs!(u)

    function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(forceout, q0*[0.0; 0.0; 1.0])
    end

    Tracfemm = FEMMBase(IntegDomain(subset(bfes, sshearl), GaussRule(2, 3)))

    println("K = stiffness(femm, geom, u)")
    @time K = stiffness(femm, geom, u)
    fi = ForceIntensity(FFlt, 3, getshr!);
    println("F =  distribloads(Tracfemm, geom, u, fi, 2);")
    @time F =  distribloads(Tracfemm, geom, u, fi, 2);

    println("K = cholesky(K)")
    K = (K + K')/2;
    @time K = cholesky(Symmetric(K))
    println("U = K\\F")
    @time U = K\F
    # println("U = cg(K, F; tol=1e-3, maxiter=2000)")
    # @time U = cg(K, F; tol=1e-3, maxiter=2000)
    scattersysvec!(u, U[:])

    Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection $utip, normalized: $(utip/uz_ref)")
    println("Solution: $(  time()-t0 )")

    println("Done: $(  time()-t0 )")
    true

end # fiber_reinf_cant_yn_strong_no_algo

function fiber_reinf_cant_zn_strong()
    println("""
    Cantilever example.  Strongly orthotropic material. Orientation "z".
    @article{
    author = {Krysl, P.},
    title = {Mean-strain 8-node hexahedron with optimized energy-sampling stabilization},
    journal = {Finite Elements in Analysis and Design},
    volume = {108}, pages = {41-53}, DOI = {10.1016/j.finel.2015.09.008}, year = {2016}
    }
    """)

    t0 = time()
    # # Orthotropic material
    E1s = 100000.0*phun("GPa")
    E2s = 1.0*phun("GPa")
    E3s = E2s
    nu23s = nu12s = nu13s = 0.25
    G12s = 0.2*phun("GPa")
    G23s = G13s = G12s
    CTE1 = 0.0
    CTE2 = 0.0
    CTE3 = 0.0
    # # Isotropic material
    # E = 1.0e9*phun("Pa")
    # nu = 0.25
    # CTE = 0.0

    # Reference value for  the vertical deflection of the tip
    uz_ref = -1.119145781010554e-05;

    a = 90.0*phun("mm") # length of the cantilever
    b = 10.0*phun("mm") # width of the cross-section
    t = 20.0*phun("mm") # height of the cross-section
    q0 = -1000.0*phun("Pa") # shear traction
    dT = 0*phun("K") # temperature rise

    tolerance = 0.00001*t

    # Generate mesh
    n = 8
    na = 8*n # number of elements lengthwise
    nb = n # number of elements through the wwith
    nt = n # number of elements through the thickness
    xs = collect(linearspace(0.0, a, na+1))
    ys = collect(linearspace(0.0, b, nb+1))
    ts = collect(linearspace(0.0, t, nt+1))
    fens,fes = H8blockx(xs, ys, ts)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearl = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    # material = MatDeforElastIso(MR,
    #   0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = zeros(3, 3)
    rotmat3!(csmat, -45.0/180.0*pi*[0,0,1])

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = GaussRule(3, 2)

    region = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, gr), CSys(3, 3, updatecs!), material))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)

    ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
    ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )

    function getshr!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(forceout, q0*[0.0; 0.0; 1.0])
    end

    Trac = FDataDict("traction_vector"=>getshr!, "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region],
    "essential_bcs"=>[ex01, ex02, ex03],
    "traction_bcs"=>[Trac],
    "temperature_change"=>FDataDict("temperature"=>dT)
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    Tipl = selectnode(fens, box=[a a b b 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection $utip, normalized: $(utip/uz_ref)")
    println("Solution: $(  time()-t0 )")

    # File =  "NAFEMS-R0031-2-plate.vtk"
    # vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20;
    #     scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    # @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"fiber_reinf_cant_yn_strong", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>5)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    # modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
    # File = modeldata["postprocessing"]["exported"][1]["file"]
    # @async run(`"paraview.exe" $File`)

    println("Done: $(  time()-t0 )")
    true

end # fiber_reinf_cant_zn_strong

function allrun()
    println("#####################################################")
    println("# fiber_reinf_cant_iso ")
    fiber_reinf_cant_iso()
    println("#####################################################")
    println("# fiber_reinf_cant_iso_stresses ")
    fiber_reinf_cant_iso_stresses()
    println("#####################################################")
    println("# fiber_reinf_cant_iso_stresses_MST10 ")
    fiber_reinf_cant_iso_stresses_MST10()
    println("#####################################################")
    println("# fiber_reinf_cant_iso_stresses_T10 ")
    fiber_reinf_cant_iso_stresses_T10()
    println("#####################################################")
    println("# fiber_reinf_cant_yn_strong ")
    fiber_reinf_cant_yn_strong()
    println("#####################################################")
    println("# fiber_reinf_cant_yn_strong_no_algo ")
    fiber_reinf_cant_yn_strong_no_algo()
    return true
end # function allrun

end # module fiber_reinf_cant_examples


text:
using DataDeps

register(DataDep(
  "ozone-level-detection-eight",
  "https://archive.ics.uci.edu/ml/datasets/Ozone+Level+Detection",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/ozone/eighthr.data",
  "9a24c3e5c43570148869fec02edcd44dd1cb06e97c2bbc0b6ab24ebdd21de37e",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=74,
        feature_indices=2:73,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Pkg
Pkg.activate(".")
using Printf
push!(LOAD_PATH,"../src/")
using QuantumComputer

# configurable runner is below this function
function grover(qubit_count, target_input, sample_size)
    @printf("searching %d inputs for target of %d using Grover's algorithm...\n", 2^qubit_count, target_input)

    println("building circuit...")
    classical_register = QuantumComputer.ClassicalRegister(qubit_count)
    register = QuantumComputer.Register(qubit_count + 1, 1)
    superposition = QuantumComputer.Superposition(register.qubits)

    circuit_oracle = QuantumComputer.Circuits.identification_oracle(qubit_count + 1, target_input)
    gate_oracle = QuantumComputer.circuit_convert_to_gate(circuit_oracle)
    circuit_grover = QuantumComputer.Circuits.grover(gate_oracle, false)
    measurement = QuantumComputer.Measurement(Array(1:qubit_count), Array(1:qubit_count), sample_size)

    circuit = QuantumComputer.Circuit()
    QuantumComputer.add_subcircuit_to_circuit!(circuit, circuit_grover)
    QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

    @printf("sampling (%d times)...\n", sample_size)
    QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

    @printf("classical_register holds value %d\n", classical_register.value)
    classical_register.value == target_input || throw(ErrorException("Grover has failed!"))
    println("success!")
end

# this one is probably fine.
sample_size = 128

qubit_count = length(ARGS) > 0 ? parse(Int64, ARGS[1]) : 10
if length(ARGS) > 1
    QuantumComputer.gate_set_cache_base(ARGS[2])
end
target_input = length(ARGS) > 2 ? parse(Int64, ARGS[3]) : rand(0:(2^qubit_count-1))

grover(qubit_count, target_input, sample_size)


text:
TCA_condition(k, ℋₓ, τₓ′) = (abs(k / (ℋₓ * τₓ′)) < 0.1) & (abs(τₓ′) > 10.0)

# NOTE: NO NEUTRINOS 𝒩
function hierarchy!(du, u, p::AbstractCosmoParams, x)
    ℓᵧ, Ω_r, Ω_b, Ω_m = p.ℓᵧ, p.Ω_r, p.Ω_b, p.Ω_m
    H₀² = bg.H₀^2
    ℋₓ, ℋₓ′ = bg.ℋ(x), bg.ℋ′(x)
    τₓ′, τₓ′′ = ih.τ′(x), ih.τ′′(x)
    a = x2a(x)
    R = 4Ω_r / (3Ω_b * a)

    # get array views of photon perturbations
    Θ = OffsetVector(view(u, 1:(ℓᵧ+1)), 0:ℓᵧ)  # indicies 0 through ℓᵧ
    Θ′ = OffsetVector(view(du, 1:(ℓᵧ+1)), 0:ℓᵧ)
    Θᵖ = OffsetVector(view(u, (ℓᵧ+2):(2ℓᵧ+2)), 0:ℓᵧ)  # indicies 0 through ℓᵧ
    Θᵖ′ = OffsetVector(view(du, (ℓᵧ+2):(2ℓᵧ+2)), 0:ℓᵧ)
    Φ, δ, v, δ_b, v_b = u[(2ℓᵧ+3):(2ℓᵧ+7)]

    # metric perturbations
    Ψ = -Φ - 12H₀² / k^2 / a^2 * (Ω_r * Θ[2])
    Φ′ = Ψ - k^2 / (3ℋₓ^2) * Φ + H₀² / (2ℋₓ^2) * (
        Ω_m * a^(-1) * δ + Ω_b * a^(-1) * δ_b + 4Ω_r * a^(-2) * Θ[0])

    # matter
    δ′ = k / ℋₓ * v - 3Φ′
    v′ = -v - k / ℋₓ * Ψ
    δ_b′ = k / ℋₓ * v_b - 3Φ′

    Θ′[0] = -k / ℋₓ * Θ[1] - Φ′
    if TCA_condition(k, ℋₓ, τₓ′)
        Θ′[2] = 0.0  # could be solved for
        term_3Θ′_vb′ = (
            -((1-R)*τₓ′ + (1+R)*τₓ′′) * (3Θ[1] + v_b) - k * Ψ / ℋₓ
            + (1 - ℋₓ′ / ℋₓ) * (k / ℋₓ) * (-Θ[0] + 2Θ[2]) + k / ℋₓ * (-Θ′[0] + 2Θ′[2])
        ) / ((1+R)*τₓ′ + ℋₓ′ / ℋₓ - 1)
        v_b′ = (-v_b - k * Ψ / ℋₓ + R * (
            term_3Θ′_vb′ + k / ℋₓ * (-Θ[0] + 2Θ[2]) - k / ℋₓ * Ψ
        )) / (1+R)
        Θ′[1] = (term_3Θ′_vb′ - v_b′) / 3
    else
        v_b′ = -v_b - k / ℋₓ * Ψ + τₓ′ * R * (3Θ[1] + v_b)
        Θ′[1] = k / (3ℋₓ) * Θ[0] - 2k / (3ℋₓ) * Θ[2] + k / (3ℋₓ) * Ψ + τₓ′ * (Θ[1] + v_b/3)
    end

    # photons
    Π = Θ[2] + Θᵖ[2] + Θᵖ[0]
    for ℓ in 2:(ℓᵧ-1)
        Θ′[ℓ] = ℓ * k / ((2ℓ+1) * ℋₓ) * Θ[ℓ-1] -
            (ℓ+1) * k / ((2ℓ+1) * ℋₓ) * Θ[ℓ+1] + τₓ′ * (Θ[ℓ] - Π * δ_kron(ℓ, 2) / 10)
    end
    # polarized photons
    Θᵖ′[0] = -k / ℋₓ * Θᵖ[1] + τₓ′ * (Θᵖ[0] - Π / 2)
    for ℓ in 1:(ℓᵧ-1)
        Θᵖ′[ℓ] = ℓ * k / ((2ℓ+1) * ℋₓ) * Θᵖ[ℓ-1] -
            (ℓ+1) * k / ((2ℓ+1) * ℋₓ) * Θᵖ[ℓ+1] + τₓ′ * (Θᵖ[ℓ] - Π * δ_kron(ℓ, 2) / 10)
    end
    # photon hierarchy boundary conditions
    Θ′[ℓᵧ] = k / ℋₓ * Θ[ℓᵧ-1] - (ℓᵧ + 1) / (ℋₓ * bg.η(x)) + τₓ′ * Θ[ℓᵧ]
    Θᵖ′[ℓᵧ] = k / ℋₓ * Θᵖ[ℓᵧ-1] - (ℓᵧ + 1) / (ℋₓ * bg.η(x)) + τₓ′ * Θᵖ[ℓᵧ]

    du[(2ℓᵧ+3):(2ℓᵧ+7)] .= Φ′, δ′, v′, δ_b′, v_b′  # put non-photon perturbations back in
    return nothing
end


text:
using Dorayaki
using Test

@testset "functions" begin
    include("test_arithmetic.jl")
    include("test_tensor.jl")
end

@testset "gradient" begin
    include("test_gradient.jl")
end


text:
# This file is a part of PropDicts.jl, licensed under the MIT License (MIT).

using PropDicts
using Test

@testset "propdict" begin
    da = Dict("foo" => 11, "bar" => Dict("baz" => 42), "44" => "abc")
    @test @inferred(PropDict(da)) isa PropDict

    pa = PropDict(da)
    @test parent(pa) === PropDicts._dict(pa)

    @test convert(Dict, pa) === PropDicts._dict(pa)
    @test convert(Dict{Union{Symbol,Int}}, pa) === PropDicts._dict(pa)
    @test convert(Dict{Union{Symbol,Int},Any}, pa) === PropDicts._dict(pa)

    @test sort(@inferred propertynames(pa)) == [:bar, :foo]
    @test sort(propertynames(pa, true)) == [:_internal_dict, :bar, :foo]

    pb = PropDict(Dict("foo" => 13, :bar => PropDict(:baz => raw"$somevar")))

    @test pa == deepcopy(pa)
    @test pa != pb

    @test pa.foo == pa[:foo] == 11
    @test pa.bar == pa[:bar] == PropDict(:baz => 42)
    @test pa[44] == "abc"

    @test @inferred(merge(pa, pb)) isa PropDict
    pc = merge(pa, pb)
    @test pc.foo == 13
    @test pc.bar == PropDict(:baz => raw"$somevar")
    @test pc[44] == "abc"

    @test PropDicts.contains_vars(raw"fo\\$o") == true
    @test PropDicts.substitute_vars(raw"foo $bar ${baz} y", ignore_missing = true) == raw"foo $bar ${baz} y"
    PropDicts.substitute_vars(raw"foo $bar baz", Dict("bar" => "xyz"))
    s = "xyz"
    PropDicts.substitute_vars(raw"foo$(bar)x$HOME,baz", Dict("bar" => "xyz"), use_env = true)

    PropDicts.substitute_vars!(parent(pc), Dict("somevar" => "xyz"))
end


text:
# Functions that take in preference lists and return modified versions without ties.

"""
    singletiebreaking(arr)

Given schools' ranked preference lists, which may contain ties,
break ties using the single tiebreaking rule by generating
a column of floats, adding this column to each column of `arr`,
and ranking the result columnwise.
"""
function singletiebreaking(arr::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}})
    add = repeat(rand(Float64, size(arr)[1]), 1, size(arr)[2])

    return mapslices(argsort, arr + add, dims=1)
end

"""
    multipletiebreaking(arr)

Given schools' ranked preference lists, which may contain ties,
break ties using the multiple tiebreaking rule by adding to `arr`
a column of random floats having the same shape, then ranking the
result columnwise.
"""
function multipletiebreaking(arr::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}})
    # Given schools' ranked preference lists, which contain ties,
    # breaks ties using the multiple tiebreaking rule by adding a
    # float to each entry and ranking the result columnwise.
    add = rand(Float64, size(arr))

    return mapslices(argsort, arr + add, dims=1)
end

"""
    hybridtiebreaking(arr, blend; return_add)

Given schools' ranked preference lists, which may contain ties,
break ties using a hybrid tiebreaking rule as indicated by
entries of `blend`. `blend` should be a row vector with one entry
on ``[0, 1]`` for each column in `arr`. 0 means that school will use STB,
1 means MTB, and a value in between yields a convex combination
of the rules, which produces interesting results but has not yet
been theoretically analyzed. If `blend` is a scalar, use the same value
at all schools. Undefined behavior for values outside
the ``[0, 1]`` interval.

`return_add` is a `Bool` indicating whether to return the tiebreaking numbers
(lottery numbers) as second entry of output tuple.
"""
function hybridtiebreaking(arr        ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                           blend      ::Union{<:AbstractFloat, AbstractArray{<:AbstractFloat}}=0.;
                           return_add ::Bool=false)
    @assert size(blend) == () || size(blend) == (1, size(arr)[2]) "Dim mismatch between blend and arr"

    add_STB = repeat(rand(Float64, size(arr)[1]), 1, size(arr)[2])
    add_MTB = rand(Float64, size(arr))
    add = (1 .- blend) .* add_STB + blend .* add_MTB

    if return_add
        return mapslices(argsort, arr + add, dims=1), add
    else
        return mapslices(argsort, arr + add, dims=1)
    end
end

"""
    choiceaugmented(arr, targets, blend_target=0, blend_others=0; return_add)

Tiebreaker function for the choice-augmented deferred acceptance
mechanism described by Abdulkadiroğlu et al. (2015). Primary tiebreaking
is accomplished by allowing students to signal a target school; secondary
ties are broken by independent STB lotteries in the target and other
schools (by default). Or, indicate another mechanism by configuring
`blend_target` and `blend_others`, following `?hybridtiebreaking`.

`return_add` is a `Bool` indicating whether to return the tiebreaking numbers
(lottery numbers) as second entry of output tuple.
"""
function choiceaugmented(arr           ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                         targets       ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}},
                         blend_target  ::Union{<:AbstractFloat, AbstractArray{<:AbstractFloat}}=0.,
                         blend_others  ::Union{<:AbstractFloat, AbstractArray{<:AbstractFloat}}=0.;
                         return_add    ::Bool=false)
    @assert (size(arr)[1], ) == size(targets)           "Dim mismatch between arr and targets"
    @assert size(blend_target) == () ||
            size(blend_target) == (1, size(arr)[2])     "Dim mismatch between blend_target and arr"
    @assert size(blend_others) == () ||
            size(blend_others) == (1, size(arr)[2])     "Dim mismatch between blend_others and arr"

    add_STB_target = repeat(rand(Float64, size(arr)[1]), 1, size(arr)[2])
    add_MTB_target = rand(Float64, size(arr))
    add_target = (1 .- blend_target) .* add_STB_target +
                 blend_target .* add_MTB_target
    add_target = mapslices(argsort, add_target, dims=1) / size(arr)[1]

    add_STB_others = repeat(rand(Float64, size(arr)[1]), 1, size(arr)[2])
    add_MTB_others = rand(Float64, size(arr))
    add_others = (1 .- blend_others) .* add_STB_others +
                 blend_others .* add_MTB_others
    add_others /= size(arr)[1]

    add = (1 .+ copy(add_others)) / 2
    for (s, c) in enumerate(targets)
        add[s, c] = add_target[s, c] / 2
    end

    if return_add
        return mapslices(argsort, arr + add, dims=1), add
    else
        return mapslices(argsort, arr + add, dims=1)
    end
end

"""
    welfaretiebreaking(students, schools, blend; equity, return_add)

Given schools' ranked preference lists, which contain ties,
first break ties using student welfare, then break subsequent
ties using a hybrid tiebreaking rule indicated by entries of `blend`.
Or in `equity=true` mode, break ties by minimizing student welfare, which
gives priority in DA to students whose current assignment is poor.
See `?hybridtiebreaking` for an explanation of how to configure `blend`.

`blend=0` is known as the Boston mechanism.

`return_add` is a `Bool` indicating whether to return the tiebreaking numbers
(lottery numbers) as second entry of output tuple.
"""
function welfaretiebreaking(students   ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                            schools    ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                            blend      ::Union{<:AbstractFloat, AbstractArray{<:AbstractFloat}}=0.;
                            equity     ::Bool=false,
                            return_add ::Bool=false)

    @assert size(schools) == size(students')        "Dim mismatch between students and schools"
    @assert size(blend) == () ||
            size(blend) == (1, size(students)[1])   "Dim mismatch between blends and arr"

    add_welfare = equity ? -1 * students' / size(schools)[1] : students' / size(schools)[1]
    add_STB = (1 / size(schools)[1]) *
              repeat(rand(Float64, size(schools)[1]), 1, size(schools)[2])
    add_MTB = (1 / size(schools)[1]) *
              rand(Float64, size(schools))
    add = add_welfare +
          (1 .- blend) .* add_STB + blend .* add_MTB

    if return_add
        return mapslices(argsort, schools + add, dims=1), add
    else
        return mapslices(argsort, schools + add, dims=1)
    end
end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
using BinaryBuilder

name = "LibGit2"
version = v"1.2.3"

# Collection of sources required to build libgit2
sources = [
    GitSource("https://github.com/libgit2/libgit2.git",
              "7f4fa178629d559c037a1f72f79f79af9c1ef8ce"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libgit2*/

atomic_patch -p1 $WORKSPACE/srcdir/patches/libgit2-agent-nonfatal.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/libgit2-hostkey.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/libgit2-continue-zlib.patch

BUILD_FLAGS=(
    -DCMAKE_BUILD_TYPE=Release
    -DTHREADSAFE=ON
    -DUSE_BUNDLED_ZLIB=ON
    "-DCMAKE_INSTALL_PREFIX=${prefix}"
    "-DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}""
)

# Special windows flags
if [[ ${target} == *-mingw* ]]; then
    BUILD_FLAGS+=(-DWIN32=ON -DMINGW=ON -DBUILD_CLAR=OFF)
    if [[ ${target} == i686-* ]]; then
        BUILD_FLAGS+=(-DCMAKE_C_FLAGS="-mincoming-stack-boundary=2")
    fi

    # For some reason, CMake fails to find libssh2 using pkg-config.
    BUILD_FLAGS+=(-Dssh2_RESOLVED=${bindir}/libssh2.dll)
elif [[ ${target} == *linux* ]] || [[ ${target} == *freebsd* ]]; then
    # If we're on Linux or FreeBSD, explicitly ask for mbedTLS instead of OpenSSL
    BUILD_FLAGS+=(-DUSE_HTTPS=mbedTLS -DSHA1_BACKEND=CollisionDetection -DCMAKE_INSTALL_RPATH="\$ORIGIN")
fi
export CFLAGS="-I${prefix}/include"

mkdir build && cd build

cmake .. "${BUILD_FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgit2", :libgit2),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("MbedTLS_jll", v"2.24.0"),
    Dependency("LibSSH2_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function criticalValues(ecdet, typ)
	if ecdet=="none"
		cvals = [6.5, 12.91, 18.9, 24.78, 30.84, 36.25, 42.06, 48.43, 54.01, 59, 65.07, 8.18, 14.9, 21.07, 27.14, 33.32, 39.43, 44.91, 51.07, 57, 62.42, 68.27, 11.65, 19.19, 25.75, 32.14, 38.78, 44.59, 51.3, 57.07, 63.37, 68.61, 74.36, 6.5, 15.66, 28.71, 45.23, 66.49, 85.18, 118.99, 151.38, 186.54, 226.34, 269.53, 8.18, 17.95, 31.52, 48.28, 70.6, 90.39, 124.25, 157.11, 192.84, 232.49, 277.39, 11.65, 23.52, 37.22, 55.43, 78.87, 104.2, 136.06, 168.92, 204.79, 246.27, 292.65]
	elseif ecdet=="const"
		cvals = [7.52, 13.75, 19.77, 25.56, 31.66, 37.45, 43.25, 48.91, 54.35, 60.25, 66.02, 9.24, 15.67, 22, 28.14, 34.4, 40.3, 46.45, 52, 57.42, 63.57, 69.74, 12.97, 20.2, 26.81, 33.24, 39.79, 46.82, 51.91, 57.95, 63.71, 69.94, 76.63, 7.52, 17.85, 32, 49.65, 71.86, 97.18, 126.58, 159.48, 196.37, 236.54, 282.45, 9.24, 19.96, 34.91, 53.12, 76.07, 102.14, 131.7, 165.58, 202.92, 244.15, 291.4, 12.97, 24.6, 41.07, 60.16, 84.45, 111.01, 143.09, 177.2, 215.74, 257.68, 307.64]
	elseif ecdet=="trend"
		cvals = [10.49, 16.85, 23.11, 29.12, 34.75, 40.91, 46.32, 52.16, 57.87, 63.18, 69.26, 12.25, 18.96, 25.54, 31.46, 37.52, 43.97, 49.42, 55.5, 61.29, 66.23, 72.72, 16.26, 23.65, 30.34, 36.65, 42.36, 49.51, 54.71, 62.46, 67.88, 73.73, 79.23, 10.49, 22.76, 39.06, 59.14, 83.2, 110.42, 141.01, 176.67, 215.17, 256.72, 303.13, 12.25, 25.32, 42.44, 62.99, 87.31, 114.9, 146.76, 182.82, 222.21, 263.42, 310.81, 16.26, 30.45, 48.45, 70.05, 96.58, 124.75, 158.49, 196.08, 234.41, 279.07, 327.45]
	end
	cvals = reshape(cvals, 11, 3, 2)
	if typ=="trace"
		return convert(Matrix, cvals[:,:,1])
	elseif typ=="eigen"
		return convert(Matrix, cvals[:,:,2])
	end
end

function selector(i, n)
	v = zeros(n)
	v[i] = 1
	return v
end

function cointegrationJohansen(data, ecdet, K, spec; season = "no", dumvar = "no")
	@assert K >= 2 "K must be at least 2."
	N, P = size(data)
	season == "no" ? s = 0 : s = season - 1
	@assert N * P > P + s * P + K * P^2 + P * (P + 1)/2 "Insufficient degrees of freedom."

	@assert ecdet in ["none"; "const"; "trend"] "Unsupported structure of equation, ecdet must be one of: none, const, trend"
	@assert spec in ["longrun"; "transitory"] "Specification of equation not supported, spec must be one of longrun, transitory"

	nanindexes = !mapslices(x->any(isnan(x)), data, [2])

	if season != "no"
		seasondums = (hcat([selector(mod(x,season)+1,season) for x=1:N]...)')[:,2:season]
		seasondums = seasondums[nanindexes, :]
	end

	if dumvar != "no"
		@assert size(dumvar)[1]==N "Unequal row length between dummy variables and data matrix."
	end

	if ecdet == "trend"
		trend = collect(1:N)
		trend = trend[nanindexes, :]
	end

	data = data[nanindexes, :]
    N, P = size(data)
    Z = lagData(vcat(mapslices(x->diff(x), data, [1]), zeros(P)'), K, N)
    Z0 = Z[:, 1:P]

    spec == "longrun" ? Lnotation = K : Lnotation = 1

    if (ecdet == "none") && (spec == "longrun")
    	ZK = data[1:(N - K), :]
    	Z1 = Z[:, (P+1):(size(Z)[2])]
        Z1 = hcat([1.0 for i in 1:size(Z1)[1]], Z1)
        idx = collect(0:(P - 1))
        model = "with linear trend"
    end
    if (ecdet == "none") && (spec == "transitory")
		ZK = data[1:(N-1), :][K:(N - 1), :]
		Z1 = Z[:, (P+1):(size(Z)[2])]
        Z1 = hcat([1.0 for i in 1:size(Z1)[1]], Z1)
        idx = collect(0:(P - 1))
        model = "with linear trend"
    end
    if (ecdet == "const") && (spec == "longrun")
    	ZK = hcat(data[1:(N - K), :], [1 for i in 1:(N-K)])
    	Z1 = Z[:, (P+1):(size(Z)[2])]
    	P = P + 1
        idx = collect(0:(P - 2))
        model = "without linear trend and constant in cointegration"
    end
    if (ecdet == "const") && (spec == "transitory")
    	ZK = hcat(data[1:(N-1), :], [1 for i in 1:(N-1)])[K:(N - 1), :]
    	Z1 = Z[:, (P+1):(size(Z)[2])]
    	P = P + 1
        idx = collect(0:(P - 2))
        model = "without linear trend and constant in cointegration"
    end
    if (ecdet == "trend") && (spec == "longrun")
		ZK = hcat(data[1:(N - K), :], trend[1:(N - K)])
		Z1 = Z[:, (P+1):(size(Z)[2])]
        Z1 = hcat([1.0 for i in 1:size(Z1)[1]], Z1)
        P = P + 1
        idx = collect(0:(P - 2))
        model = "with linear trend in cointegration"
    end
    if (ecdet == "trend") && (spec == "transitory")
    	ZK = hcat(data[1:(N-1), :], trend[1:(N-1)])[K:(N - 1), :]
    	Z1 = Z[:, (P+1):(size(Z)[2])]
        Z1 = hcat([1.0 for i in 1:size(Z1)[1]], Z1)
        P = P + 1
        idx = collect(0:(P - 2))
        model = "with linear trend in cointegration"
    end
    N = size(Z0)[1]

    if season != "no"
    	if ecdet == "const"
    		Z1 = hcat(seasondums[(K+1):size(seasondums)[1],:], Z1)
    	else
    		Z1 = hcat(Z1[:, 1], seasondums[(K+1):size(seasondums)[1],:], Z1[:, 2:size(Z1)[2]])
    	end
    end

    if dumvar != "no"
    	if ecdet == "const"
    		Z1 = hcat(dumvar[(K+1):size(dumvar)[1],:], Z1)
    	else
    		Z1 = hcat(Z1[:, 1], dumvar[(K+1):size(dumvar)[1],:], Z1[:, 2:size(Z1)[2]])
    	end
    end

    M00 = (Z0'*Z0)/N
    M11 = (Z1'*Z1)/N
    MKK = (ZK'*ZK)/N
    M01 = (Z0'*Z1)/N
    M0K = (Z0'*ZK)/N
    MK0 = (ZK'*Z0)/N
    M10 = (Z1'*Z0)/N
    M1K = (Z1'*ZK)/N
    MK1 = (ZK'*Z1)/N
    M11inv = inv(M11)
    R0 = Z0 - (M01 * M11inv * Z1')'
    RK = ZK - (MK1 * M11inv * Z1')'
    S00 = M00 - M01 * M11inv * M10
    S0K = M0K - M01 * M11inv * M1K
    SK0 = MK0 - MK1 * M11inv * M10
    SKK = MKK - MK1 * M11inv * M1K
    Ctemp = cholfact(SKK, :U, Val{true})
    pivot = sortperm(Ctemp.piv)
    C = (convert(Matrix,Ctemp[:U])[:, pivot])'
    Cinv = inv(C)
    S00inv = inv(S00)
    λ, e = eig(Cinv * SK0 * S00inv * S0K * Cinv')
    V = Cinv' * e
    Vorg = V
    V = hcat([V[:, x]./V[1, x] for x in 1:P]...)
    W = S0K * V * inv(V' * SKK * V)
    Π = S0K * inv(SKK)
    Δ = S00 - S0K * V * inv(V' * SKK * V) * V' * SK0
    Γ = M01 * M11inv - Π * MK1 * M11inv

    return (data, Z0, Z1, ZK, dumvar, λ, Vorg, V, W, Π, Δ, Γ, R0, RK, spec, season)
end

function johansenTeststat(cointEst, typ)
	@assert typ in ["eigen"; "trace"] "Unsupported type, typ must be one of: eigen, trace"
	P = size(cointEst.ZK)[2]
	N = size(cointEst.data)[1]

	if (ecdet == "none") && (spec == "longrun")
        idx = collect(0:(P - 1))
    end
    if (ecdet == "none") && (spec == "transitory")
        idx = collect(0:(P - 1))
    end
    if (ecdet == "const") && (spec == "longrun")
        idx = collect(0:(P - 2))
    end
    if (ecdet == "const") && (spec == "transitory")
        idx = collect(0:(P - 2))
    end
    if (ecdet == "trend") && (spec == "longrun")
        idx = collect(0:(P - 2))
    end
    if (ecdet == "trend") && (spec == "transitory")
        idx = collect(0:(P - 2))
    end

	if P > 11
		warning("Too many variables, critical values cannot be computed.")
		cvals = []
	else
		cvals = criticalValues(cointEst.ecdet, typ)
	end

	if typ == "trace"
    	teststat = reverse(map(x -> -N * sum(log(1-sort(cointEst.λ, rev = true)[x:P])), idx + 1))
    elseif typ == "eigen"
    	teststat = reverse( (-N*log(1-sort(cointEst.λ, rev = true)))[idx + 1])
    end
    return (teststat, cvals)
end

type CointegrationEstimate
	data::Matrix{Float64}
	Z0::Matrix{Float64}
	Z1::Matrix{Float64}
	ZK::Matrix{Float64}
    dumvar
    λ::Vector{Float64}
    Vorg::Matrix{Float64}
    V::Matrix{Float64}
    W::Matrix{Float64}
    PI::Matrix{Float64}
    DELTA::Matrix{Float64}
    GAMMA::Matrix{Float64}
    R0::Matrix{Float64}
    RK::Matrix{Float64}
    spec::AbstractString
    season
    ecdet::AbstractString

	function CointegrationEstimate(data::Matrix{Float64}, ecdet::AbstractString, K::Int, spec::AbstractString; season = "no", dumvar = "no")
		data, Z0, Z1, ZK, dumvar, λ, Vorg, V, W, PI, DELTA, GAMMA, R0, RK, spec, season = cointegrationJohansen(data, ecdet, K, spec; season = season, dumvar = dumvar)
		new(data, Z0, Z1, ZK, dumvar, λ, Vorg, V, W, PI, DELTA, GAMMA, R0, RK, spec, season, ecdet)
	end
end

function getVARrepresentation(cointObject, r::Int)
	@assert typeof(cointObject)==CointegrationEstimate "The first argument should be cointegration estimate."
	@assert r < size(cointObject.data)[2] "The cointegration rank---the second argument---can be no bigger than number of variables in the system."

	N, vars = size(cointObject.data)

	etc = cointObject.ZK * cointObject.V[:, 1:r]
	coeffs = hcat(etc, cointObject.Z1) \ cointObject.Z0

	Π = cointObject.W[:, 1:r] * cointObject.V[:,1:r]'

	if cointObject.ecdet == "const"
		detcoeffs = Π[:, vars + 1]
		Π = Π[:, 1:vars]
		rhs = hcat([1 for i in 1:size(cointObject.ZK)[1]], cointObject.Z1)
	elseif cointObject.ecdet == "none"
		detcoeffs = coeffs[r+1,:]'
		rhs = cointObject.Z1
	elseif cointObject.ecdet == "trend"
		detcoeffs = hcat(coeffs[r+1,:]',Π[:,vars + 1])
		Π = Π[:, 1:vars]
		rhs = hcat([1 for i in 1:size(cointObject.ZK)[1]], cointObject.ZK[:, vars + 1], cointObject.Z1[:, 2:size(cointObject.Z1)[2]])
	end

	cointObject.ecdet == "const" ? a = 0 : a = 1
	if cointObject.season != "no"
		s = cointObject.season - 1
		coefs_from_coeffs = coeffs[(r+a) + 1:s, :]
		detcoeffs = vcat(detcoeffs, coefs_from_coeffs)
	else
		s = 0
	end
	if cointObject.dumvar != "no"
		d = size(cointObject.dumvar)[2]
		coefs_from_coeffs = coeffs[(r+a+s) + 1:d, :]
		detcoeffs = vcat(detcoeffs, coefs_from_coeffs)
	else
		d = 0
	end

	detcoeffs = detcoeffs'

	lags = N-size(cointObject.ZK)[1]
	Γ = coeffs[((-(lags-1)*vars+1):0)+size(coeffs)[1],:]'
	A = zeros(vars, vars, lags)

	i = hcat([vars*i + 1 for i=0:(lags-2)], [vars + vars*i for i=0:(lags-2)])
	if cointObject.spec == "transitory"
		A[:,:,1] = Γ[:, i[1,1]:i[1,2]] + Π + eye(vars)
		for j=2:(lags-1)
			A[:,:,j] = Γ[:, i[j,1]:i[j,2]] - Γ[:,i[j - 1,1]:i[j - 1,2]]
		end
		A[:,:,lags] = -Γ[:, i[lags-1,1]:i[lags-1, 2]]
	elseif cointObject.spec == "longrun"
		A[:,:,1] = Γ[:, i[1,1]:i[1,2]] + eye(vars)
		for j=2:(lags-1)
			A[:,:,j] = Γ[:, i[j,1]:i[j,2]] - Γ[:,i[j - 1,1]:i[j - 1,2]]
		end
		A[:,:,lags] = Π-Γ[:, i[lags-1,1]:i[lags-1, 2]]
	end

	determin_data = rhs[:,1:size(detcoeffs)[1]]
	data = lagData(vcat(cointObject.data, zeros(vars)'), lags, N)
	Y = (cointObject.data)[(lags+1):size(cointObject.data)[1],:]
	C = vcat(detcoeffs, vcat([A[:,:,i]' for i=1:size(A)[3]]...))
	X = hcat(determin_data, data)

	resids = Y - X * C

	Σ = cov(resids)*(N-1)/(N-lags*vars-size(detcoeffs)[1])
	return (lags, size(detcoeffs)[1], vars, N, X, Y, C, Σ)
end

type varRepresentationVECM <: VARRepresentation
	lags::Int
	typ::AbstractString
	ntyp::Int
	vars::Int
	obs::Int
	X::Matrix{Float64}
	Y::Matrix{Float64}
	C::Matrix{Float64}
	Σ::Matrix{Float64}
	r::Int
	HPsi::Int
	HPhi::Int
	Phi::Array{Float64, 3}
	Psi::Array{Float64, 3}
	seriesNames::Vector{AbstractString}

	function varRepresentationVECM(cointegrationEstimate, r::Int)
		(lags, ntyp, vars, obs, X, Y, C, Σ) = getVARrepresentation(cointegrationEstimate, r::Int)
		new(lags, "varRepresentation", ntyp, vars, obs, X, Y, C, Σ, r, 0, 0, zeros(vars,vars,1), zeros(vars,vars,1), ["" for i=1:vars])
	end
end


text:
using Iterators
import Iterators: PeekIter
import Base: start, next, done

immutable UniqueIter
    inner::PeekIter
    UniqueIter(inner::PeekIter) = new(inner)
    UniqueIter(inner) = new(PeekIter(inner))
end

# TODO: eltype, iteratorsize, length, size

function start(it::UniqueIter)
    return start(it.inner)
end

function next(it::UniqueIter, s)
    hd, s = next(it.inner, s)
    while !done(it.inner, s) && hd == get(peek(it.inner, s))
        hd, s = next(it.inner, s)
    end
    return hd, s
end

function done(it::UniqueIter, s)
    return done(it.inner, s)
end

function unique(it)
    UniqueIter(it)
end


text:
function get_timeseries(sol)
    sol_copy = deepcopy(sol)

    for i = eachindex(sol_copy.u)
        while length(sol_copy.u[i]) < length(sol_copy.u[end])
            ind = length(sol_copy.u[i]) - sol_copy.prob.p.M + 1
            insert!(sol_copy.u[i], ind, 0.0)
        end
    end

    return(hcat(sol_copy.u...)')
end

function get_Nsp(sol, abstol::Float64 = 1e-10)
    series = get_timeseries(sol)
    return(mapslices(x -> sum(x .> abstol), series[: , 1:sol.prob.p.N], dims = 2))
end

#respiration callback
function save_func(u, t, integrator)
    integrator.p.Rm .* integrator.u[(integrator.p.M+1 : end)]
end

###
#Interactions
###

#get resource explicit competition
function get_a_ij_comp(y,u,l_sum,i,j,k)
    y[k] * u[i,k] * (u[j,k] * (1 - l_sum[k]))
end

function get_a_ij_facil(y,u,l,l_sum,i,j,k,M)
    sum_var = 0.0
    for l_ = 1:M
        sum_var += y[l_] * u[i,l_] * l[l_,k]
    end
    return(sum_var * (u[j,k] * (1 - l_sum[k])))
end

function get_resource_a(sol, t, type::Symbol)
    #assert the interaction type is valid
    @assert type ∈ [:competition, :facilitation, :both] "type of interaction must be competition, facilitation or both"
    #get resource enviroment and parameter object
    p = sol.prob.p
    y = sol(t)[1:p.M]

    #loop over all species pairs and resources
    a = zeros(p.N,p.N)
    for i = 1:p.N, j = 1:p.N
        for k = 1:p.M
            if type == :competition
                a[i,j] -= get_a_ij_comp(y,p.u,p.l_sum,i,j,k)
            elseif type == :facilitation
                a[i,j] += get_a_ij_facil(y,p.u,p.l,p.l_sum,i,j,k,p.M)
            elseif type == :both
                a[i,j] += get_a_ij_facil(y,p.u,p.l,p.l_sum,i,j,k,p.M) - get_a_ij_comp(y,p.u,p.l_sum,i,j,k)
            end

        end
    end
    return(a)
end

#naive interaction measures
#competition - degree to which the uptake of i overlaps with the uptake of j
function get_overlap_competition(sol)
    p = sol.prob.p
    a = zeros(p.N,p.N)
    for i = 1:p.N, j = 1:p.N
        t = 0.0
        b = 0.0
        for k = 1:p.M
            t += min(p.u[i,k],p.u[j,k])
            b += max(p.u[i,k],p.u[j,k])
        end
        a[i,j] = t / b
    end
    return(a)
end

#facilitation - degree to which the leaked uptake of i overlaps with the uptake of j
function get_overlap_facilitation(sol)
    p = sol.prob.p
    a = zeros(p.N,p.N)
    for i = 1:p.N, j = 1:p.N
        t = 0.0
        b = 0.0
        for k = 1:p.M
            t += min(sum(p.u[i,:] .* p.l[:,k]), p.u[j,k])
            b += max(sum(p.u[i,:] .* p.l[:,k]), p.u[j,k])
        end
        a[i,j] = t / b
    end
    return(a)
end

#get respiration
function get_R(sol,t)
    sum(sol(t)[(sol.prob.p.M+1):end] .* sol.prob.p.Rm)
end

#get outflow
function get_ω(sol,t)
    sum(sol(t)[1:sol.prob.p.M] .* sol.prob.p.ω)
end


text:
# Point in 2D space

export Point2D

mutable struct Point2D <: Point
    x :: Float64
    y :: Float64
end

point(x,y) = Point2D(x,y)
zero(::Type{Point2D}) = Point2D(0.0, 0.0)
function rand(r::AbstractRNG, ::Random.SamplerType{Point2D})
    t = 2pi*rand(r)
    Point2D(cos(t),sin(t))
end

+(lhs :: Point2D, rhs :: Point2D) = Point2D(lhs.x + rhs.x, lhs.y + rhs.y)
-(lhs :: Point2D, rhs :: Point2D) = Point2D(lhs.x - rhs.x, lhs.y - rhs.y)
-(p :: Point2D) = Point2D(-p.x, -p.y)
*(a :: Real, p :: Point2D) = Point2D(a*p.x, a*p.y)
/(p :: Point2D, a :: Real) = Point2D(p.x/a, p.y/a)

abs(p :: Point2D) = hypot(p.x, p.y)
abs2(p :: Point2D) = p.x*p.x + p.y*p.y

plot_str(p :: Point2D) = string(p.x, " ", p.y)

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
constants=[:me≖9.1094e-31,:qe≖1.6022e-19,:kB≖1.3807e-23,:ε0≖8.8542e-12]

# Solar wind
ne=:ne≖10.0^7
Te=:Te≖10.0*:qe
variables=[ne,Te]
# Debye length lD
lD=:λD≖sqrt(:ε0*:Te/(:ne*:qe^2)) #valid when Te>>Ti
# plasma frequency wp
wp=:ωp≖sqrt(:ne*:qe^2/(:ε0*:me))
# particles in a debye cube ND
ND=:ND≖:ne*:λD^3

@test (lD&variables&constants).rhs ≈ 7.433892903446526
@test (wp&lD&variables&constants).rhs ≈ 178400.95597166813
@test (ND&lD&variables&constants).rhs ≈ 4.108174668936227e9

# Drifts
ve=:ve≖Cross(:E⊥,:B)/norm(:B)^2
ex=[1,0,0];ey=[0,1,0];ez=[0,0,1]
B=:B≖:Bt*Vec(ez)
Bt=:Bt≖1.5
r=:r≖1
dBt=:dBt≖0.3
Eperp=:E⊥≖-:r/2*:dBt*Vec(ey)
vars=[B,Bt,Eperp,r,dBt]

@test (ve&vars).rhs[1][1] ≈ -0.1
@test (ve&vars).rhs[1][2]==Vec([1,0,0])

# Temperatures
equ=Equation
Ekin=equ(:Ekin,0.5*:Tpa+:Tpe)
parfac=1
Tpa1=equ(:Tpa1,parfac*:T0) #is parallel temperature affected?
Tpe1=equ(:Tpe1,2*:T0)
Ek1=equ(:Ek1,0.5*:Tpa1+:Tpe1)
Ek2=equ(:Ek1,0.5*:T2+:T2)
Ek1s=Ek1&Tpa1&Tpe1 #s=solved
T2=(Ek2&Ek1s)/1.5
T2=equ(T2.rhs,T2.lhs)
Tpa3=equ(:Tpa3,1/parfac*:T2) #is it?
Tpe3=equ(:Tpe3,:T2/2)
Ek3=equ(:Ek3,0.5*:Tpa3+:Tpe3)
Ekf=[equ(:Ekf,0.5*:Tf+:Tf),equ(:Ekf,:Ek3)]
Tf=(Ekf[1]&Ekf[2]&Ek3&Tpa3&Tpe3&T2)/1.5

@test T2.rhs[1][1] ≈ 1.6666666666666665
@test Tf.lhs[1][1] ≈ 1.111111111111111


text:
using Dates
using ManpowerPlanning
using ManpowerPlanningIO
using Plots
using Test
using XLSX

versionMPIO()
println()

MPsim = ManpowerSimulation
MPIO = ManpowerPlanningIO

tStart = now()

@testset "Manpower Simulation I/O tests" begin

@testset "Excel configuration tests" begin

include( "base.jl" )
include( "attrition.jl" )
include( "attribute.jl" )
include( "basenode.jl" )
include( "compoundnode.jl" )
include( "recruitment.jl" )
include( "transition.jl" )
include( "retirement.jl" )
include( "initpop.jl" )

end  # @testset "Excel configuration tests"

@testset "Excel run tests" begin

include( "randomseed.jl" )

end  # @testset "Excel run tests"

@testset "Database configuration tests" begin

include( "db/complete.jl" )

end  # @testset "Database configuration tests" begin

include( "report/simprep.jl" )

@testset "Simulation report tests" begin

rm.( string.( "report/", filter( fname -> endswith( fname, ".xlsx" ),
    readdir( "report" ) ) ) )
include( "report/transitionflux.jl" )
include( "report/nodeflux.jl" )
include( "report/nodepop.jl" )
include( "report/nodecomposition.jl" )
include( "report/subpop.jl" )
include( "report/initpop.jl" )

end  # @testset "Simulation report tests"

@testset "Simulation plot tests" begin

gui( plot() )

rm.( string.( "plot/", filter( fname -> endswith( fname, ".svg" ),
    readdir( "plot" ) ) ) )
include( "plot/transitionflux.jl" )
include( "plot/nodeflux.jl" )
include( "plot/nodepop.jl" )
include( "plot/nodecomposition.jl" )
include( "plot/subpop.jl" )
include( "plot/initpop.jl" )

end  # @testset "Simulation plot tests"

@testset "Excel plot request tests" begin

rm( "excelrequest/plotRequest", recursive=true, force=true )
include( "excelrequest/transitionflux.jl" )
include( "excelrequest/nodepop.jl" )
include( "excelrequest/initpop.jl" )

end  # @testset "Excel plot request tests"

end  # @testset "Manpower Simulation I/O tests"

tElapsed = ( now() - tStart ).value / 1000
@info string( "Unit tests completed in ", tElapsed, " seconds." )


text:
module SIMD

using Base: @propagate_inbounds

export Vec, vload, vloada, vloadnt, vloadx, vstore, vstorea, vstorent, vstorec,
       vgather, vgathera, vscatter, vscattera, shufflevector, vifelse, valloc,
       VecRange

const VE         = Base.VecElement
const LVec{N, T} = NTuple{N, VE{T}}

const IntTypes      = Union{Int8, Int16, Int32, Int64} # Int128 and UInt128 does not get passed as LLVM vectors
const BIntTypes     = Union{IntTypes, Bool}
const UIntTypes     = Union{UInt8, UInt16, UInt32, UInt64}
const IntegerTypes  = Union{IntTypes, UIntTypes}
const BIntegerTypes = Union{IntegerTypes, Bool}
const FloatingTypes = Union{Float32, Float64} # Float16 support is non-native in Julia and gets passed as an i16
const ScalarTypes   = Union{IntegerTypes, FloatingTypes}
const VecTypes      = Union{ScalarTypes, Ptr, Bool}

if VERSION < v"1.6.0-DEV.674"
    include("LLVM_intrinsics_old.jl")
else
    include("LLVM_intrinsics.jl")
end

include("simdvec.jl")
include("arrayops.jl")

end


text:
class Point2 {
  int x;
  int y;

  void move (int dx, int dy) {
     x = x + dx;
     y = y + dy;
  }

  int getX () { return x; }

  int getY () { return y; }
}

class Point3 extends Point2 {
  int z;

  void moveZ (int dz) {
    z = z + dz;
  }

  int getZ () { return z; }

}

class Point4 extends Point3 {
  int w;

  void moveW (int dw) {
    w = w + dw;
  }

  int getW () { return w; }

}

int main () {
  Point2 p = new Point3;

  Point3 q = new Point3;

  Point4 r = new Point4;

  q.move(2,4);
  q.moveZ(7);
  p = q;

  p.move(3,5);

  r.move(1,3);
  r.moveZ(6);
  r.moveW(2);

  printInt(p.getX());
  printInt(p.getY());
  printInt(q.getZ());
  printInt(r.getW());
  return 0;

}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
abstract type MSetSolver end
function build!(solver::MSetSolver, reward::Matrix{Distribution})
  nothing
end

struct MSet <: CombinatorialInstance{Int}
  # Mandatory properties.
  n_arms::Int
  optimal_average_reward::Float64
  m::Int

  # Probability distributions for the arm rewards.
  reward::Vector{Distribution}

  # Internal solver.
  solver::MSetSolver

  function MSet(reward::Vector{Distribution}, m::Int, solver::MSetSolver)
    # Basically, any combination of parameters is allowable, even if m is
    # larger than the number of arms.

    # Prepare the solver to be used (if required).
    build!(solver, m, length(reward))

    avg_reward = Dict{Int, Float64}(i => mean(reward[i...]) for i in all_arm_indices(reward))
    opt_sol = solve_linear(solver, avg_reward)
    opt = sum(avg_reward[i] for i in opt_sol)

    # Done!
    return new(length(reward), opt, m, reward, solver)
  end
end

function is_feasible(instance::MSet, arms::Vector{Int})
  if length(arms) > instance.m
    return false
  end

  if length(arms) != length(unique(arms))
    return false
  end

  return true
end

function is_partially_acceptable(instance::MSet, arms::Vector{Int})
  return is_feasible(instance, arms)
end


text:
#= -------------------------------------------------------------------------
# @file practica1_problema3.jl
#
# @date 09/22/17 13:41:52
# @author Martin Noblia
# @email martin.noblia@openmailbox.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using DifferentialEquations
using Plots; pyplot()

#---------------------------------------------------------------------------
# constants
#---------------------------------------------------------------------------
const ϵ = 0.8
const Rₑ = 0.4
const D = 1e-3
const T_alr = 300
const T_inf = 300
const h = 100
const σ = 5.67e-8
const ρ = 8933
const c = 385
I₁(t) = (t ≥ 0) ? 5.2 : 0.0 # step function
I₂(t) = t
I₃(t) = sin(t)
#differential equation
f(t, T) = (I₁(t)^2 * Rₑ - π * D * h * (T - T_inf) - π * D * ϵ * σ *(T^4 - T_alr^4)) / (ρ * c *(π * (D^2)/4))
T₀ = 25.0 + 273.0
tsimu = (0.0, 240.0)
prob  = ODEProblem(f, T₀, tsimu)
solution = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

plot(solution.t, solution - 273)


text:
using Test
using OMEinsum
using OMEinsum: get_size_dict

@testset "unspecified index sizes" begin
    v = randn(5)
    a = randn(6,6)
    xs = (v, a)
    ixs = (('j',), ('k', 'k'))
    @test einsum(ein"j,kk->j", xs, get_size_dict(ixs, xs)) ≈ ein"j,kk->j"(xs...)
    res = ein"j,kk->iij"(xs...; size_info=Dict('i'=>9))
    @test size(res) == (9,9,5)
    @test ein"ijk,
    ijk->
    ikl" == ein"ijk,ijk->ikl"
end


text:
# Test uniform contexts
d, n = 5, 10
csampler = UniformContexts(repeat([10.0 11.0], d, 1))
X = csampler(n)
@test size(X) == (d, n)
@test all(10 .<= X .<= 11)


text:
#!/usr/bin/env julia

using HTTP, JSON3, Scratch, TimeZones, Dates, Tar, Gzip_jll

# Get list of Pkg servers
pkg_servers = JSON3.read(String(HTTP.get("https://pkg.julialang.org/meta/siblings").body))
#storage_servers = JSON3.read(String(HTTP.get("https://pkg.julialang.org/meta/parents").body))
storage_servers = ["https://us-east.storage.juliahub.com", "https://kr.storage.juliahub.com"]

function update_git_repo(repo_url::AbstractString, repo_path::AbstractString; branch::AbstractString = "master")
    mkpath(repo_path)
    st = stat(joinpath(repo_path,"HEAD"))
    if !isfile(st)
        run(`git clone --mirror $(repo_url) $(repo_path)`)
    end

    # Return the latest fetched treehash
    run(`git -C $(repo_path) fetch --all`)
    return readchomp(`git -C $(repo_path) rev-parse "$(branch)^{tree}"`)
end

# Clone/update General registry
gen_uuid = Base.UUID("23338594-aafe-5451-b93e-139f81909106")
reg_dir = get_scratch!("Registry-General", Base.UUID("eac38ba3-4627-46d4-b1a1-bcb86ba22f8b"))
latest_treehash = update_git_repo("https://github.com/JuliaRegistries/General", reg_dir)

function commit_for_tree(repo_path, tree)
    p = open(`git -C $(repo_path) log --pretty='%H %T'`; read=true)
    for line in eachline(p.out)
        commit_hash, commit_tree = split(line)[1:2]
        if commit_tree == tree
            return commit_hash
        end
    end
    return nothing
end
function commit_time(repo_path, commit)
    time_str = readchomp(`git -C $(repo_path) show --quiet --pretty='%cI' $(commit)`)
    @show time_str
    ZonedDateTime(time_str, "yyyy-mm-ddTHH:MM:SSzzzzz")
end

hash_map = Dict()
function get_general_time(treehash)
    if haskey(hash_map, treehash)
        return hash_map[treehash]
    end

    commit = commit_for_tree(reg_dir, treehash)
    if commit === nothing
        # If we couldn't find this commit, it's probably because the Registry contains an empty directory or something.
        # This causes us to generate a different hash than git would.  We get the git-compatible hash by downloading the
        # registry then using `Tar` to generate a yskip hash:
        local skip_treehash
        mktempdir() do dir
            @info("Commit for treehash $(treehash) could not be found, determining skip treehash...")
            tarball = joinpath(dir, treehash)
            tarball = HTTP.download("https://pkg.julialang.org/registry/$(gen_uuid)/$(treehash)", tarball)
            gzip() do gzip
                p = open(pipeline(`cat $(tarball)`, `gzip -d`), read=true)
                skip_treehash = Tar.tree_hash(p.out; skip_empty=true)
                close(p.in)
                @info("skip_treehash: $(skip_treehash)")
            end
            prepare_for_deletion(dir)
        end

        commit = commit_for_tree(reg_dir, skip_treehash)
        if commit === nothing
            error("Unable to find commit for General hash $(treehash)/$(skip_treehash)")
        end
    end

    hash_map[treehash] = commit_time(reg_dir, commit)
    return hash_map[treehash]
end

latest_time = get_general_time(latest_treehash)
function get_server_registry_lag(server)
    local registries
    try
        registries = split(String(HTTP.get("$(server)/registries").body), "\n")
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        @warn("Unable to fetch `/registries` from $(server)")
    end
    treehash = basename(first(filter(r -> startswith(r, "/registry/$(gen_uuid)/"), registries)))
    time = get_general_time(treehash)
    return latest_time - time
end

# Get the current registry hash for each server:
@info("Storage Servers:")
for server in storage_servers
    @info(server, lag=get_server_registry_lag(server))
end

@info("Pkg Servers:")
for server in pkg_servers
    @info(server, lag=get_server_registry_lag(server))
end