text:
# indexing

export global_size, synchronize_threads, linear_index

# thread indexing functions
for sym in (:x, :y, :z)
    for f in (:blockidx, :blockdim, :threadidx, :griddim)
        fname = Symbol(string(f, '_', sym))
        @eval $fname(state)::Int = error("Not implemented") # COV_EXCL_LINE
        @eval export $fname
    end
end

"""
    global_size(state)

Global size == blockdim * griddim == total number of kernel execution
"""
@inline function global_size(state)
    # TODO nd version
    griddim_x(state) * blockdim_x(state)
end

"""
    linear_index(state)

linear index corresponding to each kernel launch (in OpenCL equal to get_global_id).

"""
@inline function linear_index(state)
    (blockidx_x(state) - 1) * blockdim_x(state) + threadidx_x(state)
end

"""
    linearidx(A, statesym = :state)

Macro form of `linear_index`, which calls return when out of bounds.
So it can be used like this:

    ```julia
    function kernel(state, A)
        idx = @linear_index A state
        # from here on it's save to index into A with idx
        @inbounds begin
            A[idx] = ...
        end
    end
    ```
"""
macro linearidx(A, statesym = :state)
    quote
        x1 = $(esc(A))
        i1 = linear_index($(esc(statesym)))
        i1 > length(x1) && return
        i1
    end
end

"""
    cartesianidx(A, statesym = :state)

Like [`@linearidx(A, statesym = :state)`](@ref), but returns an N-dimensional `NTuple{ndim(A), Int}` as index
"""
macro cartesianidx(A, statesym = :state)
    quote
        x = $(esc(A))
        i2 = @linearidx(x, $(esc(statesym)))
        gpu_ind2sub(x, i2)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    MaternVariogram(range=r, sill=s, nugget=n, order=ν)
    MaternVariogram(ball; sill=s, nugget=n, order=ν)

A Matérn variogram with range `r`, sill `s` and nugget `n`.
The parameter `ν` is the order of the Bessel function.
Optionally, use a custom metric `ball`.
"""
struct MaternVariogram{T,B} <: Variogram
  sill::T
  nugget::T
  order::T
  ball::B
end

MaternVariogram(ball; sill=1.0, nugget=0.0, order=1.0) =
  MaternVariogram(sill, nugget, order, ball)

MaternVariogram(; range=1.0, sill=1.0, nugget=0.0, order=1.0) =
  MaternVariogram(sill, nugget, order, MetricBall(range))

function (γ::MaternVariogram)(h::T) where {T}
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget
  ν = γ.order

  # shift lag by machine precision to
  # avoid explosion at the origin
  h′ = sqrt(2ν)*(h + eps(T)) / r
  Β  = besselk(ν, h′)
  Γ  = gamma(ν)

  (s - n) * (1 - 2^(1 - ν) / Γ * h′^ν * Β) + (h′ > 0) * n
end

isstationary(::Type{<:MaternVariogram}) = true

ext: jl


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

  E = 1.0;
  nu = 1.0/3;
  width = 48.0; height = 44.0; thickness  = 1.0;
  free_height  = 16.0;
  Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
  magn = 1.0/(free_height*thickness);# Density of applied load
  convutip = 23.97;
  n = 30;# number of elements per side
  tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

  fens,fes = T3block(width, height, n, n)

  # Reshape into a trapezoidal panel
  for i = 1:count(fens)
      fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
  end

  # Clamped edge of the membrane
  l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
  ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
  ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

  # Traction on the opposite edge
  boundaryfes =  meshboundary(fes);
  Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
  el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 2), thickness))
  flux1 = FDataDict("traction_vector"=>[0.0,+magn],
      "femm"=>el1femm
      )

  # Make the region
  MR = DeforModelRed2DStress
  material = MatDeforElastOrtho(MR,  0.0, E, nu, 0.0)
  region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
      IntegDomain(fes, TriRule(1), thickness), material))

  modeldata = FDataDict("fens"=>fens,
   "regions"=>[region1],
   "essential_bcs"=>[ess1, ess2],
   "traction_bcs"=>[flux1]
   )

  # Call the solver
  modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

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

  # Extract the solution
  nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
            inflate=tolerance);
  theutip = u.values[nl,:]
  # println("displacement =$(theutip[2]) as compared to converged $convutip")
  @test abs(theutip[2]-23.79623002934245) < 1.0e-3

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew",
     "quantity"=>:Cauchy, "component"=>:xy)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of Cauchy_xy = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.07275125002229098, 0.1309644027374448]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-vm",
     "quantity"=>:vm, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of vm = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [0.014136291979824203, 0.4181381117075297]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-pressure",
     "quantity"=>:pressure, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of pressure = $((minimum(fld.values), maximum(fld.values)))")
  # display(norm([minimum(fld.values), maximum(fld.values)] - [-0.12996180159464202, 0.2436183072544255]))
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.08664120106309468, 0.16241220483628366]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ1",
     "quantity"=>:princCauchy, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of princCauchy Max = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.10280794467415574, 0.24831137383158813]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ3",
     "quantity"=>:princCauchy, "component"=>2)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of princCauchy Min = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.4398236425818378, 0.022575693063719465]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  AE = AbaqusExporter("Cookstress_algo_stress");
  HEADING(AE, "Cook trapezoidal panel, plane stress");
  COMMENT(AE, "Converged free mid-edge displacement = 23.97");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  COMMENT(AE, "We are assuming three node triangles in plane-stress");
  COMMENT(AE, "CPS3 are pretty poor-accuracy elements, but here we don't care about it.");
@test nodesperelem(modeldata["regions"][1]["femm"].integdomain.fes) == 3
  ELEMENT(AE, "CPS3", "AllElements", connasarray(modeldata["regions"][1]["femm"].integdomain.fes))
  NSET_NSET(AE, "clamped", modeldata["essential_bcs"][1]["node_list"])
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", thickness);
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 1)
  BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 2)
  bfes = modeldata["traction_bcs"][1]["femm"].integdomain.fes
  COMMENT(AE, "Concentrated loads: we are assuming that the elements on the boundary");
  COMMENT(AE, "have two nodes each and also that they are the same length.");
  COMMENT(AE, "Then the concentrated loads below will be correctly lumped.");
  nl = connectednodes(bfes)
  F = zeros(count(modeldata["fens"]))
  for ix = 1:count(bfes)
    for jx = 1:2
      F[bfes.conn[ix][jx]] += 1.0/n/2/thickness
    end
  end
  for ixxxx = 1:length(F)
    if F[ixxxx] != 0.0
      CLOAD(AE, "ASSEM1.INSTNC1.$(ixxxx)", 2, F[ixxxx])
    end
  end
  END_STEP(AE)
  close(AE)

  nlines = 0
  open(AE.filename) do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 2886
  rm(AE.filename)
end
end
using .mmmCookmstressmorthommm
mmmCookmstressmorthommm.test()

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

  E = 1.0;
  nu = 1.0/3;
  width = 48.0; height = 44.0; thickness  = 1.0;
  free_height  = 16.0;
  Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
  magn = 1.0/(free_height*thickness);# Density of applied load
  convutip = 23.97;
  n = 30;# number of elements per side
  tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

  fens,fes = T3block(width, height, n, n)

  # Reshape into a trapezoidal panel
  for i = 1:count(fens)
      fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
  end

  # Clamped edge of the membrane
  l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
  ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
  ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

  # Traction on the opposite edge
  boundaryfes =  meshboundary(fes);
  Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
  el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 2), thickness))
  flux1 = FDataDict("traction_vector"=>[0.0,+magn],
      "femm"=>el1femm
      )

  # Make the region
  MR = DeforModelRed2DStress
  material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
  region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
      IntegDomain(fes, TriRule(1), thickness), material))

  modeldata = FDataDict("fens"=>fens,
   "regions"=>[region1],
   "essential_bcs"=>[ess1, ess2],
   "traction_bcs"=>[flux1]
   )

  # Call the solver
  modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

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

  # Extract the solution
  nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
            inflate=tolerance);
  theutip = u.values[nl,:]
  # println("displacement =$(theutip[2]) as compared to converged $convutip")
  @test abs(theutip[2]-23.79623002934245) < 1.0e-3

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew",
     "quantity"=>:Cauchy, "component"=>:xy)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of Cauchy_xy = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.07275125002229098, 0.1309644027374448]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-vm",
     "quantity"=>:vm, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of vm = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [0.014136291979824203, 0.4181381117075297]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-pressure",
     "quantity"=>:pressure, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of pressure = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.08664120106309468, 0.16241220483628366]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ1",
     "quantity"=>:princCauchy, "component"=>1)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of princCauchy Max = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.10280794467415574, 0.24831137383158813]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ3",
     "quantity"=>:princCauchy, "component"=>2)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  fld = modeldata["postprocessing"]["exported"][1]["field"]
  # println("range of princCauchy Min = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-0.4398236425818378, 0.022575693063719465]) < 1.0e-5
  # File = modeldata["postprocessing"]["exported"][1]["file"]
  # @async run(`"paraview.exe" $File`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  AE = AbaqusExporter("Cookstress_algo_stress");
  HEADING(AE, "Cook trapezoidal panel, plane stress");
  COMMENT(AE, "Converged free mid-edge displacement = 23.97");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  COMMENT(AE, "We are assuming three node triangles in plane-stress");
  COMMENT(AE, "CPS3 are pretty poor-accuracy elements, but here we don't care about it.");
@test  nodesperelem(modeldata["regions"][1]["femm"].integdomain.fes) == 3
  ELEMENT(AE, "CPS3", "AllElements", connasarray(modeldata["regions"][1]["femm"].integdomain.fes))
  NSET_NSET(AE, "clamped", modeldata["essential_bcs"][1]["node_list"])
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", thickness);
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 1)
  BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 2)
  bfes = modeldata["traction_bcs"][1]["femm"].integdomain.fes
  COMMENT(AE, "Concentrated loads: we are assuming that the elements on the boundary");
  COMMENT(AE, "have two nodes each and also that they are the same length.");
  COMMENT(AE, "Then the concentrated loads below will be correctly lumped.");
  nl = connectednodes(bfes)
  F = zeros(count(modeldata["fens"]))
  for ix = 1:count(bfes)
    for jx = 1:2
      F[bfes.conn[ix][jx]] += 1.0/n/2/thickness
    end
  end
  for ixxxx = 1:length(F)
    if F[ixxxx] != 0.0
      CLOAD(AE, "ASSEM1.INSTNC1.$(ixxxx)", 2, F[ixxxx])
    end
  end
  END_STEP(AE)
  close(AE)

  nlines = 0
  open(AE.filename) do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 2886
  rm(AE.filename)
end
end
using .mmmCookmstressisommm
mmmCookmstressisommm.test()

module mmLE10expimpmm
using FinEtools
using FinEtoolsDeforLinear
using FinEtools.MeshExportModule
using FinEtools.MeshImportModule
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")
    tolerance = Thickness/1000.; # Geometrical tolerance

    INP_file = """
    *HEADING
    NAFEMS TEST NLE10, COARSE MESH, C3D10M ELEMENTS
    *NODE
           1,     2.83277,       1.348
           2,     2.48116,     1.04837
           3,     2.12955,    0.748738
           4,     3.14146,    0.704779
           5,     2.33893,    0.399071
           6,     2.68977,    0.374369
           7,        3.25,          0.
           8,      2.8335,          0.
           9,       2.417,          0.
          10,     2.83277,       1.348,        0.15
          11,     2.48116,     1.04837,        0.15
          12,     2.12955,    0.748738,        0.15
          13,     2.33841,    0.400381,        0.15
          14,     3.14128,     0.70533,        0.15
          15,        3.25,          0.,        0.15
          16,      2.8335,          0.,        0.15
          17,       2.417,          0.,        0.15
          18,     2.83277,       1.348,         0.3
          19,     2.48116,     1.04837,         0.3
          20,     2.12955,    0.748738,         0.3
          21,     2.62488,       0.674,         0.3
          22,     2.33893,    0.399071,         0.3
          23,     3.14146,    0.704779,         0.3
          24,        3.25,          0.,         0.3
          25,      2.8335,          0.,         0.3
          26,       2.417,          0.,         0.3
          27,     2.83277,       1.348,        0.45
          28,     2.48116,     1.04837,        0.45
          29,     2.12955,    0.748738,        0.45
          30,     2.33841,    0.400381,        0.45
          31,     3.14128,     0.70533,        0.45
          32,        3.25,          0.,        0.45
          33,      2.8335,          0.,        0.45
          34,       2.417,          0.,        0.45
          35,     2.83277,       1.348,         0.6
          36,     2.48116,     1.04837,         0.6
          37,     2.12955,    0.748738,         0.6
          38,     3.14146,    0.704779,         0.6
          39,     2.33893,    0.399071,         0.6
          40,     2.68977,    0.374369,         0.6
          41,        3.25,          0.,         0.6
          42,      2.8335,          0.,         0.6
          43,       2.417,          0.,         0.6
          45,     1.95628,    0.600869
          46,     1.78302,       0.453
          47,     2.06477,    0.374369
          48,     1.93715,    0.248725
          51,      2.2085,          0.
          52,          2.,          0.
          54,     1.95628,    0.600869,        0.15
          55,     1.78302,       0.453,        0.15
          56,     1.93661,    0.249767,        0.15
          59,      2.2085,          0.,        0.15
          60,          2.,          0.,        0.15
          62,     1.95628,    0.600869,         0.3
          63,     1.78302,       0.453,         0.3
          64,     1.93715,    0.248725,         0.3
          66,     2.10001,      0.2265,         0.3
          68,      2.2085,          0.,         0.3
          69,          2.,          0.,         0.3
          71,     1.95628,    0.600869,        0.45
          72,     1.78302,       0.453,        0.45
          74,     1.93661,    0.249767,        0.45
          76,      2.2085,          0.,        0.45
          77,          2.,          0.,        0.45
          79,     1.95628,    0.600869,         0.6
          80,     1.78302,       0.453,         0.6
          81,     2.06477,    0.374369,         0.6
          83,     1.93715,    0.248725,         0.6
          85,      2.2085,          0.,         0.6
          86,          2.,          0.,         0.6
          87,       1.783,     2.29921
          88,     1.57618,     1.80182
          89,     1.36937,     1.30443
          90,     1.95627,     1.52397
          91,     2.36495,     1.88628
          92,     1.78146,     1.06985
          96,       1.783,     2.29921,        0.15
          97,     1.57618,     1.80182,        0.15
          98,     1.36937,     1.30443,        0.15
          99,     1.78071,     1.07038,        0.15
         100,     2.36449,     1.88669,        0.15
         104,       1.783,     2.29921,         0.3
         105,     1.57618,     1.80182,         0.3
         106,     1.36937,     1.30443,         0.3
         107,     2.36495,     1.88628,         0.3
         108,     1.78146,     1.06985,         0.3
         109,     2.10107,     1.32621,         0.3
         113,       1.783,     2.29921,        0.45
         114,     1.57618,     1.80182,        0.45
         115,     1.36937,     1.30443,        0.45
         116,     2.36449,     1.88669,        0.45
         117,     1.78071,     1.07038,        0.45
         121,       1.783,     2.29921,         0.6
         122,     1.57618,     1.80182,         0.6
         123,     1.36937,     1.30443,         0.6
         124,     1.95627,     1.52397,         0.6
         125,     2.36495,     1.88628,         0.6
         126,     1.78146,     1.06985,         0.6
         131,     1.26718,     1.05863
         132,       1.165,    0.812831
         134,     1.49321,    0.665266
         135,     1.64727,    0.780784
         140,     1.26718,     1.05863,        0.15
         141,       1.165,    0.812831,        0.15
         143,      1.4924,    0.665723,        0.15
         148,     1.26718,     1.05863,         0.3
         149,       1.165,    0.812831,         0.3
         150,     1.57619,    0.878714,         0.3
         152,     1.49321,    0.665266,         0.3
         157,     1.26718,     1.05863,        0.45
         158,       1.165,    0.812831,        0.45
         160,      1.4924,    0.665723,        0.45
         165,     1.26718,     1.05863,         0.6
         166,       1.165,    0.812831,         0.6
         168,     1.49321,    0.665266,         0.6
         169,     1.64727,    0.780784,         0.6
         173,          0.,       1.583
         174,          0.,      1.2915
         175,          0.,          1.
         176,      0.5825,     1.19792
         177,    0.699442,     1.51527
         178,    0.590531,    0.955415
         182,          0.,       1.583,        0.15
         183,          0.,      1.2915,        0.15
         184,          0.,          1.,        0.15
         185,    0.698861,     1.51538,        0.15
         186,    0.590076,    0.955486,        0.15
         190,          0.,       1.583,         0.3
         191,          0.,      1.2915,         0.3
         192,          0.,          1.,         0.3
         193,    0.699442,     1.51527,         0.3
         194,    0.590531,    0.955415,         0.3
         195,    0.684684,     1.15221,         0.3
         199,          0.,       1.583,        0.45
         200,          0.,      1.2915,        0.45
         201,          0.,          1.,        0.45
         202,    0.698861,     1.51538,        0.45
         203,    0.590076,    0.955486,        0.45
         207,          0.,       1.583,         0.6
         208,          0.,      1.2915,         0.6
         209,          0.,          1.,         0.6
         210,      0.5825,     1.19792,         0.6
         211,    0.699442,     1.51527,         0.6
         212,    0.590531,    0.955415,         0.6
         216,          0.,        2.75
         217,          0.,      2.1665
         219,    0.920251,     2.63745
         221,      0.8915,      1.9411
         225,          0.,        2.75,        0.15
         226,          0.,      2.1665,        0.15
         228,    0.919707,     2.63759,        0.15
         233,          0.,        2.75,         0.3
         234,          0.,      2.1665,         0.3
         236,    0.684684,     2.02721,         0.3
         237,    0.920251,     2.63745,         0.3
         242,          0.,        2.75,        0.45
         243,          0.,      2.1665,        0.45
         245,    0.919707,     2.63759,        0.45
         250,          0.,        2.75,         0.6
         251,          0.,      2.1665,         0.6
         253,    0.920251,     2.63745,         0.6
         255,      0.8915,      1.9411,         0.6
    **
    **
    *ELEMENT, TYPE=C3D10M, ELSET=EALL
           1,       1,       7,      18,       3,       4,      14,      10,
           2,       6,      11
           2,      24,      18,       7,      26,      23,      14,      15,
          25,      21,      16
           3,       9,       3,      26,       7,       5,      13,      17,
           8,       6,      16
           4,      20,      26,       3,      18,      22,      13,      12,
          19,      21,      11
           5,       3,       7,      18,      26,       6,      14,      11,
          13,      16,      21
           6,      24,      41,      18,      26,      32,      31,      23,
          25,      33,      21
           7,      35,      18,      41,      37,      27,      31,      38,
          36,      28,      40
           8,      20,      37,      26,      18,      29,      30,      22,
          19,      28,      21
           9,      43,      26,      37,      41,      34,      30,      39,
          42,      33,      40
          10,      18,      26,      41,      37,      21,      33,      31,
          28,      30,      40
          11,       9,      26,       3,      52,      17,      13,       5,
          51,      59,      47
          12,      20,       3,      26,      63,      12,      13,      22,
          62,      54,      66
          13,      46,      63,      52,       3,      55,      56,      48,
          45,      54,      47
          14,      69,      52,      63,      26,      60,      56,      64,
          68,      59,      66
          15,       3,      52,      26,      63,      47,      59,      13,
          54,      56,      66
          16,      20,      26,      37,      63,      22,      30,      29,
          62,      66,      71
          17,      43,      37,      26,      86,      39,      30,      34,
          85,      81,      76
          18,      69,      63,      86,      26,      64,      74,      77,
          68,      66,      76
          19,      80,      86,      63,      37,      83,      74,      72,
          79,      81,      71
          20,      63,      26,      37,      86,      66,      30,      71,
          74,      76,      81
          21,       1,      18,      87,       3,      10,     100,      91,
           2,      11,      90
          22,     104,      87,      18,     106,      96,     100,     107,
         105,      97,     109
          23,      89,     106,       3,      87,      98,      99,      92,
          88,      97,      90
          24,      20,       3,     106,      18,      12,      99,     108,
          19,      11,     109
          25,      87,       3,      18,     106,      90,      11,     100,
          97,      99,     109
          26,     104,      18,     121,     106,     107,     116,     113,
         105,     109,     114
          27,      35,     121,      18,      37,     125,     116,      27,
          36,     124,      28
          28,      20,     106,      37,      18,     108,     117,      29,
          19,     109,      28
          29,     123,      37,     106,     121,     126,     117,     115,
         122,     124,     114
          30,     106,      18,     121,      37,     109,     116,     114,
         117,      28,     124
          31,      89,       3,     106,     132,      92,      99,      98,
         131,     135,     140
          32,      20,     106,       3,      63,     108,      99,      12,
          62,     150,      54
          33,      46,     132,      63,       3,     134,     143,      55,
          45,     135,      54
          34,     149,      63,     132,     106,     152,     143,     141,
         148,     150,     140
          35,     132,       3,     106,      63,     135,      99,     140,
         143,      54,     150
          36,      20,      37,     106,      63,      29,     117,     108,
          62,      71,     150
          37,     123,     106,      37,     166,     115,     117,     126,
         165,     157,     169
          38,     149,     166,      63,     106,     158,     160,     152,
         148,     157,     150
          39,      80,      63,     166,      37,      72,     160,     168,
          79,      71,     169
          40,     106,      63,      37,     166,     150,      71,     117,
         157,     160,     169
          41,      89,     106,     173,     132,      98,     185,     177,
         131,     140,     176
          42,     190,     173,     106,     192,     182,     185,     193,
         191,     183,     195
          43,     175,     192,     132,     173,     184,     186,     178,
         174,     183,     176
          44,     149,     132,     192,     106,     141,     186,     194,
         148,     140,     195
          45,     173,     132,     106,     192,     176,     140,     185,
         183,     186,     195
          46,     190,     106,     207,     192,     193,     202,     199,
         191,     195,     200
          47,     123,     207,     106,     166,     211,     202,     115,
         165,     210,     157
          48,     149,     192,     166,     106,     194,     203,     158,
         148,     195,     157
          49,     209,     166,     192,     207,     212,     203,     201,
         208,     210,     200
          50,     192,     106,     207,     166,     195,     202,     200,
         203,     157,     210
          51,     216,      87,     233,     173,     219,     228,     225,
         217,     221,     226
          52,     104,     233,      87,     106,     237,     228,      96,
         105,     236,      97
          53,      89,     173,     106,      87,     177,     185,      98,
          88,     221,      97
          54,     190,     106,     173,     233,     193,     185,     182,
         234,     236,     226
          55,     173,      87,     233,     106,     221,     228,     226,
         185,      97,     236
          56,     104,     121,     233,     106,     113,     245,     237,
         105,     114,     236
          57,     250,     233,     121,     207,     242,     245,     253,
         251,     243,     255
          58,     190,     207,     106,     233,     199,     202,     193,
         234,     243,     236
          59,     123,     106,     207,     121,     115,     202,     211,
         122,     114,     255
          60,     233,     106,     121,     207,     236,     114,     245,
         243,     202,     255
    """
    write("NLE10.inp", INP_file)

    output = MeshImportModule.import_ABAQUS("NLE10.inp")
    fens, fes = output["fens"], output["fesets"][1]

    try rm("NLE10.inp"); catch end

    # 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, facing=true,
        direction=[1.0, 1.0, 0.0], dotmin= 0.001);
    topbfl = selectelem(fens, bdryfes, box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate=tolerance);

    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), TriRule(3)))
    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 = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(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[:])

    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]")
    @test norm(thecorneru - [-0.0268854 0.0 -0.0919955]) < 1.0e-5

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :invdistance, reportat = :meanonly)#
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yP = $(sigma_yP/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -2.2616980060965024) < 1.0e-3

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extrapmean)#
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yP = $(sigma_yP/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -2.382478776709117) < 1.0e-3

    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("$(fld.values[nl,1][1]/phun("MPa"))")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -5.470291697493607) < 1.0e-3

    File =  "LE10NAFEMS_MST10_sigmay.vtk"
    vtkexportmesh(File, connasarray(fes), geom.values,
                   FinEtools.MeshExportModule.VTK.T10; vectors=[("u", u.values)],
                   scalars=[("sigmay", fld.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    AE = AbaqusExporter("LE10NAFEMS_MST10");
    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, "c3d10", "AllElements", 1, connasarray(femm.integdomain.fes))
    ELEMENT(AE, "SFM3D6", "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.L1", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.L2", 2)
    BOUNDARY(AE, "ASSEM1.INSTNC1.L3", 3)
    BOUNDARY(AE, "ASSEM1.INSTNC1.L12", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.L12", 2)
    DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec([0.0, 0.0, -qmagn]))
    END_STEP(AE)
    close(AE)

    output = MeshImportModule.import_ABAQUS(AE.filename)
    fens, fes = output["fens"], output["fesets"][1]
try rm(AE.filename) catch end

    # 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, facing=true, direction=[1.0, 1.0, 0.0], dotmin= 0.001);
    topbfl = selectelem(fens, bdryfes, box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate=tolerance);

    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), TriRule(3)))
    # 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 = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(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[:])

    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]")
    @test norm(thecorneru - [-0.0268854 0.0 -0.0919955]) < 1.0e-5

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; reportat = :meanonly)#
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yP = $(sigma_yP/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -2.2616980060965024) < 1.0e-3

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extrapmean)#
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yP = $(sigma_yP/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -2.382478776709117) < 1.0e-3

    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("$(fld.values[nl,1][1]/phun("MPa"))")
    @test abs(fld.values[nl,1][1]/phun("MPa") - -5.470291697493607) < 1.0e-3

    File =  "LE10NAFEMS_MST10_sigmay.vtk"
    vtkexportmesh(File, connasarray(fes), geom.values,
                   FinEtools.MeshExportModule.VTK.T10; vectors=[("u", u.values)],
                   scalars=[("sigmay", fld.values)])
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end
end
end
using .mmLE10expimpmm
mmLE10expimpmm.test()

module mmtruncatedmfreem1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import Arpack: eigs
import LinearAlgebra: norm, cholesky, cross
function test()
    # println("""
    # Vibration modes of truncated cylindrical shell.
    # """)

    # t0 = time()

    E = 205000*phun("MPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 7850*phun("KG*M^-3");# mass density
    OmegaShift = (2*pi*100) ^ 2; # to resolve rigid body modes
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 5; nl  = 12; nc = 40;
    tolerance = h/nh/100;
    neigvs = 20;

    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2*pi,nh,nl,nc)
    # Shape into a cylinder
    R = zeros(3, 3)
    for i = 1:count(fens)
        x, y, z = fens.xyz[i,:];
        rotmat3!(R, [0, z, 0])
        Q = [cos(psi*pi/180) sin(psi*pi/180) 0;
            -sin(psi*pi/180) cos(psi*pi/180) 0;
            0 0 1]
        fens.xyz[i,:] = reshape([x+Rmed-h/2, y-l/2, 0], 1, 3)*Q*R;
    end
    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);

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

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)), material)
    femm = associategeometry!(femm, geom)
    K =stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M =mass(femm, geom, u)

    # eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
    # v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
    # of iterations niter and the number of matrix vector multiplications nmult, as
    # well as the final residual vector resid.

    if true
        d,v,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
        d[:] = d .- OmegaShift;
        fs = real(sqrt.(complex(d)))/(2*pi)
        # println("Eigenvalues: $fs [Hz]")
        @test norm(sort(fs)[1:8] - [0.0, 0.0, 0.0, 0.0, 0.000166835, 0.000182134, 517.147, 517.147]) < 2.0e-2
        # mode = 7
        # scattersysvec!(u, v[:,mode])
        # File =  "unit_cube_modes.vtk"
        # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        # @async run(`"paraview.exe" $File`)
    end

    if true
        solver = AlgoDeforLinearModule.ssit
        v0 = [i==j ? one(FFlt) : zero(FFlt) for i=1:size(K,1), j=1:2*neigvs]
        tol = 1.0e-2
        maxiter = 20
        lamb, v, nconv, niter, nmult, lamberr =
            solver(K+OmegaShift*M, M; nev=neigvs, v0=v0, tol=tol, maxiter=maxiter)
        @test nconv != neigvs
        # if nconv < neigvs
            # println("NOT converged")
        # end
        broadcast!(+, lamb, lamb, - OmegaShift);
        fs = real(sqrt.(complex(lamb)))/(2*pi)
        # println("Eigenvalues: $fs [Hz]")
        # println("$(sort(fs))")
        @test norm(sort(fs)[1:8] - [0.0, 0.0, 0.0, 0.0, 7.9048e-5, 0.0, 517.147, 517.147]) < 2.0e-2
        # println("Eigenvalue errors: $lamberr [ND]")
        # mode = 7
        # scattersysvec!(u, v[:,mode])
        # File =  "unit_cube_modes.vtk"
        # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        # @async run(`"paraview.exe" $File`)
    end

end
end
using .mmtruncatedmfreem1
mmtruncatedmfreem1.test()

module mmFV32mm1
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    # println("""
    # FV32: Cantilevered tapered membrane
    # This is a test recommended by the National Agency for Finite Element Methods and
    # Standards (U.K.): Test FV32 from NAFEMS publication TNSB, Rev. 3, “The Standard
    # NAFEMS Benchmarks,” October 1990.
    #
    # Reference solution: 44.623	130.03	162.70	246.05	379.90	391.44 for the first
    # six modes.
    # """)

    t0 = time()

    E = 200*phun("GPA");
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    L = 10*phun("M");
    W0 = 5*phun("M");
    WL = 1*phun("M");
    H = 0.05*phun("M");
    nL, nW, nH = 8, 4, 1;# How many element edges per side?
    neigvs=20                   # how many eigenvalues
    Reffs = [44.623	130.03	162.70	246.05	379.90	391.44]

    fens,fes =H20block(1.0, 2.0, 1.0, nL, nW, nH)
    for i = 1:count(fens)
        xi, eta, theta = fens.xyz[i,:];
        eta = eta - 1.0
        fens.xyz[i,:] = [xi*L eta*(1.0 - 0.8*xi)*W0/2 theta*H/2];
    end
    # File =  "mesh.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)),
      material), "femm_mass"=>FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)),
      material))

    nl1 = selectnode(fens; plane=[1.0 0.0 0.0 0.0], thickness=H/1.0e4)
    ebc1 = FDataDict("node_list"=>nl1, "component"=>1, "displacement"=>0.0)
    ebc2 = FDataDict("node_list"=>nl1, "component"=>2, "displacement"=>0.0)
    ebc3 = FDataDict("node_list"=>nl1, "component"=>3, "displacement"=>0.0)

    nl4 = selectnode(fens; plane=[0.0 0.0 1.0 0.0], thickness=H/1.0e4)
    ebc4 = FDataDict("node_list"=>nl4, "component"=>3, "displacement"=>0.0)

    # Make model data
    modeldata =  FDataDict(
      "fens"=> fens, "regions"=>  [region1], "essential_bcs"=>[ebc1 ebc2 ebc3 ebc4],
      "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtoolsDeforLinear.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
     #println("Eigenvalues: $fs [Hz]")
    # println("Percentage frequency errors: $((vec(fs[1:6]) - vec(Reffs))./vec(Reffs)*100)")
    @test norm((vec(fs[1:6]) - vec(Reffs))./vec(Reffs)*100 - [0.0162775, 0.0623384, 0.00799148, 0.151669, 0.376663, 0.0191388]) < 1.0e-6

    modeldata["postprocessing"] = FDataDict("file"=>"FV32-modes", "mode"=>1:10)
    modeldata=FinEtoolsDeforLinear.AlgoDeforLinearModule.exportmode(modeldata)
    # @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)
    try rm(modeldata["postprocessing"]["file"]*"1.vtk") catch end
end
end
using .mmFV32mm1
mmFV32mm1.test()

module mmtruncatedmfreem2
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import Arpack: eigs
import LinearAlgebra: norm, cholesky, cross
function test()
    # println("""
    # Vibration modes of truncated cylindrical shell.
    # """)

    # t0 = time()

    E = 205000*phun("MPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 7850*phun("KG*M^-3");# mass density
    OmegaShift = (2*pi*100) ^ 2; # to resolve rigid body modes
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 5; nl  = 12; nc = 40;
    tolerance = h/nh/100;
    neigvs = 20;

    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2*pi,nh,nl,nc)
    # Shape into a cylinder
    R = zeros(3, 3)
    for i = 1:count(fens)
        x, y, z = fens.xyz[i,:];
        rotmat3!(R, [0, z, 0])
        Q = [cos(psi*pi/180) sin(psi*pi/180) 0;
            -sin(psi*pi/180) cos(psi*pi/180) 0;
            0 0 1]
        fens.xyz[i,:] = reshape([x+Rmed-h/2, y-l/2, 0], 1, 3)*Q*R;
    end
    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);

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

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)), material)
    femm = associategeometry!(femm, geom)
    K =stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M =mass(femm, geom, u)

    # eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
    # v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
    # of iterations niter and the number of matrix vector multiplications nmult, as
    # well as the final residual vector resid.

    if true
        d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
        broadcast!(+, d, d, - OmegaShift);
        fs = real(sqrt.(complex(d)))/(2*pi)
        # println("Eigenvalues: $fs [Hz]")
        @test norm(sort(fs)[1:8] - [0.0, 0.0, 0.0, 0.0, 0.000166835, 0.000182134, 517.147, 517.147]) < 2.0e-2
        # mode = 7
        # scattersysvec!(u, v[:,mode])
        # File =  "unit_cube_modes.vtk"
        # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        # @async run(`"paraview.exe" $File`)
    end

    if true
        solver = AlgoDeforLinearModule.ssit
        v0 = [i==j ? one(FFlt) : zero(FFlt) for i=1:size(K,1), j=1:2*neigvs]
        tol = 1.0e-2
        maxiter = 20
        lamb, v, nconv, niter, nmult, lamberr =
            solver(K+OmegaShift*M, M; nev=neigvs, v0=v0, tol=tol, maxiter=maxiter, withrr=true)
        @test nconv == neigvs
        # if nconv < neigvs
            # println("NOT converged")
        # end
        lamb = lamb .- OmegaShift;
        fs = real(sqrt.(complex(lamb)))/(2*pi)
        # println("Eigenvalues: $fs [Hz]")
        # println("$(sort(fs))")
        @test norm(sort(fs)[1:8] - [0.0, 0.0, 0.0, 0.0, 7.9048e-5, 0.0, 517.147, 517.147]) < 2.0e-2
        # println("Eigenvalue errors: $lamberr [ND]")
        # mode = 7
        # scattersysvec!(u, v[:,mode])
        # File =  "unit_cube_modes.vtk"
        # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        # @async run(`"paraview.exe" $File`)
    end

end
end
using .mmtruncatedmfreem2
mmtruncatedmfreem2.test()


text:
using FredData
using Base.Test

# Normal usage - API key must be present in ENV
if haskey(ENV, "FRED_API_KEY")
    include("test_with_key.jl")
end


text:
# The uncertainty handling framework in this file will be added
# as part of a 1.X release. Can be ignored for now.

import UncertainData:
    resample,
    UncertainDataset,
    UncertainIndexDataset,
    UncertainValueDataset,
    UncertainIndexValueDataset
import .s_measure
import .jdd
import .transferentropy; export transferentropy
import .crossmap
import .ccm
import .predictive_asymmetry
import HypothesisTests.OneSampleTTest

##################################################
# Basic resampling for `UncertainDataset`s
##################################################
const UT = Union{UncertainValueDataset, UncertainIndexDataset, UncertainDataset}

s_measure(s::UT, t::UT, args...; kwargs...) =
    s_measure(resample(s), resample(t), args...; kwargs...)

jdd(s::UT, t::UT; kwargs...) =
    jdd(resample(s), resample(t); kwargs...)

jdd(test::OneSampleTTest, s::UT, t::UT; kwargs...) =
    jdd(test, resample(s), resample(t); kwargs...)

mutualinfo(s::UT, t::UT, method; kwargs...) =
    mutualinfo(resample(s), resample(t), method; kwargs...)

info_methods = [
    :VisitationFrequency, :TransferOperator,
    :SymbolicPermutation, :SymbolicAmplitudeAwarePermutation, :SymbolicWeightedPermutation,
    :NaiveKernel,
    :Kraskov,
    :Kraskov1,
    :Kraskov2,
    :KozachenkoLeonenko,
    :Hilbert,
    :TimeScaleMODWT
]

for method in info_methods
    @eval transferentropy(s::UT, t::UT, $(method); kwargs...) =
        transferentropy(resample(s), resample(t), method; kwargs...)

    @eval transferentropy(s::UT, t::UT, c::UT, $(method); kwargs...) =
        transferentropy(resample(s), resample(t), resample(c), method; kwargs...)
end

# transferentropy(s::UT, t::UT, method; kwargs...) =
#     transferentropy(resample(s), resample(t), method; kwargs...)

# transferentropy(s::UT, t::UT, c::UT, method; kwargs...) =
#     transferentropy(resample(s), resample(t), resample(c), method; kwargs...)

predictive_asymmetry(s::UT, t::UT, method; kwargs...) =
    predictive_asymmetry(resample(s), resample(t), method; kwargs...)

predictive_asymmetry(s::UT, t::UT, c::UT, method; kwargs...) =
    predictive_asymmetry(resample(s), resample(t), resample(c), method; kwargs...)

crossmap(s::UT, t::UT, args...; kwargs...) =
    crossmap(resample(s), resample(t), args...; kwargs...)

ccm(s::UT, t::UT, args...; kwargs...) =
    ccm(resample(s), resample(t), args...; kwargs...)

##########################################################################
# Basic resampling for `UncertainIndexValueDataset` (no constraints)
##########################################################################
const UIVD = UncertainIndexValueDataset

# TODO: warn about potential index reversals?
#
# function warn_about_sampling(s::V, t::W)
#     if s isa UIVD
#         @warn "`s` isa UncertainIndexValueDataset. Index reversals may occur. Consider constrained resampling."
#     end

#     if t isa UIVD
#         @warn "`t` isa UncertainIndexValueDataset. Index reversals may occur. Consider constrained resampling."
#     end
# end

# function warn_about_sampling(s::V, t::W, c::X)
#     warn_about_sampling(s, t)
#     if c isa UIVD
#         @warn "`c` isa UncertainIndexValueDataset. Index reversals may occur. Consider constrained resampling."
#     end
# end

s_measure(s::UIVD, t::UIVD, args...; kwargs...) =
    s_measure(resample(s.values), resample(t.values), args...; kwargs...)

jdd(s::UIVD, t::UIVD; kwargs...) =
    jdd(resample(s.values), resample(t.values); kwargs...)

jdd(test::OneSampleTTest, s::UIVD, t::UIVD; kwargs...) =
    jdd(test, resample(s), resample(t); kwargs...)

mutualinfo(s::UIVD, t::UIVD, method; kwargs...) =
    mutualinfo(resample(s.values), resample(t.values), method; kwargs...)

transferentropy(s::UIVD, t::UIVD, method; kwargs...) =
    transferentropy(resample(s.values), resample(t.values), method; kwargs...)

transferentropy(s::UIVD, t::UIVD, c::UIVD, method; kwargs...) =
    transferentropy(resample(s.values), resample(t.values), resample(c.values), method; kwargs...)

predictive_asymmetry(s::UIVD, t::UIVD, method; kwargs...) =
    predictive_asymmetry(resample(s.values), resample(t.values), method; kwargs...)

predictive_asymmetry(s::UIVD, t::UIVD, c::UIVD, method; kwargs...) =
    predictive_asymmetry(resample(s.values), resample(t.values), resample(c.values), method; kwargs...)

crossmap(s::UIVD, t::UIVD, args...; kwargs...) =
    crossmap(resample(s.values), resample(t.values), args...; kwargs...)

ccm(s::UIVD, t::UIVD, args...; kwargs...) =
    ccm(resample(s.values), resample(t.values), args...; kwargs...)


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

module Primes

using Base.Iterators: repeated

import Base: iterate, eltype, IteratorSize, IteratorEltype
using Base: BitSigned
using Base.Checked: checked_neg

export isprime, primes, primesmask, factor, ismersenneprime, isrieselprime,
       nextprime, nextprimes, prevprime, prevprimes, prime, prodfactors, radical, totient,
       SegmentedSieve

include("factorization.jl")
include("segmented_sieve/SegmentedSieve.jl")

# Primes generating functions
#     https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
#     https://en.wikipedia.org/wiki/Wheel_factorization
#     http://primesieve.org
#     Jonathan Sorenson, "An analysis of two prime number sieves", Computer Science Technical Report Vol. 1028, 1991
const wheel         = [4,  2,  4,  2,  4,  6,  2,  6]
const wheel_primes  = [7, 11, 13, 17, 19, 23, 29, 31]
const wheel_indices = [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7]

@inline function wheel_index(n)
    d, r = divrem(n - 1, 30)
    return 8d + wheel_indices[r + 2]
end
@inline function wheel_prime(n)
    d, r = (n - 1) >>> 3, (n - 1) & 7
    return 30d + wheel_primes[r + 1]
end

function _primesmask(limit::Int)
    limit < 7 && throw(ArgumentError("The condition limit ≥ 7 must be met."))
    n = wheel_index(limit)
    m = wheel_prime(n)
    sieve = ones(Bool, n)
    @inbounds for i = 1:wheel_index(isqrt(limit))
        if sieve[i]
            p = wheel_prime(i)
            q = p^2
            j = (i - 1) & 7 + 1
            while q ≤ m
                sieve[wheel_index(q)] = false
                q += wheel[j] * p
                j = j & 7 + 1
            end
        end
    end
    return sieve
end

function _primesmask(lo::Int, hi::Int)
    7 ≤ lo ≤ hi || throw(ArgumentError("The condition 7 ≤ lo ≤ hi must be met."))
    lo == 7 && return _primesmask(hi)
    wlo, whi = wheel_index(lo - 1), wheel_index(hi)
    m = wheel_prime(whi)
    sieve = ones(Bool, whi - wlo)
    hi < 49 && return sieve
    small_sieve = _primesmask(isqrt(hi))
    @inbounds for i = 1:length(small_sieve)  # don't use eachindex here
        if small_sieve[i]
            p = wheel_prime(i)
            j = wheel_index(2 * div(lo - p - 1, 2p) + 1)
            r = widemul(p, wheel_prime(j + 1))
            r > m && continue # use widemul to avoid r <= m caused by overflow
            j = j & 7 + 1
            q = Int(r)
            # q < 0 indicates overflow when incrementing q inside loop
            while 0 ≤ q ≤ m
                sieve[wheel_index(q) - wlo] = false
                q += wheel[j] * p
                j = j & 7 + 1
            end
        end
    end
    return sieve
end

"""
    primesmask([lo,] hi)

Returns a prime sieve, as a `BitArray`, of the positive integers (from `lo`, if specified)
up to `hi`. Useful when working with either primes or composite numbers.
"""
function primesmask(lo::Int, hi::Int)
    0 < lo ≤ hi || throw(ArgumentError("The condition 0 < lo ≤ hi must be met."))
    sieve = falses(hi - lo + 1)
    lo ≤ 2 ≤ hi && (sieve[3 - lo] = true)
    lo ≤ 3 ≤ hi && (sieve[4 - lo] = true)
    lo ≤ 5 ≤ hi && (sieve[6 - lo] = true)
    hi < 7 && return sieve
    wheel_sieve = _primesmask(max(7, lo), hi)
    lsi = lo - 1
    lwi = wheel_index(lsi)
    @inbounds for i = 1:length(wheel_sieve)   # don't use eachindex here
        sieve[wheel_prime(i + lwi) - lsi] = wheel_sieve[i]
    end
    return sieve
end
primesmask(lo::Integer, hi::Integer) = lo ≤ hi ≤ typemax(Int) ? primesmask(Int(lo), Int(hi)) :
    throw(ArgumentError("Both endpoints of the interval to sieve must be ≤ $(typemax(Int)), got $lo and $hi."))

primesmask(limit::Int) = primesmask(1, limit)
primesmask(n::Integer) = n ≤ typemax(Int) ? primesmask(Int(n)) :
    throw(ArgumentError("Requested number of primes must be ≤ $(typemax(Int)), got $n."))

"""
    primes([lo,] hi)

Returns a collection of the prime numbers (from `lo`, if specified) up to `hi`.
"""
function primes(lo::Int, hi::Int)
    lo ≤ hi || throw(ArgumentError("The condition lo ≤ hi must be met."))
    list = Int[]
    lo ≤ 2 ≤ hi && push!(list, 2)
    lo ≤ 3 ≤ hi && push!(list, 3)
    lo ≤ 5 ≤ hi && push!(list, 5)
    hi < 7 && return list
    lo = max(2, lo)
    sizehint!(list, 5 + floor(Int, hi / (log(hi) - 1.12) - lo / (log(lo) - 1.12 * (lo > 7))) ) # http://projecteuclid.org/euclid.rmjm/1181070157
    sieve = _primesmask(max(7, lo), hi)
    lwi = wheel_index(lo - 1)
    @inbounds for i = 1:length(sieve)   # don't use eachindex here
        sieve[i] && push!(list, wheel_prime(i + lwi))
    end
    return list
end
primes(n::Int) = primes(1, n)

const PRIMES = primes(2^16)

"""
    isprime(n::Integer) -> Bool

Returns `true` if `n` is prime, and `false` otherwise.

```julia
julia> isprime(3)
true
```
"""
function isprime(n::Integer)
    # Small precomputed primes + Miller-Rabin for primality testing:
    #     https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
    #     https://github.com/JuliaLang/julia/issues/11594
    for m in (2, 3, 5, 7, 11, 13, 17, 19, 23)
        n % m == 0 && return n == m
    end
    n < 841 && return n > 1
    s = trailing_zeros(n - 1)
    d = (n - 1) >>> s
    for a in witnesses(n)::Tuple{Vararg{Int}}
        x = powermod(a, d, n)
        x == 1 && continue
        t = s
        while x != n - 1
            (t -= 1) ≤ 0 && return false
            x = oftype(n, widemul(x, x) % n)
            x == 1 && return false
        end
    end
    return true
end

"""
    isprime(x::BigInt, [reps = 25]) -> Bool

Probabilistic primality test. Returns `true` if `x` is prime with high probability (pseudoprime);
and `false` if `x` is composite (not prime). The false positive rate is about `0.25^reps`.
`reps = 25` is considered safe for cryptographic applications (Knuth, Seminumerical Algorithms).

```julia
julia> isprime(big(3))
true
```
"""
isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp),
                                    Cint, (Any, Cint), x, reps) > 0
# TODO: Change `Any` to `Ref{BigInt}` when 0.6 support is dropped.
# The two have the same effect but `Ref{BigInt}` won't be optimized on 0.6.

# Miller-Rabin witness choices based on:
#     http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers
#     http://primes.utm.edu/prove/merged.html
#     http://miller-rabin.appspot.com
#     https://github.com/JuliaLang/julia/issues/11594
#     Forišek and Jančina, "Fast Primality Testing for Integers That Fit into a Machine Word", 2015
#         (in particular, see function FJ32_256, from which the hash and bases were taken)
const bases = UInt16[
    15591,  2018,   166,  7429,  8064, 16045, 10503,  4399,  1949,  1295,  2776,  3620,
      560,  3128,  5212,  2657,  2300,  2021,  4652,  1471,  9336,  4018,  2398, 20462,
    10277,  8028,  2213,  6219,   620,  3763,  4852,  5012,  3185,  1333,  6227,  5298,
     1074,  2391,  5113,  7061,   803,  1269,  3875,   422,   751,   580,  4729, 10239,
      746,  2951,   556,  2206,  3778,   481,  1522,  3476,   481,  2487,  3266,  5633,
      488,  3373,  6441,  3344,    17, 15105,  1490,  4154,  2036,  1882,  1813,   467,
     3307, 14042,  6371,   658,  1005,   903,   737,  1887,  7447,  1888,  2848,  1784,
     7559,  3400,   951, 13969,  4304,   177,    41, 19875,  3110, 13221,  8726,   571,
     7043,  6943,  1199,   352,  6435,   165,  1169,  3315,   978,   233,  3003,  2562,
     2994, 10587, 10030,  2377,  1902,  5354,  4447,  1555,   263, 27027,  2283,   305,
      669,  1912,   601,  6186,   429,  1930, 14873,  1784,  1661,   524,  3577,   236,
     2360,  6146,  2850, 55637,  1753,  4178,  8466,   222,  2579,  2743,  2031,  2226,
     2276,   374,  2132,   813, 23788,  1610,  4422,  5159,  1725,  3597,  3366, 14336,
      579,   165,  1375, 10018, 12616,  9816,  1371,   536,  1867, 10864,   857,  2206,
     5788,   434,  8085, 17618,   727,  3639,  1595,  4944,  2129,  2029,  8195,  8344,
     6232,  9183,  8126,  1870,  3296,  7455,  8947, 25017,   541, 19115,   368,   566,
     5674,   411,   522,  1027,  8215,  2050,  6544, 10049,   614,   774,  2333,  3007,
    35201,  4706,  1152,  1785,  1028,  1540,  3743,   493,  4474,  2521, 26845,  8354,
      864, 18915,  5465,  2447,    42,  4511,  1660,   166,  1249,  6259,  2553,   304,
      272,  7286,    73,  6554,   899,  2816,  5197, 13330,  7054,  2818,  3199,   811,
      922,   350,  7514,  4452,  3449,  2663,  4708,   418,  1621,  1171,  3471,    88,
    11345,   412,  1559,   194
]

function _witnesses(n::UInt64)
    i = xor((n >> 16), n) * 0x45d9f3b
    i = xor((i >> 16), i) * 0x45d9f3b
    i = xor((i >> 16), i) & 255 + 1
    @inbounds return (Int(bases[i]),)
end
witnesses(n::Integer) =
        n < 4294967296      ? _witnesses(UInt64(n)) :
        n < 2152302898747   ? (2, 3, 5, 7, 11) :
        n < 3474749660383   ? (2, 3, 5, 7, 11, 13) :
                              (2, 325, 9375, 28178, 450775, 9780504, 1795265022)

isprime(n::UInt128) =
    n ≤ typemax(UInt64) ? isprime(UInt64(n)) : isprime(BigInt(n))
isprime(n::Int128) = n < 2 ? false :
    n ≤ typemax(Int64)  ? isprime(Int64(n))  : isprime(BigInt(n))

# Trial division of small (< 2^16) precomputed primes +
# Pollard rho's algorithm with Richard P. Brent optimizations
#     https://en.wikipedia.org/wiki/Trial_division
#     https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
#     http://maths-people.anu.edu.au/~brent/pub/pub051.html
#
function factor!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
    # check for special cases
    if n < 0
        h[-1] = 1
        if isa(n, BitSigned) && n == typemin(T)
            h[2] = 8 * sizeof(T) - 1
            return h
        else
            return factor!(checked_neg(n), h)
        end
    elseif n == 1
        return h
    elseif n == 0 || isprime(n)
        h[n] = 1
        return h
    end

    local p::T
    for p in PRIMES
        if n % p == 0
            h[p] = get(h, p, 0) + 1
            n = div(n, p)
            while n % p == 0
                h[p] = get(h, p, 0) + 1
                n = div(n, p)
            end
            n == 1 && return h
            isprime(n) && (h[n] = 1; return h)
        end
    end
    T <: BigInt || widemul(n - 1, n - 1) ≤ typemax(n) ? pollardfactors!(n, h) : pollardfactors!(widen(n), h)
end

"""
    factor(n::Integer) -> Primes.Factorization

Compute the prime factorization of an integer `n`. The returned
object, of type `Factorization`, is an associative container whose
keys correspond to the factors, in sorted order. The value associated
with each key indicates the multiplicity (i.e. the number of times the
factor appears in the factorization).

```julia
julia> factor(100)
2^2 ⋅ 5^2
```

For convenience, a negative number `n` is factored as `-1*(-n)` (i.e. `-1` is considered
to be a factor), and `0` is factored as `0^1`:

```julia
julia> factor(-9)
-1 ⋅ 3^2

julia> factor(0)
0

julia> collect(factor(0))
1-element Array{Pair{Int64,Int64},1}:
 0=>1
```
"""
factor(n::T) where {T<:Integer} = factor!(n, Factorization{T}())

"""
    factor(ContainerType, n::Integer) -> ContainerType

Return the factorization of `n` stored in a `ContainerType`, which must be a
subtype of `AbstractDict` or `AbstractArray`, a `Set`, or an `BitSet`.

```julia
julia> factor(DataStructures.SortedDict, 100)
DataStructures.SortedDict{Int64,Int64,Base.Order.ForwardOrdering} with 2 entries:
  2 => 2
  5 => 2
```

When `ContainerType <: AbstractArray`, this returns the list
of all prime factors of `n` with multiplicities, in sorted order.

```julia
julia> factor(Vector, 100)
4-element Array{Int64,1}:
 2
 2
 5
 5

julia> prod(factor(Vector, 100)) == 100
true
```

When `ContainerType == Set`, this returns the distinct prime
factors as a set.

```julia
julia> factor(Set, 100)
Set([2,5])
```
"""
factor(::Type{D}, n::T) where {T<:Integer, D<:AbstractDict} = factor!(n, D(Dict{T,Int}()))
factor(::Type{A}, n::T) where {T<:Integer, A<:AbstractArray} = A(factor(Vector{T}, n))
factor(::Type{Vector{T}}, n::T) where {T<:Integer} =
    mapreduce(collect, vcat, [repeated(k, v) for (k, v) in factor(n)], init=Vector{T}())
factor(::Type{S}, n::T) where {T<:Integer, S<:Union{Set,BitSet}} = S(keys(factor(n)))
factor(::Type{T}, n) where {T<:Any} = throw(MethodError(factor, (T, n)))

"""
    prodfactors(factors)

Compute `n` (or the radical of `n` when `factors` is of type `Set` or
`BitSet`) where `factors` is interpreted as the result of
`factor(typeof(factors), n)`. Note that if `factors` is of type
`AbstractArray` or `Primes.Factorization`, then `prodfactors` is equivalent
to `Base.prod`.

```jldoctest
julia> prodfactors(factor(100))
100
```
"""
function prodfactors end

prodfactors(factors::AbstractDict{K}) where {K} = isempty(factors) ? one(K) : prod(p^i for (p, i) in factors)
prodfactors(factors::Union{AbstractArray, Set, BitSet}) = prod(factors)

"""
    Base.prod(factors::Primes.Factorization{T}) -> T

Compute `n` where `factors` is interpreted as the result of `factor(n)`.
"""
Base.prod(factors::Factorization) = prodfactors(factors)

"""
    radical(n::Integer)

Compute the radical of `n`, i.e. the largest square-free divisor of `n`.
This is equal to the product of the distinct prime numbers dividing `n`.

```jldoctest
julia> radical(2*2*3)
6
```
"""
radical(n) = prod(factor(Set, n))

function pollardfactors!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
    while true
        c::T = rand(1:(n - 1))
        G::T = 1
        r::K = 1
        y::T = rand(0:(n - 1))
        m::K = 1900
        ys::T = 0
        q::T = 1
        x::T = 0
        while c == n - 2
            c = rand(1:(n - 1))
        end
        while G == 1
            x = y
            for i in 1:r
                y = y^2 % n
                y = (y + c) % n
            end
            k::K = 0
            G = 1
            while k < r && G == 1
                for i in 1:(m > (r - k) ? (r - k) : m)
                    ys = y
                    y = y^2 % n
                    y = (y + c) % n
                    q = (q * (x > y ? x - y : y - x)) % n
                end
                G = gcd(q, n)
                k += m
            end
            r *= 2
        end
        G == n && (G = 1)
        while G == 1
            ys = ys^2 % n
            ys = (ys + c) % n
            G = gcd(x > ys ? x - ys : ys - x, n)
        end
        if G != n
            isprime(G) ? h[G] = get(h, G, 0) + 1 : pollardfactors!(G, h)
            G2 = div(n,G)
            isprime(G2) ? h[G2] = get(h, G2, 0) + 1 : pollardfactors!(G2, h)
            return h
        end
    end
end

"""
    ismersenneprime(M::Integer; [check::Bool = true]) -> Bool

Lucas-Lehmer deterministic test for Mersenne primes. `M` must be a
Mersenne number, i.e. of the form `M = 2^p - 1`, where `p` is a prime
number. Use the keyword argument `check` to enable/disable checking
whether `M` is a valid Mersenne number; to be used with caution.
Return `true` if the given Mersenne number is prime, and `false`
otherwise.

```jldoctest
julia> ismersenneprime(2^11 - 1)
false

julia> ismersenneprime(2^13 - 1)
true
```
"""
function ismersenneprime(M::Integer; check::Bool = true)
    if check
        d = ndigits(M, base=2)
        M >= 0 && isprime(d) && (M >> d == 0) ||
            throw(ArgumentError("The argument given is not a valid Mersenne Number (`M = 2^p - 1`)."))
    end
    M < 7 && return M == 3
    return ll_primecheck(M)
end

"""
    isrieselprime(k::Integer, Q::Integer) -> Bool

Lucas-Lehmer-Riesel deterministic test for N of the form `N = k * Q`,
with `0 < k < Q`, `Q = 2^n - 1` and `n > 0`, also known as Riesel primes.
Returns `true` if `R` is prime, and `false` otherwise or
if the combination of k and n is not supported.

```jldoctest
julia> isrieselprime(1, 2^11 - 1)  # == ismersenneprime(2^11 - 1)
false

julia> isrieselprime(3, big(2)^607 - 1)
true
```
"""
function isrieselprime(k::Integer, Q::Integer)
    n = ndigits(Q, base=2)
    0 < k < Q || throw(ArgumentError("The condition 0 < k < Q must be met."))
    if k == 1 && isodd(n)
        return n % 4 == 3 ? ll_primecheck(Q, 3) : ll_primecheck(Q)
    elseif k == 3 && (n % 4) % 3 == 0
        return ll_primecheck(Q, 5778)
    else
        # TODO: Implement a case for (k % 6) % 4 == 1 && ((k % 3) * powermod(2, n, 3)) % 3 < 2
        error("The LLR test is not currently implemented for numbers of this form.")
    end
end

# LL backend -- not for export
function ll_primecheck(X::Integer, s::Integer = 4)
    S, N = BigInt(s), BigInt(ndigits(X, base=2))
    X < 7 && throw(ArgumentError("The condition X ≥ 7 must be met."))
    for i in 1:(N - 2)
        S = (S^2 - 2) % X
    end
    return S == 0
end

"""
    totient(f::Factorization{T}) -> T

Compute the Euler totient function of the number whose prime factorization is
given by `f`. This method may be preferable to [`totient(::Integer)`](@ref)
when the factorization can be reused for other purposes.
"""
function totient(f::Factorization{T}) where T <: Integer
    if !isempty(f) && first(first(f)) == 0
        throw(ArgumentError("ϕ(0) is not defined"))
    end
    result = one(T)
    for (p, k) in f
        result *= p^(k-1) * (p - 1)
    end
    result
end

"""
    totient(n::Integer) -> Integer

Compute the Euler totient function ``ϕ(n)``, which counts the number of
positive integers less than or equal to ``n`` that are relatively prime to
``n`` (that is, the number of positive integers `m ≤ n` with `gcd(m, n) == 1`).
The totient function of `n` when `n` is negative is defined to be
`totient(abs(n))`.
"""
function totient(n::Integer)
    totient(factor(abs(n)))
end

# add: checked add (when makes sense), result of same type as first argument

add(n::BigInt, x::Int) = n + x
add(n::Integer, x::Int) = Base.checked_add(n, oftype(n, x))
sub(n::BigInt, x::Int) = n - x
sub(n::Integer, x::Int) = Base.checked_sub(n, oftype(n, x))

# add_! : "may" mutate the Integer argument (only for BigInt currently)

# modify a BigInt in-place
function add_!(n::BigInt, x::Int)
    if x < 0
        ccall((:__gmpz_sub_ui, :libgmp), Cvoid, (Ref{BigInt}, Ref{BigInt}, Culong), n, n, -x)
    else
        ccall((:__gmpz_add_ui, :libgmp), Cvoid, (Ref{BigInt}, Ref{BigInt}, Culong), n, n, x)
    end
    n
end

# checked addition, without mutation
add_!(n::Integer, x::Int) = add(n, x)
sub_!(n::BigInt, x::Int) = add_!(n, -x)
sub_!(n::Integer, x::Int) = sub(n, x)

"""
    nextprime(n::Integer, i::Integer=1; interval::Integer=1)

The `i`-th smallest prime not less than `n` (in particular,
`nextprime(p) == p` if `p` is prime). If `i < 0`, this is equivalent to
prevprime(n, -i). Note that for `n::BigInt`, the returned number is
only a pseudo-prime (the function [`isprime`](@ref) is used
internally). See also [`prevprime`](@ref).

If `interval` is provided, primes are sought in increments of `interval`.
This can be useful to ensure the presence of certain divisors in `p-1`.
The range of possible values for `interval` is currently `1:typemax(Int)`.

```jldoctest
julia> nextprime(4)
5

julia> nextprime(5)
5

julia> nextprime(4, 2)
7

julia> nextprime(5, 2)
7

julia> nextprime(2^10+1; interval=1024)
12289

julia> gcd(12289 - 1, 1024) # 1024 | p - 1
1024
```
"""
function nextprime(n::Integer, i::Integer=1; interval::Integer=1)
    i == 0 && throw(DomainError(i))
    i < 0 && return prevprime(n, -i; interval=interval)
    interval < 1 && throw(DomainError(interval, "interval must be >= 1"))
    # TODO: lift the following condition
    interval > typemax(Int) && throw(DomainError(interval, "interval must be <= $(typemax(Int))"))
    interval = oftype(n, interval)
    if n < 2
        n = interval == 1 ?
            oftype(n, 2) :
            # smallest value >= 2 whose difference from n is a multiple of interval
            oftype(n, n + interval * (1 + (n-1)÷(-interval)))
    end
    if n == 2
        if i <= 1
            return n
        else
            n += interval
            i -= 1
        end
    else
        n += iseven(n) ? interval : zero(n)
    end
    # n can now be safely mutated
    # @assert n >= 3
    if iseven(n)
        @assert iseven(interval)
        throw(DomainError((n, interval),
            "`n` and `interval` should not be both even (there is then no correct answer)."))
    end
    # @assert isodd(n)
    interval = Int(interval)
    isodd(interval) && (interval = Base.checked_mul(interval, 2))

    while true
        while !isprime(n)
            n = add_!(n, interval)
        end
        i -= 1
        i <= 0 && break
        n = add_!(n, interval)
    end
    n
end

"""
    prevprime(n::Integer, i::Integer=1; interval::Integer=1)

The `i`-th largest prime not greater than `n` (in particular
`prevprime(p) == p` if `p` is prime). If `i < 0`, this is equivalent to
`nextprime(n, -i)`. Note that for `n::BigInt`, the returned number is
only a pseudo-prime (the function [`isprime`](@ref) is used internally). See
also [`nextprime`](@ref).

If `interval` is provided, primes are sought in increments of `interval`.
This can be useful to ensure the presence of certain divisors in `p-1`.
The range of possible values for `interval` is currently `1:typemax(Int)`.

```jldoctest
julia> prevprime(4)
3

julia> prevprime(5)
5

julia> prevprime(5, 2)
3
```
"""
function prevprime(n::Integer, i::Integer=1; interval::Integer=1)
    i <= 0 && return nextprime(n, -i; interval=interval)
    interval < 1 && throw(DomainError(interval, "interval must be >= 1"))
    interval > typemax(Int) && throw(DomainError(interval, "interval must be <= $(typemax(Int))"))
    interval = Int(interval)

    n += zero(n) # deep copy of n, which is mutated below

    # A bit ugly, but this lets us skip half of the isprime tests when isodd(interval)
    @inline function decrement(n)
        n = sub_!(n, interval)
        iseven(n) && n != 2 ? # n obviously not prime
            sub_!(n, interval) :
            n
    end

    while true
        while !isprime(n)
            n < 2 && throw(ArgumentError("There is no prime less than or equal to $n"))
            n = decrement(n)
        end
        i -= 1
        i <= 0 && break
        n = decrement(n)
    end
    n
end

"""
    prime(::Type{<:Integer}=Int, i::Integer)

The `i`-th prime number.

```jldoctest
julia> prime(1)
2

julia> prime(3)
5

```
"""
prime(::Type{T}, i::Integer) where {T<:Integer} = i < 0 ? throw(DomainError(i)) : nextprime(T(2), i)
prime(i::Integer) = prime(Int, i)

struct NextPrimes{T<:Integer}
    start::T
end

function iterate(np::NextPrimes, state=np.start)
    p = nextprime(state)
    (p, add(p, 1))
end

IteratorSize(::Type{<:NextPrimes}) = Base.IsInfinite()
IteratorEltype(::Type{<:NextPrimes}) = Base.HasEltype()

eltype(::Type{NextPrimes{T}}) where {T} = T

"""
    nextprimes(start::Integer)

Return an iterator over all primes greater than or equal to `start`,
in ascending order.
"""
nextprimes(start::Integer) = NextPrimes(start)

"""
    nextprimes(T::Type=Int)

Return an iterator over all primes, with type `T`.
Equivalent to `nextprimes(T(1))`.
"""
nextprimes(::Type{T}=Int) where {T<:Integer} = nextprimes(one(T))

"""
    nextprimes(start::Integer, n::Integer)

Return an array of the first `n` primes greater than or equal to `start`.

# Example

```
julia> nextprimes(10, 3)
3-element Array{Int64,1}:
 11
 13
 17
```
"""
nextprimes(start::T, n::Integer) where {T<:Integer} =
    collect(T, Iterators.take(nextprimes(start), n))

struct PrevPrimes{T<:Integer}
    start::T
end

function iterate(np::PrevPrimes, state=np.start)
    if isone(state)
        nothing
    else
        p = prevprime(state)
        (p, p-one(p))
    end
end

IteratorSize(::Type{<:PrevPrimes}) = Base.SizeUnknown()
IteratorEltype(::Type{<:PrevPrimes}) = Base.HasEltype()

eltype(::Type{PrevPrimes{T}}) where {T} = T

"""
    prevprimes(start::Integer)

Return an iterator over all primes less than or equal to `start`,
in descending order.

# Example

```
julia> collect(prevprimes(10))
4-element Array{Int64,1}:
 7
 5
 3
 2
```
"""
prevprimes(start::Integer) = PrevPrimes(max(one(start), start))

"""
    prevprimes(start::Integer, n::Integer)

Return an array of the first `n` primes less than or equal to `start`,
in descending order. When there are less than `n` primes less than or
equal to `start`, those primes are returned without an error.

# Example

```
julia> prevprimes(10, 3)
3-element Array{Int64,1}:
 7
 5
 3

julia> prevprimes(10, 10)
3-element Array{Int64,1}:
 7
 5
 3
 2
```
"""
prevprimes(start::T, n::Integer) where {T<:Integer} =
    collect(T, Iterators.take(prevprimes(start), n))

end # module


text:
using StatisticalRethinking
using DiffEqBayes, DynamicHMC, DataFrames
using BenchmarkTools
using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions

using StatsPlots
gr(size=(600,900))

ProjDir = @__DIR__
cd(ProjDir)
isdir(ProjDir*"/tmp") && rm(ProjDir*"/tmp", recursive=true)

f = @ode_def LotkaVolterraTest begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [30.0, 4.0]
tspan = (0.0, 21.0)
p = [0.55, 0.028, 0.84, 0.026]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(1,stop=20,length=20))

df = CSV.read("$(ProjDir)/../lynx_hare.csv", delim=",")

lv = "
  functions {
    real[] dz_dt(real t,       // time
                 real[] z,     // system state {prey, predator}
                 real[] theta, // parameters
                 real[] x_r,   // unused data
                 int[] x_i) {
      real u = z[1];
      real v = z[2];

      real alpha = theta[1];
      real beta = theta[2];
      real gamma = theta[3];
      real delta = theta[4];

      real du_dt = (alpha - beta * v) * u;
      real dv_dt = (-gamma + delta * u) * v;

      return { du_dt, dv_dt };
    }
  }
  data {
    int<lower = 0> N;          // number of measurement times
    real ts[N];                // measurement times > 0
    real y_init[2];            // initial measured populations
    real<lower = 0> y[N, 2];   // measured populations
  }
  parameters {
    real<lower = 0> theta[4];   // { alpha, beta, gamma, delta }
    real<lower = 0> z_init[2];  // initial population
    real<lower = 0> sigma[2];   // measurement errors
  }
  transformed parameters {
    real z[N, 2]
      = integrate_ode_rk45(dz_dt, z_init, 0, ts, theta,
                           rep_array(0.0, 0), rep_array(0, 0),
                           1e-5, 1e-3, 5e2);
  }
  model {
    theta[{1, 3}] ~ normal(1, 0.5);
    theta[{2, 4}] ~ normal(0.05, 0.05);
    sigma ~ lognormal(-1, 1);
    z_init ~ lognormal(10, 1);
    for (k in 1:2) {
      y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
      y[ , k] ~ lognormal(log(z[, k]), sigma[k]);
    }
  }
  generated quantities {
    real y_init_rep[2];
    real y_rep[N, 2];
    for (k in 1:2) {
      y_init_rep[k] = lognormal_rng(log(z_init[k]), sigma[k]);
      for (n in 1:N)
        y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
    }
  }
";

tmpdir = ProjDir*"/tmp"
sm = SampleModel("lynx_hare", lv; tmpdir=tmpdir)

N = size(df, 1) -1
lv_data = Dict(
  "N" => N,
  "ts" => collect(1:N),
  "y_init" => [30, 4],
  "y" => Array(df[2:end, [:Hare, :Lynx]])
)

@time rc = stan_sample(sm, data=lv_data)

if success(rc)

  dfa = read_samples(sm, output_format=:dataframe)
  p = Particles(dfa)
  p[Symbol("theta.1")] |> display
  p[Symbol("theta.2")] |> display
  p[Symbol("theta.3")] |> display
  p[Symbol("theta.4")] |> display
  p[Symbol("sigma.1")] |> display
  p[Symbol("sigma.2")] |> display
  p[Symbol("z_init.1")] |> display
  p[Symbol("z_init.2")] |> display

  println()
  sdf = StanSample.read_summary(sm)
  sdf[8:15, [1, 2, 4, 5, 7, 8, 9]] |> display

  original_stan_result = "
                mean se_mean    sd    10%    50%    90% n_eff Rhat
    theta[1]   0.549   0.002 0.065  0.469  0.545  0.636  1163    1
    theta[2]   0.028   0.000 0.004  0.023  0.028  0.034  1281    1
    theta[3]   0.797   0.003 0.091  0.684  0.791  0.918  1125    1
    theta[4]   0.024   0.000 0.004  0.020  0.024  0.029  1170    1
    sigma[1]   0.248   0.001 0.045  0.198  0.241  0.306  2625    1
    sigma[2]   0.252   0.001 0.044  0.201  0.246  0.310  2808    1
    z_init[1] 33.960   0.056 2.909 30.363 33.871 37.649  2684    1
    z_init[2]  5.949   0.011 0.533  5.294  5.926  6.644  2235    1
  ";

  p1 = plot()
  scatter!(vcat(0, t), df[:, :Hare], lab="Obs hare")
  scatter!(vcat(0, t), vcat(30, [mean(dfa[:, Symbol("y_rep.$i.1")]) for i in 1:20]),
    lab="Pred hare")
  hares = transpose(convert(Array, VectorOfArray(vcat(
    [hpdi(dfa[:, Symbol("z_init.1")])],
    [hpdi(dfa[:, Symbol("y_rep.$i.1")]) for i in 1:20]
  ))))
  mu = vcat(
    mean(dfa[:, Symbol("z_init.1")]),
    [mean(dfa[:, Symbol("y_rep.$i.1")]) for i in 1:20]
  )
  plot!(vcat(0, t), mu; ribbon=hares, alpha=0.4, lab="89% hpd", color=:lightgrey)

  plot!(sol)

  p2 = plot()
  scatter!(vcat(0, t), vcat(4, [mean(dfa[:, Symbol("y_rep.$i.2")]) for i in 1:20]),
    lab="Pred lynx")
  scatter!(vcat(0, t), df[:, :Lynx], lab="Obs Lynx")
  lynxs = transpose(convert(Array, VectorOfArray(vcat(
    [hpdi(dfa[:, Symbol("z_init.2")])],
    [hpdi(dfa[:, Symbol("y_rep.$i.2")]) for i in 1:20]
  ))))
  mu = vcat(
    mean(dfa[:, Symbol("z_init.2")]),
    [mean(dfa[:, Symbol("y_rep.$i.2")]) for i in 1:20]
  )
  plot!(vcat(0, t), mu; ribbon=lynxs, alpha=0.4, lab="89% hpd", color=:lightgrey)

  plot!(sol)

  plot(p1, p2, layout=(2,1))
  savefig("$(ProjDir)/fig_01.png")

end


text:
addparent(A::Type{<:AbstractArray}, P′::DataType) = A
addparent(::Type{CuArray{T,N,P}}, P′::DataType) where {T,N,P} = CuArray{T,N,P′}

const StatefulOptimiser = Union{Flux.Momentum, Flux.Nesterov, Flux.RMSProp, Flux.ADAM, Flux.RADAM, Flux.AdaMax, Flux.ADAGrad, Flux.ADADelta, Flux.AMSGrad, Flux.NADAM}

function Base.show(io::IO, optimiser::StatefulOptimiser)
    print(io, typeof(optimiser), getproperty.(Ref(optimiser), propertynames(optimiser)[1:end-1]))
end

@adjoint function reduce(::typeof(hcat), As::AbstractVector{<:AbstractVecOrMat})
   cumsizes = cumsum(size.(As, 2))
   return reduce(hcat, As), Δ -> (nothing, map((sz, A) -> Zygote.pull_block_horz(sz, Δ, A), cumsizes, As))
end

@adjoint function reduce(::typeof(vcat), As::AbstractVector{<:AbstractVecOrMat})
   cumsizes = cumsum(size.(As, 1))
   return reduce(vcat, As), Δ -> (nothing, map((sz, A) -> Zygote.pull_block_vert(sz, Δ, A), cumsizes, As))
end

"""
    tensor2vecofmats(Xs::DenseArray{<:Real,3}) -> Vector{<:DenseVecOrMat}

Given a 3D tensor `Xs` of dimensions D×T×B, constructs a vector of length T of D×B slices of `Xs` along the second dimension.
"""
tensor2vecofmats(Xs::DenseArray{<:Real,3}) = [Xs[:,t,:] for t ∈ axes(Xs, 2)]

"""
    vecofmats2tensor(xs::DenseVector{<:DenseVecOrMat}) -> DenseArray{<:Real,3}

Constructs a 3D tensor of dimensions D×T×B by concatenating along the second dimension the elements of the input vector `xs` of length T, whose each element is a D×B matrix.
"""
function vecofmats2tensor(xs::DenseVector{<:DenseVecOrMat})
   x₁ = first(xs)
   D, B = size(x₁)
   Xs = similar(x₁, D, length(xs), B)
   @inbounds for t ∈ eachindex(xs)
      Xs[:,t,:] = xs[t]
   end
   return Xs
end

function pad(xs::DenseVector{<:DenseVector}, multiplicity::Integer)
   T = length(xs)
   newT = ceil(Int, T / multiplicity)multiplicity
   z = zero(first(xs))
   return [xs; (_ -> z).(1:(newT - T))]
end

"""
    batch_inputs!(Xs::AbstractVector{<:AbstractVector{<:DenseVector}}, multiplicity::Integer, maxT::Integer = maximum(length, Xs))::DenseArray{<:Real,3}

Given a B-length collection `Xs`, whose each element is a sequence vector of D-dimensional input features, pads sequences to ensure every sequence has the same length `T` that is a multiple of `multiplicity`. Then it arranges these padded sequences into a single tensor of size D×T×B.
"""
function batch_inputs!(Xs::AbstractVector{<:AbstractVector{<:DenseVector}}, multiplicity::Integer, maxT::Integer = maximum(length, Xs))::DenseArray{<:Real,3}
   # find the smallest multiple of `multiplicity` that is no less than `maxT`
   newT = ceil(Int, maxT / multiplicity)multiplicity
   # construct padding element vector
   z = (zero ∘ first ∘ first)(Xs)
   # resize each sequence `xs` to the size `newT` by paddding it with vector z of zeros
   for xs ∈ Xs
      T = length(xs)
      resize!(xs, newT)
      xs[(T+1):end] .= Ref(z)
   end
   # first concatenate D-dimensional input features of each padded sequence of length newT into the single vector of length D*newT, then concatenate the resulting vectors along the 2nd dimension to get the D*newT×B matrix and then finally reshape the resulting matrix into D×newT×B tensor
      X = reshape(reduce(hcat, reduce.(vcat, Xs)), length(z), newT, :)
      # check: X == cat(reduce.(hcat, Xs)...; dims=3)
   return X
end

"""
    batch_targets(ys::AbstractVector{<:DenseVector{<:Integer}}, dim_out::Integer, maxT::Integer = maximum(length, ys))

Given a batch vector of target sequences `ys` returns a vector of corresponding linear indexes into the prediction `Ŷs`, which is assumed to be a tensor od dimensions D×B×T. Here D denotes the dimensionality of the output, B is the batch size and T is the maximum time length in the batch.
"""
function batch_targets(ys::AbstractVector{<:DenseVector{<:Integer}}, dim_out::Integer, maxT::Integer = maximum(length, ys))
   batch_size = length(ys)
   cartesian_indices = Vector{Vector{CartesianIndex{3}}}(undef, maxT)
   cartesian_indices_t = Vector{CartesianIndex{3}}(undef, batch_size)
   @views for time ∈ 1:maxT
      n = 0
      for (batch, yᵇ) ∈ enumerate(ys)
         if time <= length(yᵇ)
            n += 1
            cartesian_indices_t[n] = CartesianIndex(yᵇ[time], batch, time)
         end
      end
      cartesian_indices[time] = cartesian_indices_t[1:n]
   end
   cartesian_indices′ = reduce(vcat, cartesian_indices)
   linear_indices = LinearIndices((dim_out, batch_size, maxT))[cartesian_indices′]
   @assert issorted(linear_indices)
   return cartesian_indices′
end

"""
    batch_dataset(Xs::DenseVector{<:DenseVector{<:DenseVector}}, ys::DenseVector{<:DenseVector}, batch_size::Integer, multiplicity::Integer)

Arranges dataset into batches such that the number of batches approximately equals the ratio of dataset size to `batch_size`.
Batches are formed by first sorting sequences in the dataset according to their length (which minimizes the total number of elements to pad in inputs) and then partitioning the result into batches such that each batch approximately the same total number of sequence elements (this ensures that each batch takes up the same amount of memory, so as to avoid memory overflow when loading data into the GPU).
"""
function batch_dataset(Xs::DenseVector{<:DenseVector{<:DenseVector}},
                       ys::DenseVector{<:DenseVector},
                       dim_out::Integer,
                       batch_size::Integer,
                       multiplicity::Integer)

   sortidxs = sortperm(Xs; by=length)
   Xs, ys = Xs[sortidxs], ys[sortidxs]

   cumseqlengths = cumsum(length.(ys))
   nbatches = ceil(Int, length(Xs) / batch_size)
   # subtract 0.5 from the last element of the range
   # to ensure that i index inside the loop won't go out of bounds due to floating point rounding errors
   cum_n_elts_rng = range(0, last(cumseqlengths)-0.5; length = nbatches+1)[2:end]
   lastidxs = similar(sortidxs, nbatches)
   i = 1
   for (n, cum_n_elts_for_batch) ∈ enumerate(cum_n_elts_rng)
      while cumseqlengths[i] < cum_n_elts_for_batch
         i += 1
      end
      lastidxs[n] = i
   end
   firstidxs = [1; @view(lastidxs[1:(end-1)]) .+ 1]

   maxTs = length.(@view Xs[lastidxs])
   X_batches = [ batch_inputs!(Xs[firstidx:lastidx], multiplicity, maxT) for (firstidx, lastidx, maxT) ∈ zip(firstidxs, lastidxs, maxTs) ]
   indices_batches = [ batch_targets(ys[firstidx:lastidx], dim_out, maxT) for (firstidx, lastidx, maxT) ∈ zip(firstidxs, lastidxs, maxTs) ]

   batches = [(X |> gpu, indices, maxT) for (X, indices, maxT) ∈ zip(X_batches, indices_batches, maxTs)]
   return batches
end

"""
    truthprob(l::Real, total_length::Integer)::Real
    truthprob(l::Real, indices::DenseVector)::Real
    truthprob(l::Real, dataset::Vector{<:Tuple{DenseArray{<:Real,3}, DenseVector{<:Integer}, Integer}})::Real

Given a loss `l` for either a batch of length `total_length` or a batch with linear indices `indices` of correct labels or a collection of batches, `dataset`, returns mean probability of the correct prediction
"""
truthprob(l::Real, total_length::Integer)::Real = exp(-l / total_length)
truthprob(l::Real, indices::DenseVector)::Real = exp(-l / length(indices))
truthprob(l::Real, dataset::Vector{<:Tuple{DenseArray{<:Real,3}, DenseVector{<:Integer}, Integer}})::Real =
   exp(-l / sum(((_, indices, _),) -> length(indices), dataset))

function printlog(io::IO, message...)
   println(io, message...)
   flush(io)
   println(message...)
end

function fill_param_dict!(dict::AbstractDict, rm::Recur, name::AbstractString="")
   verbose = string(rm)
   short = replace(verbose, r"^Recur\(|Cell|\)$" =>"")
   name = replace(name, verbose => short)
   fill_param_dict!(dict, rm.cell, name)
end
function fill_param_dict!(dict::AbstractDict, x::AbstractArray{<:Number}, name::AbstractString)
   dict[name] = vec(x)
   return nothing
end
function fill_param_dict!(dict::AbstractDict, model::Chain, name::AbstractString)
   for (i, layer) ∈ enumerate(model.layers)
      fill_param_dict!(dict, layer, name * "/layer_$i:$layer")
   end
end
function fill_param_dict!(dict::AbstractDict, x, name::AbstractString)
   for (field, child) in pairs(Flux.trainable(x))
      fill_param_dict!(dict, child, name * "/$field")
   end
end
function param_dict(model, prefix::AbstractString="")
   dict = Dict{String,Any}()
   fill_param_dict!(dict, model, prefix)
   return sort(dict)
end


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

using Test, LinearAlgebra, SparseArrays
let ambig = detect_ambiguities(LinearAlgebra; recursive=true)
    @test isempty(ambig)
    ambig = Set{Any}(((m1.sig, m2.sig) for (m1, m2) in ambig))
    expect = []
    good = true
    while !isempty(ambig)
        sigs = pop!(ambig)
        i = findfirst(==(sigs), expect)
        if i === nothing
            println(stderr, "push!(expect, (", sigs[1], ", ", sigs[2], "))")
            good = false
            continue
        end
        deleteat!(expect, i)
    end
    @test isempty(expect)
    @test good
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using LinearAlgebra
using SparseArrays

function expmv!(
        f::AbstractVecOrMat, t::StepRangeLen, A, b::AbstractVecOrMat,
        M = nothing, norm = LinearAlgebra.norm, opnorm = LinearAlgebra.opnorm;
        precision = :double, shift = DEFAULT_SHIFT, check_positive = false
    )

    n = size(A, 2)
    T = eltype(b)
    Tr = real(T)
    q = t.len - 1

    t0 = Tr(t.ref)
    tmax = Tr(t.ref + (t.len - 1.) * t.step)

    if shift == true && !hasmethod(tr, Tuple{typeof(A)})
        shift = false
        @warn "Shift set to false as $(typeof(A)) doesn't support tr"
    end

    force_estm = !hasmethod(opnorm, Tuple{typeof(A), typeof(1)})
    Anorm_dt = (tmax - t0) * normAm(A, 1; check_positive = check_positive);

    if (precision == :single || precision == :half && Anorm_dt > 85.496) || ( precision == :double && Anorm_dt > 63.152)
       force_estm = true;
    end

    if M == nothing
        (M, alpha, unA) = select_taylor_degree(A, size(b,2); precision=precision, shift=shift, force_estm=force_estm, check_positive=check_positive)
    end

    tol =
      if precision == :double
          Tr(2.0)^(-53)
      elseif precision == :single
          Tr(2.0)^(-24)
      elseif precision == :half
          Tr(2.0)^(-10)
      end

    X = zeros(T, n, q+1);
    (m_max, p) = size(M);
    U = diagm(0 => 1:m_max);

    _, s = degree_selector(tmax - t0, M, U, p)
    h = (tmax - t0)/q;

    expmv!(view(X,:,1), t0, A, b, M, norm, opnorm;
        precision=precision, shift=shift, check_positive=check_positive);

    mu = zero(Tr)
    if shift
        mu = tr(A)/n
        A = diagshift!(A, mu) #A = A - mu * I
    end

    d = max(1, Integer(floor(q/s)))
    j = Integer(floor(q/d))
    r = q - d * j
    m_opt, _ = degree_selector(d, M, U, p)
    dr = d

    z = X[:,1]
    K = zeros(T, n, m_opt+1)
    temp = similar(b)

    for i = 1:j+1
        if i > j
            dr = r
        end
        @views K[:,1] .= z
        m = 0
        for k = 1:dr
            f = copyto!(f, z)
            c1 = norm(z, Inf)

            p = 1
            for p = 1:m_opt # outer?
                if p > m
                    temp = mul!(temp, A, view(K,:,p)) # A*K[:,p]
                    @views K[:,p+1] .= (h/p) .* temp
                end

                @views temp .= (Tr(k)^p) .* K[:,p+1]
                f .+= temp
                c2 = norm(temp, Inf)
                if c1 + c2 <= tol * norm(f, Inf)
                    break
                end
                c1 = c2
            end

            m = max(m,p)
            if shift
                @views X[:, k + (i-1)*d + 1] .= exp(k*h*mu) .* f
            else
                @views X[:, k + (i-1)*d + 1] .= f
            end
        end

        if i <= j
            @views z .= X[:, i*d + 1]
        end
    end

    return X
end

function expmv(
        t::StepRangeLen, A, b::AbstractVecOrMat,
        M = nothing, norm = LinearAlgebra.norm, opnorm = LinearAlgebra.opnorm;
        precision = :double, shift = DEFAULT_SHIFT, check_positive = false
    )

    return expmv!(
        similar(b), t, A, b, M, norm, opnorm;
        precision = precision, shift = shift, check_positive = check_positive
    )
end