text:
module mfiber_reinf_cant_yn_strong_Abaqus
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import LinearAlgebra: Symmetric, cholesky
import Statistics: mean
function test()

# 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 = 2
na = 4*n # number of elements lengthwise
nb = 2*n # number of elements through the depth
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("Mesh generation")
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

femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), 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.femmbase, 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("Stiffness")
K = stiffness(femm, geom, u)
fi = ForceIntensity(FFlt, 3, getshr!);
# println("Traction loads")
F =  distribloads(Tracfemm, geom, u, fi, 2);

# println("Factorization")
K = (K + K')/2;
K = cholesky(Symmetric(K))
# println("U = K\\F")
U = K\F
# # println("U = cg(K, F; tol=1e-3, maxiter=2000)")
# 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)")
@test abs(utip - -0.00653888266072445)/abs(utip) < 1.0e-6
# println("Solution: $(  time()-t0 )")

AE = AbaqusExporter("fiber_reinf_cant_yn_strong_Abaqus");
HEADING(AE, "fiber_reinf_cant_yn_strong_Abaqus.");
PART(AE, "part1");
END_PART(AE);
ASSEMBLY(AE, "ASSEM1");
INSTANCE(AE, "INSTNC1", "PART1");
NODE(AE, fens.xyz);
ELEMENT(AE, "c3d20r", "AllElements", connasarray(femm.integdomain.fes))
ELEMENT(AE, "SFM3D8", "TractionElements", connasarray(Tracfemm.integdomain.fes))
NSET_NSET(AE, "L1", lx0)
NSET_NSET(AE, "L2", lx0)
NSET_NSET(AE, "L3", lx0)
NSET_NSET(AE, "tip", Tipl)
ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
ORIENTATION(AE, "MaterialOrientation", vec(csmat[:,1]), vec(csmat[:,2]));
SOLID_SECTION(AE, "elasticity", "MaterialOrientation", "AllElements", "Hourglassctl");
SURFACE_SECTION(AE, "TractionElements")
END_INSTANCE(AE);
END_ASSEMBLY(AE);
MATERIAL(AE, "elasticity")
ELASTIC(AE, E1s, E2s, E3s, nu12s, nu13s, nu23s, G12s, G13s, G23s)
SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
STEP_PERTURBATION_STATIC(AE)
BOUNDARY(AE, "ASSEM1.INSTNC1.L1", 1)
BOUNDARY(AE, "ASSEM1.INSTNC1.L2", 2)
BOUNDARY(AE, "ASSEM1.INSTNC1.L3", 3)
DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec([0.0, 0.0, q0]))
NODE_PRINT(AE, "ASSEM1.INSTNC1.tip")
ENERGY_PRINT(AE)
END_STEP(AE)
close(AE)
try rm(AE.filename); catch end

# println("Done: $(  time()-t0 )")
true
end
end
using .mfiber_reinf_cant_yn_strong_Abaqus
mfiber_reinf_cant_yn_strong_Abaqus.test()

module mmorthoballoonpenaltymm
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()

    # Orthotropic balloon inflation, axially symmetric model

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    tolerance=rin/1000.

    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,pi/2,5,20);
    bdryfes = meshboundary(fes);

    icl = selectelem(fens, bdryfes, box=[0.,0.,0.,pi/2],inflate=tolerance);
    for i=1:count(fens)
        r=rin+fens.xyz[i,1]; a=fens.xyz[i,2];
        fens.xyz[i,:]=[r*cos(a) r*sin(a)];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # the symmetry plane
    ly =selectelem(fens, bdryfes; box=[0 rex 0 0], inflate = tolerance)
    # the axis of symmetry
    lx =selectelem(fens, bdryfes; box=[0 0 0 rex], inflate = tolerance)

    # No EBC
    applyebc!(u)
    numberdofs!(u)
    # println("Number of degrees of freedom = $(u.nfreedofs)")

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(1, 3), true))
    function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
      copyto!(forceout, XYZ/norm(XYZ)*p)
      return forceout
    end
    fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
    F2= distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    ##
    # The restraints of the nodes on the bounding cross-sections in the direction
    # of the normal to the plane of the cross-section  in the
    # circumferential direction are introduced using a penalty formulation.
    # For that purpose we introduce  a finite element model machine for the
    # surface  finite elements on the cross-sections.
    springcoefficient =1.0e9/ (abs(p)/E1)
    xsfemm = FEMMDeforWinkler(IntegDomain(subset(bdryfes, lx), GaussRule(1, 3), true))
    ysfemm = FEMMDeforWinkler(IntegDomain(subset(bdryfes, ly), GaussRule(1, 3), true))
    H = surfacenormalspringstiffness(xsfemm,  geom, u, springcoefficient, SurfaceNormal(3)) +
        surfacenormalspringstiffness(ysfemm,  geom, u, springcoefficient, SurfaceNormal(3))
    K =stiffness(femm, geom, u)
    U=  (K + H)\(F2)
    scattersysvec!(u,U[:])

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 3)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
@test abs(minimum(fld.values) - (-0.04635309320688638)) < 1.0e-5
@test abs(maximum(fld.values) - (0.5708149883384825)) < 1.0e-5
    # File =  "orthoballoon_penalty_sigmaz.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    #               vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mmorthoballoonpenaltymm
mmorthoballoonpenaltymm.test()

module mbar1
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.FENodeSetModule
using FinEtools.MeshExportModule
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    Area = 2.0*phun("in^2")
    E = 30e6*phun("psi") # Young's modulus
    nu = 0.0
    alpha = 5.5e-6*phun("in")/phun("in")/phun("F")
    fens = FENodeSetModule.FENodeSet([10. -5 20;
    30 25 -15]*phun("in") )
    fes = FESetL2(reshape([1,2],1,2))

    # function otherdimensionfu(loc::FFltMat,
    #   conn::CC, N::FFltMat)::FFlt where {CC<:AbstractArray{FInt}}
    #   return otherdimension::FFlt
    # end
    integdomain = IntegDomain(fes, GaussRule(1, 2), (loc, conn, N) -> Area, false)
    # display(IntegDomain)

    MR = DeforModelRed1D
    material = MatDeforElastIso(MR,  0.0, E, nu, alpha)
    # display(material )

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

    femm = FEMMDeforLinear(MR, integdomain, CSys(3, 1), material)
    K = stiffness(femm,  geom,  u)
    # println("K = $(K/(phun("lbf")/phun("in")))")
    ref_K=   1.0e+05 *[
    1.8916    2.8373   -3.3102   -1.8916   -2.8373    3.3102
    2.8373    4.2560   -4.9653   -2.8373   -4.2560    4.9653
   -3.3102   -4.9653    5.7929    3.3102    4.9653   -5.7929
   -1.8916   -2.8373    3.3102    1.8916    2.8373   -3.3102
   -2.8373   -4.2560    4.9653    2.8373    4.2560   -4.9653
    3.3102    4.9653   -5.7929   -3.3102   -4.9653    5.7929];
    @test norm(K/(phun("lbf")/phun("in")) - ref_K)/1.0e5 < 1.0e-3

    dT = NodalField(broadcast(+, zeros(size(fens.xyz, 1), 1),
        100*phun("F"))) # temperature field
    # display(dT)

    F2 = thermalstrainloads(femm, geom, u, dT)
    # println("F2 = $(F2/(phun("lbf")))")
    ref_F=  1.0e+04 *[
    -1.313449091077187
    -1.970173636615779
    2.298535909385076
    1.313449091077187
    1.970173636615779
    -2.298535909385076]
    @test norm(F2/(phun("lbf")) - ref_F) < 1.0e-2

   # K = cholesky(K)
   # U=  K\(F2)
   # scattersysvec!(u, U[:])

    # File = "playground.vtk"
    # MeshExportModule.VTK.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

end
end
using .mbar1
mbar1.test()

module mbar2
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.FENodeSetModule
using FinEtools.MeshExportModule
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    Area = 1.5
    E = 1.0e7 # Young's modulus
    nu = 0.0
    alpha = 0.0
    fens = FENodeSetModule.FENodeSet([0.0 0;
    0 40;
    40 0;
    40 40;
    80 0;
    80 40] )
    fes = FESetL2([1     3
     1     4
     2     4
     3     4
     3     5
     5     4
     6     4
     5     6])

    MR = DeforModelRed1D
    material = MatDeforElastIso(MR,  0.0, E, nu, alpha)
    # display(material )

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field
    setebc!(u, 1)
    setebc!(u, 2)
    applyebc!(u)
    numberdofs!(u)
    # display(u)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(1, 1),  (loc, conn, N) -> Area, false), CSys(2, 1), material)
    K = stiffness(femm,  geom,  u)

    fi = ForceIntensity(vec([0 -2000.0]))
    lfemm = FEMMBase(IntegDomain(FESetP1(reshape([3], 1,1)), PointRule()))
    F = distribloads(lfemm,  geom,  u,  fi,  3);
    fi = ForceIntensity(vec([+2000.0 0]))
    lfemm = FEMMBase(IntegDomain(FESetP1(reshape([5], 1,1)), PointRule()))
    F = F + distribloads(lfemm,  geom,  u,  fi,  3);
    fi = ForceIntensity(vec([+4000.0 +6000.0]))
    lfemm = FEMMBase(IntegDomain(FESetP1(reshape([6], 1,1)), PointRule()))
    F = F + distribloads(lfemm,  geom,  u,  fi,  3);

    K = cholesky(K)
    U=  K\F
    scattersysvec!(u, U[:])
    @test norm(u.values  - [ 0.0         0.0
                              0.0         0.0
                              0.0213333   0.0408366
                             -0.016       0.0461699
                              0.0426667   0.150091
                             -0.00533333  0.166091 ]) < 1.0e-4

    sfld =  elemfieldfromintegpoints(femm, geom, u, :Cauchy, 1)
    # display(sfld)
    # println("Cauchy = $(sfld.values)")
    @test norm(sfld.values - [5333.33; 3771.24; -4000.0; 1333.33; 5333.33; -5656.85; 2666.67; 4000.0]) < 1.0e-2
    vfld =  elemfieldfromintegpoints(femm, geom, u, :vm, 1)
    # display(vfld)

    File = "Planar_truss.vtk"
    MeshExportModule.VTK.vtkexportmesh(File, fens, fes;
    scalars=[("sx", sfld.values), ("vm", vfld.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

end
end
using .mbar2
mbar2.test()

module mmmLE10expiAbaqus2mmmm
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.MeshExportModule
using Test
import LinearAlgebra: norm, cholesky, cross
function test()

    # Thick elliptical plate with an elliptical hole is clamped on its exterior
    # boundary and is loaded with transverse  pressure.
    # This is a NAFEMS Benchmark, Test No. LE10.
    # The plate is discretized with solid elements.
    # Because of the symmetries of the geometry and load, only quarter of the plate is modeled.
    # The $\sigma_y=\sigma_2$ at the point $P$ is to be determined. Since the
    # target point is on the boundary of the domain it will not be an
    # integration node as we use Gauss quadrature. The reference value is -5.38 MPa.

    # println("LE10NAFEMS: Transverse deflection of elliptical plate with elliptical hole."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    qmagn = 1.0*phun("MEGA*PA");# transverse pressure
    sigma_yP = -5.38*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Ae =3.25*phun("m"); # Major radius of the exterior ellipse
    Be =2.75*phun("m"); # Minor radius of the exterior ellipse
    Ai =2.0*phun("m"); # Major radius of the interior ellipse
    Bi =1.0*phun("m"); # Minor radius of the interior ellipse
    Thickness = 0.6*phun("m")
    nc = 6; # number of elements per side
    nr = 5; # number of elements per side
    nt = 2; # number of elements through the thickness
    # nc = 26; # number of elements per side
    # nr = 25; # number of elements per side
    # nt = 18; # number of elements through the thickness
    tolerance = Thickness/nt/1000.; # Geometrical tolerance

    fens,fes = Q4block(1.0, pi/2, nr, nc)
    #
    @test  nt % 2 == 0
    fens,fes  = H8extrudeQ4(fens, fes,
      nt, (xyz, layer)->[xyz[1], xyz[2], (layer)/nt*Thickness]);

    # Select the  boundary faces, on the boundary that is clamped,  and on the part
    # of the boundary that is loaded with the transverse pressure
    bdryfes = meshboundary(fes);
    exteriorbfl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
    topbfl = selectelem(fens, bdryfes, box=[0.0, 1.0, 0.0, pi/2, Thickness, Thickness], inflate=tolerance);

    # Reshape the generated block into the elliptical plate
    for i=1:count(fens)
        r=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
        fens.xyz[i,:]=[(r*Ae+(1-r)*Ai)*cos(a) (r*Be+(1-r)*Bi)*sin(a) z];
    end

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

    l12 =connectednodes(subset(bdryfes, exteriorbfl)) # external boundary
    setebc!(u, l12, true, 1, 0.0)
    setebc!(u, l12, true, 2, 0.0)
    ll = selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness/2.0, Thickness/2.0], inflate = tolerance)
    l3 = intersect(ll, connectednodes(subset(bdryfes, exteriorbfl)))
    setebc!(u, l3, true, 3, 0.0)
    l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0) # symmetry plane X = 0
    l2 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l2,true, 2, 0.0) # symmetry plane Y = 0

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,topbfl), GaussRule(2, 2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        forceout .=  [0.0, 0.0, -qmagn]
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    material = MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    # The geometry field now needs to be associated with the FEMM
    femm = associategeometry!(femm, geom)

    K = stiffness(femm, geom, u)
    K = cholesky(K)
    U = K\(F2)
    scattersysvec!(u, U[:])

    nl = selectnode(fens, box=[Ai,Ai,0,0,Thickness,Thickness],inflate=tolerance);
    thecorneru = zeros(FFlt,1,3)
    gathervalues_asmat!(u, thecorneru, nl);
    thecorneru = thecorneru/phun("mm")
    # println("displacement =$(thecorneru) [MM] as compared to reference [-0.030939, 0, -0.10488] [MM]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extraptrend)#
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yP = $(sigma_yP/phun("MPa")) [MPa]")

    # println("$((nc, nr, nt)), $(fld.values[nl,1][1]/phun("MPa"))")
# println("$(fld.values[nl,1][1]/phun("MPa"))")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -4.627214556813842) < 1.0e-3

    # File =  "LE10NAFEMS_sigmay.vtk"
    # vtkexportmesh(File, fes.conn, geom.values,
    #                FinEtools.MeshExportModule.VTK.H8; vectors=[("u", u.values)],
    #                scalars=[("sigmay", fld.values)])
    # @async run(`"paraview.exe" $File`)
    # true

    AE = AbaqusExporter("LE10NAFEMS_H8");
    HEADING(AE, "LE10NAFEMS: Transverse deflection of elliptical plate with elliptical hole.");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(femm.integdomain.fes))
    ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(femm.integdomain.fes), connasarray(el1femm.integdomain.fes))
    NSET_NSET(AE, "l1", l1)
    NSET_NSET(AE, "l2", l2)
    NSET_NSET(AE, "l3", l3)
    NSET_NSET(AE, "l12", l12)
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    SURFACE_SECTION(AE, "TractionElements")
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1", u.is_fixed, u.fixed_values)
    DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec([0.0, 0.0, -qmagn]))
    END_STEP(AE)
    close(AE)
    lines = read(AE.filename)
    @test length(lines) - 10270 == 0
    try rm(AE.filename) catch end

end
end
using .mmmLE10expiAbaqus2mmmm
mmmLE10expiAbaqus2mmmm.test()

module mplate_w_hole_RECT_MSH8m
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.MeshExportModule
# using DataFrames
# using CSV
using Test
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()
    E = 210000*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Ri= 0.1*phun("M"); # hole radius
    Re = 2*Ri; # outer radius
    H = 0.01*phun("M") # thickness of the plate
    nRadial, nCircumferential=6, 3;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigyderrs = Dict{Symbol, FFltVec}()

    nelems = []
    for extrapolation in [:extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nelems = []
        for ref in [1]
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8elliphole(Ri, Ri, Re, Re, Thickness,
                2^ref*nCircumferential, 2^ref*nCircumferential, 2^ref*nRadial, 1)
            @test count(fes) == 144
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.VTK.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            bdryfes = meshboundary(fes);
            ixl = selectelem(fens, bdryfes, plane=[1.0, 0.0, 0.0, Re], thickness=tolerance);
            elxfemm =  FEMMBase(IntegDomain(subset(bdryfes,ixl), GaussRule(2, 2)))
            function pfunx(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxx(XYZ)
                forceout[2] = sigmaxy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfunx);
            Fx = distribloads(elxfemm, geom, u, fi, 2);
            iyl = selectelem(fens, bdryfes, plane=[0.0, 1.0, 0.0, Re], thickness=tolerance);
            elyfemm =  FEMMBase(IntegDomain(subset(bdryfes,iyl), GaussRule(2, 2)))
            function pfuny(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = -sigmaxy(XYZ)
                forceout[2] = sigmayy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfuny);
            Fy = distribloads(elyfemm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(Fx + Fy)
            scattersysvec!(u, U[:])

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            # println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            # println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            @test abs(sigyA/phun("MPa") - -0.8521990950600441) < 1.0e-3
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            # println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            @test abs(sigxB/phun("MPa") - 2.7749827820003374) < 1.0e-3
            # println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            # push!(nelems, count(fes))
            # push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.VTK.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA")),
            # ("sigmay", sigy.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    # df = DataFrame(nelems=vec(nelems),
    #     sigyderrtrend=vec(sigyderrs[:extraptrend]),
    #     sigyderrdefault=vec(sigyderrs[:extrapmean]))
    # File = "LE1NAFEMS_MSH8_convergence.CSV"
    # CSV.write(File, df)
    # @async run(`"paraview.exe" $File`)

end
end
using .mplate_w_hole_RECT_MSH8m
mplate_w_hole_RECT_MSH8m.test()

module mplate_w_hole_RECT_H20m
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.MeshExportModule
using FinEtools.MeshImportModule: import_ABAQUS
# using DataFrames
# using CSV
using Test
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()
    E = 210000*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Ri= 0.15*phun("M"); # hole radius
    Re = 2*Ri; # outer radius
    H = 0.01*phun("M") # thickness of the plate
    nRadial, nCircumferential=6, 3;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigyderrs = Dict{Symbol, FFltVec}()

    nelems = []
    for extrapolation in [:extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nelems = []
        for ref in [1]
            Thickness = H
            # Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8elliphole(Ri, Ri, Re, Re, Thickness,
            2^ref*nCircumferential, 2^ref*nCircumferential, 2^ref*nRadial, 1)
            fens,fes = H8toH20(fens,fes)
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.VTK.H20)
            # @async run(`"paraview.cexe" $File`)

            # println("My mesh=>$((count(fens), count(fes)))")
            @test count(fens) == 1131
            @test count(fes) == 144
            #
            # output = import_ABAQUS("plane_w_hole_m_debug.inp")
            # fens1,fes1 = output["fens"], output["fesets"][1]
            # println("Matlab mesh=>$((count(fens1), count(fes1[1])))")
            #
            #  fens3, newfes1, fes2 = mergemeshes(fens,fes, fens1,fes1[1], tolerance)
            #  fes3 = cat(2, newfes1)
            #  println("Merged mesh=>$((count(fens3), count(fes3)))")

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            bdryfes = meshboundary(fes);
            # ixl = selectelem(fens, bdryfes, plane=[1.0, 0.0, 0.0, Re], thickness=tolerance);
            ixl = selectelem(fens, bdryfes, box=[Re, Re, -Inf, +Inf, -Inf, +Inf], inflate = tolerance);
            elxfemm =  FEMMBase(IntegDomain(subset(bdryfes,ixl), GaussRule(2, 2)))
            function pfunx(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxx(XYZ)
                forceout[2] = sigmaxy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfunx);
            Fx = distribloads(elxfemm, geom, u, fi, 2);
            # iyl = selectelem(fens, bdryfes, plane=[0.0, 1.0, 0.0, Re], thickness=tolerance);
            iyl = selectelem(fens, bdryfes, box=[-Inf, +Inf, Re, Re, -Inf, +Inf], inflate = tolerance);
            elyfemm =  FEMMBase(IntegDomain(subset(bdryfes,iyl), GaussRule(2, 2)))
            function pfuny(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxy(XYZ)
                forceout[2] = sigmayy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfuny);
            Fy = distribloads(elyfemm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(Fx + Fy)
            scattersysvec!(u, U[:])
            # println("oof load = $(norm(Fx + Fy, 2))")
            @test abs(norm(Fx + Fy, 2) - 883.437848042617) < 1.0e-2

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, 00.0], inflate=tolerance);
            pointu = zeros(FFlt,length(nlA),3)
            gathervalues_asmat!(u, pointu, nlA);
            # println("disp@A = $(pointu/phun("mm")) [MM]")
            @test norm(pointu/phun("mm") - [0.00213238 0.0 0.0]) < 1.0e-4
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, 0.0], inflate=tolerance);
            pointu = zeros(FFlt,length(nlB),3)
            gathervalues_asmat!(u, pointu, nlB);
            # println("disp@B = $(pointu/phun("mm")) [MM]")
            @test norm(pointu/phun("mm") - [0.0 -0.000708141 0.0]) < 1.0e-4
            nlC = selectnode(fens, box=[Re, Re, Re, Re, Thickness, Thickness], inflate=tolerance);
            pointu = zeros(FFlt,length(nlC),3)
            gathervalues_asmat!(u, pointu, nlC);
            # println("disp@C = $(pointu/phun("mm")) [MM]")
            @test norm(pointu/phun("mm") - [0.00168556 -0.000455007 -1.4286e-5]) < 1.0e-4

            nlAallz = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlBallz = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlAallz,1], dims = 1)[1]
            sigyAtrue = sigmayy([Ri, 0.0, 0.0])
            # println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            @test abs(sigyA/phun("MPa") - -0.8513053526935438)/(sigyAtrue/phun("MPa")) < 1.0e-4
            sigxB = mean(sigx.values[nlBallz,1], dims = 1)[1]
            sigxBtrue = sigmaxx([0.0, Ri, 0.0])
            # println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            # @test abs(sigxB/phun("MPa") - 2.789413093796375)/3.0 < 1.0e-4
            # println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            # push!(nelems, count(fes))
            # push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            File =  "a.vtk"
            vtkexportmesh(File, connasarray(fes), geom.values,
                FinEtools.MeshExportModule.VTK.H20; vectors=[("u", u.values)],
                scalars=[("sigmax", sigx.values/phun("MEGA*PA")),
                ("sigmay", sigy.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
            try rm(File); catch end
        end
    end

    # df = DataFrame(nelems=vec(nelems),
    #     sigyderrtrend=vec(sigyderrs[:extraptrend]),
    #     sigyderrdefault=vec(sigyderrs[:extrapmean]))
    # File = "LE1NAFEMS_MSH8_convergence.CSV"
    # CSV.write(File, df)
    # @async run(`"paraview.exe" $File`)

end
end
using .mplate_w_hole_RECT_H20m
mplate_w_hole_RECT_H20m.test()

module mplate_w_hole_MST10m
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.MeshExportModule
using Test
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrapolation in [:extrapmean] # :extraptrend
        sigxderrs[extrapolation] = FFltVec[]
        sigyderrs[extrapolation] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 1:1
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.VTK.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), SimplexRule(2, 3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, SimplexRule(3, 4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            # println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
                nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            # println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            @test  abs(sigyA/phun("MPa") - -0.6705333234697736)/(sigyAtrue/phun("MPa")) < 1.0e-4
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            # println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            @test  abs(sigxB/phun("MPa") - 2.301542874107758)/3.0 < 1.0e-4
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrapolation], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrapolation], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.VTK.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

end
end
using .mplate_w_hole_MST10m
mplate_w_hole_MST10m.test()


text:
pwd()

a = 2 + 2;

b = 2 * 38;

a

b

c = a + b;

c

ext: jl

max_stars_repo_name: alexandreloures/econometrics-with-Julia

max_issues_repo_name: alexandreloures/econometrics-with-Julia

max_forks_repo_name: alexandreloures/econometrics-with-Julia


text:
import OcalAPI: check_keys, check_and_parse_values!, check_value_sizes,
                parse_parameters!

@testset "parse parameters" begin
    @testset "missing keys" begin
        p = deepcopy(RAW_PARAMETERS)
        delete!(p, "data")
        delete!(p["params"], "C")
        @test_throws ArgumentError check_keys(p)
    end

    @testset "check and parse" begin

        @testset "data" begin
            p = deepcopy(RAW_PARAMETERS)
            p["data"] = trues(2, 5)
            @test_throws ArgumentError check_and_parse_values!(p)
            p = deepcopy(RAW_PARAMETERS)
            p["data"] = [[2, 5], [5, 2, 10]]
            @test_throws ArgumentError check_and_parse_values!(p)
        end

        @testset "labels" begin
            p = deepcopy(RAW_PARAMETERS)
            p["labels"] = nothing
            @test_throws ArgumentError check_and_parse_values!(p)
        end

        @testset "params" begin
            p = deepcopy(RAW_PARAMETERS)
            p["params"] = nothing
            @test_throws ArgumentError check_and_parse_values!(p)
        end

        @testset "query_history" begin
            p = deepcopy(RAW_PARAMETERS)
            p["query_history"] = [["test"]]
            @test_throws ArgumentError check_and_parse_values!(p)
            p = deepcopy(RAW_PARAMETERS)
            p["query_history"] = [1, 2]
            @test_throws ArgumentError check_and_parse_values!(p)
        end

        @testset "subspaces" begin
            p = deepcopy(RAW_PARAMETERS)
            p["subspaces"] = [["test"]]
            @test_throws ArgumentError check_and_parse_values!(p)
            p = deepcopy(RAW_PARAMETERS)
            p["subspaces"] = [1, 2]
            @test_throws ArgumentError check_and_parse_values!(p)
        end

        @testset "subspace_grids" begin
            p = deepcopy(RAW_PARAMETERS)
            p["subspace_grids"] = []
            @test check_and_parse_values!(p) == nothing
            p = deepcopy(RAW_PARAMETERS)
            p["subspace_grids"] = p["subspace_grids"][1]
            @test_throws ArgumentError check_and_parse_values!(p)
        end
    end

    @testset "check value sizes" begin
        p = deepcopy(TEST_PARAMETERS)
        p["subspace_grids"] = []
        @test_throws ArgumentError check_value_sizes(p)
        p = deepcopy(TEST_PARAMETERS)
        p["labels"] = []
        @test_throws ArgumentError check_value_sizes(p)
    end

    @testset "encode decode" begin
        binary_parameters = transcode(GzipCompressor, JSON.json(TEST_PARAMETERS))
        r = JSON.parse(String(transcode(GzipDecompressor, binary_parameters)))
        @test parse_parameters!(r) == nothing
        @test all(collect(values(r)) .== collect(values(TEST_PARAMETERS)))
    end
end


text:
using ArrayViews
using StatsBase
using Base.Test

@test isa(weights([1, 2, 3]), WeightVec{Int})
@test isa(weights([1., 2., 3.]), WeightVec{Float64})

@test isa(weights([1 2 3; 4 5 6]), WeightVec{Int})

@test isempty(weights(Float64[]))

w  = [1., 2., 3.]
wv = weights(w)
@test eltype(wv) === Float64
@test length(wv) === 3
@test values(wv) === w
@test sum(wv) === 6.0
@test !isempty(wv)

b  = trues(3)
bv = weights(b)
@test eltype(bv) === Bool
@test length(bv) === 3
@test values(bv) === b
@test sum(bv)    === 3
@test !isempty(bv)

## wsum

x = [6., 8., 9.]
w = [2., 3., 4.]
p = [1. 2. ; 3. 4.]
q = [1., 2., 3., 4.]

@test wsum(Float64[], Float64[]) === 0.0
@test wsum(x, w) === 72.0
@test wsum(p, q) === 29.0

## wsum along dimensions

@test wsum(x, w, 1) == [72.0]

x  = rand(6, 8)
w1 = rand(6)
w2 = rand(8)

@test size(wsum(x, w1, 1)) == (1, 8)
@test size(wsum(x, w2, 2)) == (6, 1)

@test_approx_eq wsum(x, w1, 1) sum(x .* w1, 1)
@test_approx_eq wsum(x, w2, 2) sum(x .* w2', 2)

x = rand(6, 5, 4)
w1 = rand(6)
w2 = rand(5)
w3 = rand(4)

@test size(wsum(x, w1, 1)) == (1, 5, 4)
@test size(wsum(x, w2, 2)) == (6, 1, 4)
@test size(wsum(x, w3, 3)) == (6, 5, 1)

@test_approx_eq wsum(x, w1, 1) sum(x .* w1, 1)
@test_approx_eq wsum(x, w2, 2) sum(x .* w2', 2)
@test_approx_eq wsum(x, w3, 3) sum(x .* reshape(w3, 1, 1, 4), 3)

v = view(x, 2:4, :, :)

@test_approx_eq wsum(v, w1[1:3], 1) sum(v .* w1[1:3], 1)
@test_approx_eq wsum(v, w2, 2) sum(v .* w2', 2)
@test_approx_eq wsum(v, w3, 3) sum(v .* reshape(w3, 1, 1, 4), 3)

## wsum for Arrays with non-BlasReal elements

x = rand(1:100, 6, 8)
w1 = rand(6)
w2 = rand(8)

@test_approx_eq wsum(x, w1, 1) sum(x .* w1, 1)
@test_approx_eq wsum(x, w2, 2) sum(x .* w2', 2)

## wsum!

x = rand(6)
w = rand(6)

r = ones(1)
@test wsum!(r, x, w, 1; init=true) === r
@test_approx_eq r [dot(x, w)]

r = ones(1)
@test wsum!(r, x, w, 1; init=false) === r
@test_approx_eq r [dot(x, w) + 1.0]

x = rand(6, 8)
w1 = rand(6)
w2 = rand(8)

r = ones(1, 8)
@test wsum!(r, x, w1, 1; init=true) === r
@test_approx_eq r sum(x .* w1, 1)

r = ones(1, 8)
@test wsum!(r, x, w1, 1; init=false) === r
@test_approx_eq r sum(x .* w1, 1) .+ 1.0

r = ones(6)
@test wsum!(r, x, w2, 2; init=true) === r
@test_approx_eq r sum(x .* w2', 2)

r = ones(6)
@test wsum!(r, x, w2, 2; init=false) === r
@test_approx_eq r sum(x .* w2', 2) .+ 1.0

x = rand(8, 6, 5)
w1 = rand(8)
w2 = rand(6)
w3 = rand(5)

r = ones(1, 6, 5)
@test wsum!(r, x, w1, 1; init=true) === r
@test_approx_eq r sum(x .* w1, 1)

r = ones(1, 6, 5)
@test wsum!(r, x, w1, 1; init=false) === r
@test_approx_eq r sum(x .* w1, 1) .+ 1.0

r = ones(8, 1, 5)
@test wsum!(r, x, w2, 2; init=true) === r
@test_approx_eq r sum(x .* w2', 2)

r = ones(8, 1, 5)
@test wsum!(r, x, w2, 2; init=false) === r
@test_approx_eq r sum(x .* w2', 2) .+ 1.0

r = ones(8, 6)
@test wsum!(r, x, w3, 3; init=true) === r
@test_approx_eq r sum(x .* reshape(w3, (1, 1, 5)), 3)

r = ones(8, 6)
@test wsum!(r, x, w3, 3; init=false) === r
@test_approx_eq r sum(x .* reshape(w3, (1, 1, 5)), 3) .+ 1.0

## the sum and mean syntax

@test_approx_eq sum([1.0, 2.0, 3.0], weights([1.0, 0.5, 0.5])) 3.5
@test_approx_eq sum(1:3, weights([1.0, 1.0, 0.5])) 4.5

@test_approx_eq mean([1:3;], weights([1.0, 1.0, 0.5])) 1.8
@test_approx_eq mean(1:3, weights([1.0, 1.0, 0.5])) 1.8

a = reshape(1.0:27.0, 3, 3, 3)
for wt in ([1.0, 1.0, 1.0], [1.0, 0.2, 0.0], [0.2, 0.0, 1.0])
    @test_approx_eq sum(a, weights(wt), 1) sum(a.*reshape(wt, length(wt), 1, 1), 1)
    @test_approx_eq sum(a, weights(wt), 2) sum(a.*reshape(wt, 1, length(wt), 1), 2)
    @test_approx_eq sum(a, weights(wt), 3) sum(a.*reshape(wt, 1, 1, length(wt)), 3)
    @test_approx_eq mean(a, weights(wt), 1) sum(a.*reshape(wt, length(wt), 1, 1), 1)/sum(wt)
    @test_approx_eq mean(a, weights(wt), 2) sum(a.*reshape(wt, 1, length(wt), 1), 2)/sum(wt)
    @test_approx_eq mean(a, weights(wt), 3) sum(a.*reshape(wt, 1, 1, length(wt)), 3)/sum(wt)
    @test_throws ErrorException mean(a, weights(wt), 4)
end

# Weighted median tests
data = (
    [7, 1, 2, 4, 10],
    [7, 1, 2, 4, 10],
    [7, 1, 2, 4, 10, 15],
    [1, 2, 4, 7, 10, 15],
    [0, 10, 20, 30],
    [1, 2, 3, 4, 5],
    [1, 2, 3, 4, 5],
    [30, 40, 50, 60, 35],
    [2, 0.6, 1.3, 0.3, 0.3, 1.7, 0.7, 1.7, 0.4],
    [3.7, 3.3, 3.5, 2.8],
    [100, 125, 123, 60, 45, 56, 66],
    [2, 2, 2, 2, 2, 2],
    [2.3],
    [-2, -3, 1, 2, -10],
    [1, 2, 3, 4, 5],
    [5, 4, 3, 2, 1],
    [-2, 2, -1, 3, 6],
    [-10, 1, 1, -10, -10],
    [2, 4],
    [2, 2, 4, 4],
    [2, 2, 2, 4]
)
wt = (
    [1, 1/3, 1/3, 1/3, 1],
    [1, 1, 1, 1, 1],
    [1, 1/3, 1/3, 1/3, 1, 1],
    [1/3, 1/3, 1/3, 1, 1, 1],
    [30, 191, 9, 0],
    [10, 1, 1, 1, 9],
    [10, 1, 1, 1, 900],
    [1, 3, 5, 4, 2],
    [2, 2, 0, 1, 2, 2, 1, 6, 0],
    [5, 5, 4, 1],
    [30, 56, 144, 24, 55, 43, 67],
    [0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
    [12],
    [7, 1, 1, 1, 6],
    [1, 0, 0, 0, 2],
    [1, 2, -3, 4, -5],
    [0.1, 0.2, 0.3, -0.2, 0.1],
    [-1, -1, -1, -1, 1],
    [1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1]
)
median_answers = (7.0,   4.0,   8.5,
                  8.5,  10.0,   2.5,
                  5.0,  50.0,   1.7,
                  3.5, 100.0,   2.0,
                  2.3,  -2.0,   5.0,
                  2.0,  -1.0, -10.0,
                  3.0,   3.0,   2.0)
num_tests = length(data)
for i = 1:num_tests
    @test wmedian(data[i], wt[i]) == median_answers[i]
    @test wmedian(data[i], weights(wt[i])) == median_answers[i]
    @test median(data[i], weights(wt[i])) == median_answers[i]
    for j = 1:100
        # Make sure the weighted median does not change if the data
        # and weights are reordered.
        reorder = sortperm(rand(length(data[i])))
        @test median(data[i][reorder], weights(wt[i][reorder])) == median_answers[i]
    end
end
data = [4, 3, 2, 1]
wt = [0, 0, 0, 0]
@test_throws MethodError wmedian(data[1])
@test_throws ErrorException median(data, weights(wt))
@test_throws ErrorException wmedian(data, wt)
@test_throws ErrorException median((Float64)[], weights((Float64)[]))
wt = [1, 2, 3, 4, 5]
@test_throws ErrorException median(data, weights(wt))
@test_throws MethodError median([4 3 2 1 0], weights(wt))
@test_throws MethodError median([[1 2];[4 5];[7 8];[10 11];[13 14]], weights(wt))
data = [1, 3, 2, NaN, 2]
@test isnan(median(data, weights(wt)))
wt = [1, 2, NaN, 4, 5]
@test_throws ErrorException median(data, weights(wt))
data = [1, 3, 2, 1, 2]
@test_throws ErrorException median(data, weights(wt))
wt = [-1, -1, -1, -1, -1]
@test_throws ErrorException median(data, weights(wt))
wt = [-1, -1, -1, 0, 0]
@test_throws ErrorException median(data, weights(wt))

# Weighted quantile tests
data = (
    [7, 1, 2, 4, 10],
    [7, 1, 2, 4, 10],
    [7, 1, 2, 4, 10, 15],
    [1, 2, 4, 7, 10, 15],
    [0, 10, 20, 30],
    [1, 2, 3, 4, 5],
    [1, 2, 3, 4, 5],
    [30, 40, 50, 60, 35],
    [2, 0.6, 1.3, 0.3, 0.3, 1.7, 0.7, 1.7],
    [1, 2, 2],
    [3.7, 3.3, 3.5, 2.8],
    [100, 125, 123, 60, 45, 56, 66],
    [2, 2, 2, 2, 2, 2],
    [2.3],
    [-2, -3, 1, 2, -10],
    [1, 2, 3, 4, 5],
    [5, 4, 3, 2, 1],
    [-2, 2, -1, 3, 6],
    [-10, 1, 1, -10, -10],
)
wt = (
    weights([1, 1/3, 1/3, 1/3, 1]),
    weights([1, 1, 1, 1, 1]),
    weights([1, 1/3, 1/3, 1/3, 1, 1]),
    weights([1/3, 1/3, 1/3, 1, 1, 1]),
    weights([30, 191, 9, 0]),
    weights([10, 1, 1, 1, 9]),
    weights([10, 1, 1, 1, 900]),
    weights([1, 3, 5, 4, 2]),
    weights([2, 2, 5, 1, 2, 2, 1, 6]),
    weights([0.1, 0.1, 0.8]),
    weights([5, 5, 4, 1]),
    weights([30, 56, 144, 24, 55, 43, 67]),
    weights([0.1, 0.2, 0.3, 0.4, 0.5, 0.6]),
    weights([12]),
    weights([7, 1, 1, 1, 6]),
    weights([1, 0, 0, 0, 2]),
    weights([1, 2, 3, 4, 5]),
    weights([0.1, 0.2, 0.3, 0.2, 0.1]),
    weights([1, 1, 1, 1, 1]),
)
quantile_answers = (
    [1.0,3.6000000000000005,6.181818181818182,8.2,10.0],
    [1.0,2.0,4.0,7.0,10.0],
    [1.0,4.75,8.0,10.833333333333334,15.0],
    [1.0,4.75,8.0,10.833333333333334,15.0],
    [0.0,6.1387900355871885,11.600000000000001,15.912500000000001,30.0],
    [1.0,1.5365853658536586,2.5999999999999996,4.405405405405405,5.0],
    [1.0,4.239377950569287,4.492918633712858,4.746459316856429,5.0],
    [30.0,38.75,45.714285714285715,52.85714285714286,60.0],
    [0.3,0.6903846153846154,1.484,1.7,2.0],
    [1.0,2.0,2.0,2.0,2.0],
    [2.8,3.3361111111111112,3.4611111111111112,3.581578947368421,3.7],
    [45.0,59.88593155893536,100.08846153846153,118.62115384615385,125.0],
    [2.0,2.0,2.0,2.0,2.0],
    [2.3,2.3,2.3,2.3,2.3],
    [-10.0,-5.52,-2.5882352941176467,-0.9411764705882351,2.0],
    [1.0,1.75,4.25,4.625,5.0],
    [1.0,1.625,2.3333333333333335,3.25,5.0],
    [-2.0,-0.5384615384615388,1.5384615384615383,2.6999999999999997,6.0],
    [-10.0,-10.0,-10.0,1.0,1.0]
)
p = [0.0, 0.25, 0.5, 0.75, 1.0]

srand(10)
for i = 1:length(data)
    @test_approx_eq quantile(data[i], wt[i], p) quantile_answers[i]
    for j = 1:10
        # order of p does not matter
        reorder = sortperm(rand(length(p)))
        @test_approx_eq quantile(data[i], wt[i], p[reorder]) quantile_answers[i][reorder]
    end
    for j = 1:10
        # order of w does not matter
        reorder = sortperm(rand(length(data[i])))
        @test_approx_eq quantile(data[i][reorder], weights(wt[i][reorder]), p) quantile_answers[i]
    end
end
# w = 1 corresponds to base quantile
for i = 1:length(data)
    @test_approx_eq quantile(data[i], weights(ones(Int64, length(data[i]))), p) quantile(data[i], p)
    for j = 1:10
        prandom = rand(4)
        @test_approx_eq quantile(data[i], weights(ones(Int64, length(data[i]))),  prandom) quantile(data[i], prandom)
    end
end

# other syntaxes
v = [7, 1, 2, 4, 10]
w = [1, 1/3, 1/3, 1/3, 1]
answer = 6.181818181818182
@test_approx_eq  quantile(data[1], weights(w), 0.5) answer
@test_approx_eq  wquantile(data[1], weights(w), [0.5]) [answer]
@test_approx_eq  wquantile(data[1], weights(w), 0.5) answer
@test_approx_eq  wquantile(data[1], w, [0.5]) [answer]
@test_approx_eq  wquantile(data[1], w, 0.5)  answer


text:
module ToyFHE

export encrypt, decrypt, keygen, keyswitch, PolyCRTEncoding, CKKSEncoding
export BFVParams, BGVParams, CKKSParams
export PubKey, PrivKey, EvalKey, EvalMultKey, GaloisKey, GaloisKeys, KeyPair, CipherText
export modswitch, modswitch_drop
export ℛ_cipher, ℛ_key, ℛ_plain

export coefftype

using Random
using Distributions
using Nemo
using Hecke
using GaloisFields
using Mods
using BitIntegers
using StructArrays

include("poly.jl")
include("signedmod.jl")
include("pow2_cyc_rings.jl")
using .NTT
export NegacyclicRing
include("utils.jl")
using .Utils: @fields_as_locals
include("rlwe_she.jl")
include("cryptparams.jl")
include("crt.jl")
include("bgv.jl")
include("bfv.jl")
include("ckks.jl")
include("ckksencoding.jl")
include("nemo.jl")
include("encoding.jl")
include("polycrtencoding.jl")
include("insecuredebug.jl")
include("modulusraising.jl")

import .Utils: plaintext_space
export plaintext_space

if VERSION < v"1.4.0-DEV.208"
    include("div_hacks.jl")
end

end


text:
mutable struct CrossValidationResult <: ModelSelection.ModelSelectionResult

    k::Int64
    s::Float64

    ttest
    datanames
    average_data
    median_data
    data

    function CrossValidationResult(
            k,
            s,
            ttest,
            datanames,
            average_data,
            median_data,
            data
        )
        new(k, s, ttest, datanames, average_data, median_data, data)
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct VirtualInertia <: ActivePowerControl
        Ta::Float64
        kd::Float64
        kω::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
    end

Parameters of a Virtual Inertia with SRF using VSM for active power controller

# Arguments
- `Ta::Float64`: VSM inertia constant, validation range: (0, nothing)
- `kd::Float64`: VSM damping constant, validation range: (0, nothing)
- `kω::Float64`: frequency droop gain, validation range: (0, nothing)
- `P_ref::Float64`: Reference Power Set-point, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
"""
mutable struct VirtualInertia <: ActivePowerControl
    "VSM inertia constant"
    Ta::Float64
    "VSM damping constant"
    kd::Float64
    "frequency droop gain"
    kω::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
end

function VirtualInertia(Ta, kd, kω, P_ref=1.0, ext=Dict{String, Any}(), )
    VirtualInertia(Ta, kd, kω, P_ref, ext, [:ω_oc, :θ_oc], 2, )
end

function VirtualInertia(; Ta, kd, kω, P_ref=1.0, ext=Dict{String, Any}(), )
    VirtualInertia(Ta, kd, kω, P_ref, ext, )
end

# Constructor for demo purposes; non-functional.
function VirtualInertia(::Nothing)
    VirtualInertia(;
        Ta=0,
        kd=0,
        kω=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get VirtualInertia Ta."""
get_Ta(value::VirtualInertia) = value.Ta
"""Get VirtualInertia kd."""
get_kd(value::VirtualInertia) = value.kd
"""Get VirtualInertia kω."""
get_kω(value::VirtualInertia) = value.kω
"""Get VirtualInertia P_ref."""
get_P_ref(value::VirtualInertia) = value.P_ref
"""Get VirtualInertia ext."""
get_ext(value::VirtualInertia) = value.ext
"""Get VirtualInertia states."""
get_states(value::VirtualInertia) = value.states
"""Get VirtualInertia n_states."""
get_n_states(value::VirtualInertia) = value.n_states

"""Set VirtualInertia Ta."""
set_Ta!(value::VirtualInertia, val::Float64) = value.Ta = val
"""Set VirtualInertia kd."""
set_kd!(value::VirtualInertia, val::Float64) = value.kd = val
"""Set VirtualInertia kω."""
set_kω!(value::VirtualInertia, val::Float64) = value.kω = val
"""Set VirtualInertia P_ref."""
set_P_ref!(value::VirtualInertia, val::Float64) = value.P_ref = val
"""Set VirtualInertia ext."""
set_ext!(value::VirtualInertia, val::Dict{String, Any}) = value.ext = val
"""Set VirtualInertia states."""
set_states!(value::VirtualInertia, val::Vector{Symbol}) = value.states = val
"""Set VirtualInertia n_states."""
set_n_states!(value::VirtualInertia, val::Int64) = value.n_states = val


text:
export CollectGems

mutable struct CollectGems{R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Gem}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    num_gem_init::Int
    num_gem_current::Int
    gem_reward::Float64
    reward::Float64
    rng::R
end

function CollectGems(;n=8, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT, rng=Random.GLOBAL_RNG)
    objects = (EMPTY, WALL, GEM)
    w = GridWorldBase(objects, n, n)

    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true

    num_gem_init = n - 1
    num_gem_current = num_gem_init

    gem_reward = 1.0
    reward = 0.0

    env = CollectGems(w, agent_start_pos, Agent(dir=agent_start_dir), num_gem_init, num_gem_current, gem_reward, reward, rng)

    reset!(env)

    return env

end

function (w::CollectGems)(::MoveForward)
    dir = get_dir(w.agent)
    dest = dir(w.agent_pos)
    w.reward = 0.0
    if !w.world[WALL, dest]
        w.agent_pos = dest
        if w.world[GEM, dest]
            w.world[GEM, dest] = false
            w.world[EMPTY, dest] = true
            w.num_gem_current = w.num_gem_current - 1
            w.reward = w.gem_reward
        end
    end
    w
end

function (w::CollectGems)(action::Union{TurnRight, TurnLeft})
    w.reward = 0.0
    agent = get_agent(w)
    set_dir!(agent, action(get_dir(agent)))
    w
end

RLBase.get_terminal(w::CollectGems) = w.num_gem_current <= 0

RLBase.get_reward(w::CollectGems) = w.reward

function RLBase.reset!(w::CollectGems)

    n = size(w.world)[end]
    w.world[EMPTY, 2:n-1, 2:n-1] .= true
    w.world[GEM, 1:n, 1:n] .= false
    w.num_gem_current = w.num_gem_init

    w.reward = 0.0
    w.agent_pos = CartesianIndex(2, 2)
    w.agent.dir = RIGHT

    gem_placed = 0
    while gem_placed < w.num_gem_init
        gem_pos = CartesianIndex(rand(w.rng, 2:n-1), rand(w.rng, 2:n-1))
        if (gem_pos == w.agent_pos) || (w.world[GEM, gem_pos] == true)
            continue
        else
            w.world[GEM, gem_pos] = true
            w.world[EMPTY, gem_pos] = false
            gem_placed = gem_placed + 1
        end
    end

    return w
end


text:
using StanSample, MCMCChains, CSV

df = filter(row -> row[:age] >= 18,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "Howell1.csv"), DataFrame))

# Define the Stan language model

m4_3s = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] height; // Predictor
 vector[N] weight; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 height ~ normal(alpha + weight * beta , sigma);
}

generated quantities {
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m_4_3s = SampleModel("m4.3s", m4_3s);

# Input data for cmdstan

m4_3_data = Dict("N" => size(df, 1), "height" => df[!, :height],
  "weight" => df[!, :weight]);

# Sample using cmdstan

rc = stan_sample(m_4_3s, data=m4_3_data)

# Describe the draws

if success(rc)
  chn = read_samples(m_4_3s; output_format=:mcmcchains)
  #chn = set_names(chn, Dict("mu" => "μ", "sigma" => "σ"))
  describe(chn)
end


text:
function Renumbering(Grid)
for iF=1:Grid.NumFaces
  Grid.Faces[iF]=RenumberingFace4(Grid.Faces[iF],Grid);
end
for iE=1:Grid.NumEdges
  Grid.Edges[iE]=PosEdgeInFace(Grid.Edges[iE],Grid);
end
return Grid
end

function PosEdgeInFace(Edge,Grid)
Edge.FE=zeros(1,2);
for iF=1:size(Edge.F,2)
  F=Edge.F[iF];
  for iE=1:4
    if Edge.EI==Grid.Edges[Grid.Faces[F].E[iE]].EI
      Edge.FE[iF]=iE;
      break
    end
  end
end
if size(Edge.F,2)>1
  if Edge.FE[1] > Edge.FE[2]
    iTemp=Edge.FE[1];
    Edge.FE[1]=Edge.FE[2];
    Edge.FE[2]=iTemp;
    iTemp=Edge.F[1];
    Edge.F[1]=Edge.F[2];
    Edge.F[2]=iTemp;
  end
end
return Edge
end
function RenumberingFace4(Face,Grid)
local iN
for iN_in=1:4
  iN = iN_in
  N=Face.N[iN];
  num=0;
  for iE=1:4
    if N==Grid.Edges[Face.E[iE]].N[1]
      num=num+1;
    end
    if num==2
      break
    end
  end
  if num==2
    break
  end
end
if iN>1
  NTemp=[Face.N Face.N];
  ETemp=[Face.E Face.E];
  PTemp=[Face.P Face.P];
  for i=1:4
    Face.N[i]=NTemp[i+iN-1];
    Face.E[i]=ETemp[i+iN-1];
    Face.P[:,i]=PTemp[:,i+iN-1];
  end
end
OrientL=zeros(4,1);
OrientL[1]=-1;
OrientL[2]=1;
OrientL[3]=1;
OrientL[4]=-1;
Face.OrientE = Vector{Int}(undef, 0)
for i=1:4
  if Grid.Edges[Face.E[i]].N[1]==Face.N[i]
    push!(Face.OrientE, 1*OrientL[i]);
  else
    push!(Face.OrientE, -1*OrientL[i]);
  end
end
return Face
end


text:
#==============================================================================================
Here we keep some macros for doing boilerplate code generation...
==============================================================================================#
# TODO deal with page crossings!!!

opfuncname(opcode::UInt8) = Symbol(string("op", hexstring(opcode), "!"))
opstring(fname::Symbol) = strip(string(fname), '!')

function assemblydict(opname::Symbol, mode::Symbol, opcode::UInt8)
    :(ASSEMBLY_DICT[($(opstring(opname)),$(QuoteNode(mode)))] = $opcode)
end

function opdict(opcode::UInt8, nbytes::Int, ncycles::Int)
    :(OPCODES[$opcode+1] = ($(opfuncname(opcode)), $nbytes, $ncycles))
end

function opdef_Implicit(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu)
    end
end

function opdef_Immediate(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, opargs(cs, 1))
    end
end

function opdef_ZeroPage(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, Π8(opargs(cs, 1)))
    end
end

function opdef_ZeroPageX(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, DirectX, Π8(opargs(cs, 1)))
    end
end

function opdef_ZeroPageY(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, DirectY, Π8(opargs(cs, 1)))
    end
end

function opdef_Absolute(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, Direct, Π16(opargs(cs,1), opargs(cs,2)))
    end
end

function opdef_AbsoluteX(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, DirectX, Π16(opargs(cs,1), opargs(cs,2)))
    end
end

function opdef_AbsoluteY(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, DirectY, Π16(opargs(cs,1), opargs(cs,2)))
    end
end

function opdef_Indirect(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, Indirect, Π16(opargs(cs,1), opargs(cs,2)))
    end
end

function opdef_IndirectX(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, IndirectX, Π8(opargs(cs,1)))
    end
end

function opdef_IndirectY(opname::Symbol, opcode::UInt8)
    fname = opfuncname(opcode)
    quote
        $fname(cs::Chipset) = $opname(cs.cpu, cs.ram, IndirectY, Π8(opargs(cs,1)))
    end
end

function opdef(mode::Symbol, opname::Symbol, opcode::UInt8)
    if mode == :Implicit
        return opdef_Implicit(opname, opcode)
    elseif mode == :Immediate
        return opdef_Immediate(opname, opcode)
    elseif mode == :ZeroPage
        return opdef_ZeroPage(opname, opcode)
    elseif mode == :ZeroPageX
        return opdef_ZeroPageX(opname, opcode)
    elseif mode == :ZeroPageY
        return opdef_ZeroPageY(opname, opcode)
    elseif mode == :Absolute
        return opdef_Absolute(opname, opcode)
    elseif mode == :AbsoluteX
        return opdef_AbsoluteX(opname, opcode)
    elseif mode == :AbsoluteY
        return opdef_AbsoluteY(opname, opcode)
    elseif mode == :Indirect
        return opdef_Indirect(opname, opcode)
    elseif mode == :IndirectX
        return opdef_IndirectX(opname, opcode)
    elseif mode == :IndirectY
        return opdef_IndirectY(opname, opcode)
    elseif mode == :Relative
        return opdef_Relative(opname, opcode)
    else
        throw(ArgumentError("Invalid oparation mode $mode."))
    end
end

macro opdef(opname::Symbol, defblock::Expr)
    if @capture(defblock, begin defs__ end)
        defs = [ex.args for ex ∈ defs]
        funcdefs = [opdef(d[1], opname, d[2]) for d ∈ defs]
        dictries = [opdict(d[2], d[3], d[4]) for d ∈ defs]
        assemblies = [assemblydict(opname, d[1], d[2]) for d ∈ defs]
    else
        throw(AssertionError("Improper op definition for $opname."))
    end
    esc(quote
        $(funcdefs...)
        $(dictries...)
        $(assemblies...)
    end)
end


text:
module MbedTLS

using Random, Sockets, MbedTLS_jll
import Sockets: TCPSocket

export
# Message digests
    MD_NONE,
    MD_MD2,
    MD_MD4,
    MD_MD5,
    MD_SHA1,
    MD_SHA224,
    MD_SHA256,
    MD_SHA384,
    MD_SHA,
    digest,
    digest!,

# Symmetric encryption,
    encrypt,
    decrypt,
    Cipher,
    CIPHER_AES,
    CIPHER_DES,
    CIPHER_3DES,
    CIPHER_CAMELLIA,
    CIPHER_BLOWFISH,
    CIPHER_ARC4,

# Types
    CtrDrbg,
    RSA,
    SSLConfig,
    Entropy,
    CRT

import Base: show

const MBED_SUCCESS = Cint(0)

include("constants.jl")
include("error.jl")
include("md.jl")
include("cipher.jl")
include("rsa.jl")
include("entropy.jl")
include("ctr_drbg.jl")
include("pk.jl")
include("x509_crt.jl")
include("ssl.jl")

function __init__()
    __ctr_drbg__init__()
    __sslinit__()
    __entropyinit__()
end

tls_dbg(level, filename, number, msg) = @warn "MbedTLS emitted debug info: $msg in $filename:$number"

secrets_log = ""
secrets_log_file = ""
secrets_log_state = :incomplete

function tls_dbg_log_secrets(level, filename, number, msg)
    global secrets_log
    global secrets_log_file
    global secrets_log_state
    if secrets_log_state == :complete
        return
    elseif msg == "dumping 'client hello, random bytes' (32 bytes)\n"
        secrets_log = "CLIENT_RANDOM "
        secrets_log_state = :client_random
    elseif msg == "dumping 'master secret' (48 bytes)\n"
        secrets_log_state = :master_secret
        secrets_log *= " "
    elseif secrets_log_state in (:client_random, :master_secret)
        if occursin(r"^00[0-2]0:", msg)
            secrets_log *= join(split(msg)[2:17])
        else
            if secrets_log_state == :master_secret
                secrets_log *= "\n"
                open(secrets_log_file, append=true) do io
                    write(io, secrets_log)
                end
                secrets_log_state = :incomplete
                secrets_log = ""
            else
                secrets_log_state = :incomplete
            end
        end
    end
end

"""
    SSLConfig(cert_file, key_file)

Initialise an SSLConfig from a certificate and key file on disk.
This is probably the method you want if you are the server in an HTTPS connection or similar.

# Example
```
julia> MbedTLS.SSLConfig("self-signed-certificate.pem", "keyfile.pem")
MbedTLS.SSLConfig()
```
"""
# already defined SSLConfig and SSLContext types in ssl.jl
function SSLConfig(cert_file, key_file)
    ssl_cert = MbedTLS.crt_parse_file(cert_file)
    key = MbedTLS.parse_keyfile(key_file)
    conf = MbedTLS.SSLConfig()
    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.config_defaults!(conf, endpoint=MbedTLS.MBEDTLS_SSL_IS_SERVER)
    MbedTLS.seed!(rng, entropy)
    MbedTLS.rng!(conf, rng)
    MbedTLS.own_cert!(conf, ssl_cert, key)
    MbedTLS.dbg!(conf, tls_dbg)
    return conf
end

"""
    SSLConfig(verify::Bool; log_secrets=nothing)

Initialise a client SSLConfig for connecting to a server.

If `verify` is false, do not check that certificates are valid.

If `log_secrets` is a string, save connection secrets to a file with that name.
This is useful for decrypting traffic captured with Wireshark when debugging.
"""
function SSLConfig(verify::Bool; log_secrets=nothing)
    conf = MbedTLS.SSLConfig()
    MbedTLS.config_defaults!(conf)

    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)
    MbedTLS.rng!(conf, rng)

    MbedTLS.authmode!(conf,
      verify ? MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED : MbedTLS.MBEDTLS_SSL_VERIFY_NONE)
    if log_secrets !== nothing
        global secrets_log_file = log_secrets
        set_dbg_level(4)
        MbedTLS.dbg!(conf, tls_dbg_log_secrets)
    else
        MbedTLS.dbg!(conf, tls_dbg)
    end
    MbedTLS.ca_chain!(conf)
    conf
end

# already defined setup! in ssl.jl
associate!(tls::SSLContext, tcp::TCPSocket) = set_bio!(tls, tcp)
handshake!(tls::SSLContext) = handshake(tls)

end # module


text:
@testset "Op" begin
    using NNlib: gather, gather!
    src = Float32[3, 4, 5, 6, 7]
    index = Int32[
        1 2 3 4;
        4 2 1 3;
        3 5 5 3;
    ]
    output = Float32[
        3 4 5 6;
        6 4 3 5;
        5 7 7 5;
    ]

    oa = OneHotArray(10, index)
    @test gather(src, oa) == gather(src, index)
    @test gather!(similar(src, size(index)), src, oa) == output

    w = randn(3, 10)
    @test_throws DimensionMismatch randn(3,3) * oa
    @test w * oa == gather(w, oa)

end


text:
function dispatch_ethernet(visitor, data::DenseVector{UInt8})
    ep = decode_ethernet(data)
    if ep.header.ethertype == ETHERTYPE_IPV4
        ipp = decode_ipv4(ep.payload)
        if ipp.header.protocol == IPPROTOCOL_UDP
            udpp = NetworkProtocols.decode_udp(ipp.payload)
            return visitor(ep.header, ipp.header, udpp.header, udpp.payload)
        elseif ipp.header.protocol == IPPROTOCOL_TCP
        elseif ipp.header.protocol == IPPROTOCOL_IGMP
        end
    end
    nothing
end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols


text:
# requires:

primitive type EdgeType 8 end

const EDGETYPE_TO_UINT8 = fill( 0x07, (4,4) )
      EDGETYPE_TO_UINT8[4,2] = 0x00 # 'SL' = 0x00; Tx Start
      EDGETYPE_TO_UINT8[4,3] = 0x01 # 'SR' = 0x01; Tandem Last Exon
      EDGETYPE_TO_UINT8[2,4] = 0x02 # 'LS' = 0x02; Tandem First Exon
      EDGETYPE_TO_UINT8[3,4] = 0x03 # 'RS' = 0x03; Tx End
      EDGETYPE_TO_UINT8[2,3] = 0x04 # 'LR' = 0x04; Intron 5'->3' SS
      EDGETYPE_TO_UINT8[2,2] = 0x05 # 'LL' = 0x05; Alt- 5'SS
      EDGETYPE_TO_UINT8[3,3] = 0x06 # 'RR' = 0x06; Alt- 3'SS

const INDEX_TO_EDGETYPE = transpose(reshape([[0x03,0x01,0x03,0x03];
                                             [0x06  for _ in 1:4 ];
                                             [0x05,0x04,0x05,0x02];
                                             [0x00  for _ in 1:4 ] ], (4,4)))

const INDEX_TO_EDGETYPE_NODE = transpose(reshape([[0x03  for _ in 1:4 ];
                                                  [0x06  for _ in 1:4 ];
                                                  [0x05  for _ in 1:4 ];
                                                  [0x00  for _ in 1:4 ] ], (4,4)))

function Base.convert( ::Type{EdgeType}, one::UInt8, two::UInt8 )
   @assert( 5 <= one <= 7 && 5 <= one <= 7 )
   EDGETYPE_TO_UINT8[one-3,two-3]
end

Base.convert( ::Type{EdgeType}, one::DNA, two::DNA ) =
              Base.convert( EdgeType, convert(UInt8, trailing_zeros(one)), convert(UInt8, trailing_zeros(two)) )
Base.convert( ::Type{EdgeType}, edge::UInt8 ) = reinterpret(EdgeType, edge)
Base.convert( ::Type{UInt8}, edge::EdgeType ) = reinterpret(UInt8, edge)
Base.convert( ::Type{I}, edge::EdgeType) where {I <: Integer} = Base.convert(I, Base.convert(UInt8, edge))

const EDGETYPE_SL = convert(EdgeType, 0b000)
const EDGETYPE_SR = convert(EdgeType, 0b001)
const EDGETYPE_LS = convert(EdgeType, 0b010)
const EDGETYPE_RS = convert(EdgeType, 0b011)

const EDGETYPE_LR = convert(EdgeType, 0b100)
const EDGETYPE_LL = convert(EdgeType, 0b101)
const EDGETYPE_RR = convert(EdgeType, 0b110)

istxstart( edge::EdgeType ) = edge == EDGETYPE_SL || edge == EDGETYPE_LS ? true : false
istxstop(  edge::EdgeType ) = edge == EDGETYPE_SR || edge == EDGETYPE_RS ? true : false

isexonstart( edge::EdgeType ) = edge in (EDGETYPE_LR, EDGETYPE_RR, EDGETYPE_SR, EDGETYPE_LS, EDGETYPE_SL) ? true : false
isexonstop(  edge::EdgeType ) = edge in (EDGETYPE_LR, EDGETYPE_LL, EDGETYPE_SR, EDGETYPE_LS, EDGETYPE_RS) ? true : false
ishardedge(  edge::EdgeType ) = edge in (EDGETYPE_LR, EDGETYPE_LS, EDGETYPE_SR) ? true : false

# Function takes coordinate types for node boundaries
# and returns an EdgeType
function get_edgetype( minidx::Int, secidx::Int, isnode::Bool, strand::Bool )
   if isnode
      ret = INDEX_TO_EDGETYPE_NODE[minidx,secidx]
   else
      ret = INDEX_TO_EDGETYPE[minidx,secidx]
   end
   ret = EdgeType(ret)
   if !strand
      ret = invert_edgetype( ret )
   end
   ret
end

# Invert EdgeTypes using xor
# 'LR' does not need to be inverted
function invert_edgetype( edge::EdgeType )
   if edge == EdgeType(0x04) # 'LR'
      return edge
   else
      return convert(EdgeType, convert(UInt8, edge) ⊻ 0b011)
   end
end

struct GCBiasParams
   length::Int64           # length offset
   width::Float64          # gc bin size
end

const ExpectedGC = Vector{Float16}

function ExpectedGC( seq::BioSequence{A}; gc_param = GCBiasParams(50, 0.05) ) where A <: BioSequences.Alphabet
   bins = zeros(Float16, Int(div(1.0, gc_param.width)+1))
   for i in 1:length(seq)-gc_param.length+1
      gc = Int(div(gc_content(seq[i:i+gc_param.length-1]), gc_param.width)+1)
      @inbounds bins[gc] += 1.0
   end
   # normalize columns
   binsum = sum(bins)
   binsum == 0.0 && (return bins)
   for i in 1:length(bins)
      @fastmath bins[i] /= binsum
   end
   return bins
end

# This holds a representation of the splice graph
# which is a directed multigraph
struct SpliceGraph{K}
   nodeoffset::Vector{CoordInt} # SG offset
   nodecoord::Vector{CoordInt}  # Genome offset
   nodelen::Vector{CoordInt}
   edgetype::Vector{EdgeType}
   edgeleft::Vector{SGKmer{K}}
   edgeright::Vector{SGKmer{K}}
   annopath::Vector{IntSet}
   annoname::Vector{String}
   annobias::Vector{ExpectedGC}
   seq::SGSequence
end
# All positive strand oriented sequences--->
# Node array: txStart| 1 |   2   | 3 |    4    |5| 6 |txEnd
# Edge array:        1   2       3   4         5 6   7
# Node coord:  chr   100 200     300 400      500 600 700

#= empty constructor
SpliceGraph(k::Int) = SpliceGraph( Vector{CoordInt}(), Vector{CoordInt}(),
                                   Vector{CoordInt}(), Vector{EdgeType}(),
                                   Vector{SGKmer{k}}(),Vector{SGKmer{k}}(), dna"" )
=#

# Main constructor
# Build splice graph here.
function SpliceGraph( gene::RefGene, genome::SGSequence, k::Int )
   # splice graph variables
   nodeoffset = Vector{CoordInt}()
   nodecoord  = Vector{CoordInt}()
   nodelen    = Vector{CoordInt}()
   edgetype   = Vector{EdgeType}()
   seq        = DNASequence()

   strand = gene.info.strand  # Bool

   # initialize iterators
   a,alen = 1,length(gene.acc)
   d,dlen = 1,length(gene.don)
   s,slen = 1,length(gene.txst)
   p,plen = 1,length(gene.txen)

   idx = [a, d, s, p]

   # return a tuple containing the min coordinate's idx,val
   # each of 4 categories, txstart, donor, acceptor, txend
   function getmin_ind_val( gene::RefGene, idx; def=typemax(CoordInt) )
      retarr = CoordInt[ get(gene.txen, idx[1], def),
                         get(gene.acc,  idx[2], def),
                         get(gene.don,  idx[3], def),
                         get(gene.txst, idx[4], def) ]
      indmin( retarr ), min( retarr... )
   end

   while( idx[1] <= alen || idx[2] <= dlen || idx[3] <= slen || idx[4] <= plen )
      # iterate through donors, and acceptors
      # left to right. '-' strand gets revcomp and unshift
      minidx,minval = getmin_ind_val( gene, idx )
      idx[minidx] += 1
      secidx,secval = getmin_ind_val( gene, idx )

      # last coordinate in the gene:
      # if secidx == 1 && get(gene.txst, idx[secidx], Inf) == Inf
      if secval == typemax(CoordInt)
         termedge = strand ? EdgeType(0x03) : EdgeType(0x00)
         stranded_push!(edgetype, termedge, strand)
         break
      end

      # now should we make a node?
      if issubinterval( gene.exons, Interval{CoordInt}(minval,secval) ) ||
         (minidx == 2 && minval in gene.novelacc) ||
         (secidx == 3 && secval in gene.noveldon)
         #(usebam && isretainedintron( )) # TODO (not yet implemented)
         leftadj  = (minidx == 1 || minidx == 3) && minval != secval ? 1 : 0
         rightadj = (secidx == 2 || secidx == 4) && minval != secval ? 1 : 0
         nodesize = Int(secval - minval) - leftadj - rightadj + 1
         nodeseq  = copy(genome[(Int(minval)+leftadj):(Int(secval)-rightadj)]) # collect slice
         edge     = get_edgetype( minidx, secidx, true, strand ) # determine EdgeType
         pushval  = minval + leftadj
         thridx = 0
      else # don't make a node, this is a sequence gap, make edge and inc+=2
         idx[secidx] += 1 #skip ahead again
         thridx,thrval = getmin_ind_val( gene, idx )
         rightadj = (thridx == 2 || thridx == 4) && secval != thrval ? 1 : 0
         nodesize = Int(thrval - secval) - rightadj + 1
         nodeseq  = copy(genome[Int(secval):(Int(thrval)-rightadj)])
         edge     = get_edgetype( minidx, secidx, false, strand )
         pushval  = secval
      end

      if strand
         seq *= nodeseq
      else # '-' strand
         seq = reverse_complement(nodeseq) * seq
      end
      stranded_push!(nodecoord, pushval,  strand)
      stranded_push!(nodelen,   nodesize, strand)
      stranded_push!(edgetype,  edge,     strand)

   end # end while

   # calculate node offsets-->
   curoffset = 1
   for n in nodelen
      push!(nodeoffset, curoffset)
      curoffset += n
   end

   eleft  = Vector{SGKmer{k}}(length(edgetype))
   eright = Vector{SGKmer{k}}(length(edgetype))

   paths = build_paths_edges( nodecoord, nodelen, gene )
   names = map( x->x.info.name, gene.reftx )
   expgc = map( x->ExpectedGC(path_to_seq(x, nodeoffset, nodelen, seq)), paths )

   return SpliceGraph( nodeoffset, nodecoord, nodelen, edgetype, eleft, eright, paths, names, expgc, seq )
end

# re-orient - strand by using unshift! instead of push!
function stranded_push!( collection, value, strand::Bool )
   if strand
      push!( collection, value )
   else # '-' strand
      unshift!( collection, value )
   end
end

# This function looks specifically for an intersection
# between one interval, and the interval tree, such that
# some interval in the tree fully contains the sub interval
#  itree contains: (low,             high)
#  subint must be:     (low+x, high-y)
#  where x and y are >= 0
# Returns: bool
function issubinterval( itree::T, subint::I ) where {T <: IntervalTree,
                                         I <: AbstractInterval}
   for i in intersect(itree, subint)
      if i.first <= subint.first && i.last >= subint.last
         return true
      end
   end
   false
end

# hacky extension of Base.get for tuples which isn't supported natively
# because tuples aren't really a type of collection, fair enough
function Base.get( collection::T, key::I, def ) where {T <: Tuple, I <: Integer}
   if 1 <= key <= length(collection)
      return collection[key]
   else
      return def
   end
end

# Functions for iso.jl annotated edges feature
# Take a RefTx and produce an IntSet through the nodes of a SpliceGraph
# returns: IntSet
function build_annotated_path( nodecoord::Vector{CoordInt},
                               nodelen::Vector{CoordInt},
                               tx::RefTx, strand::Bool )
   path = IntSet()
   # this may be `poor form`, but 256 is too big for default!
   resize!(path.bits, 64) # Deprecated:  = zeros(UInt32,64>>>5)
   for i in 1:length(tx.acc)
      ind = collect(searchsorted( nodecoord, tx.acc[i], rev=!strand ))[end]
      push!( path, ind )
      cur = ind + (strand ? 1 : -1)
      while 1 <= cur <= length(nodecoord) &&
            nodecoord[cur] <= tx.don[i]
         push!( path, cur )
         cur = cur + (strand ? 1 : -1 )
      end
   end
   path
end

function build_paths_edges( nodecoord::Vector{CoordInt},
                            nodelen::Vector{CoordInt},
                            gene::RefGene )
   paths = Vector{IntSet}()
   for tx in gene.reftx
      curpath = build_annotated_path( nodecoord, nodelen, tx, gene.info.strand )
      push!( paths, curpath )
   end
   paths
end

function path_to_seq( path::IntSet, nodeoffset::Vector{CoordInt},
                      nodelen::Vector{CoordInt}, seq::SGSequence )
   SGSequence(join( map( x->seq[nodeoffset[x]:(nodeoffset[x]+nodelen[x]-1)], collect(path) ) ))
end