text:
module m111ocylpull14nnn # From miscellaneous
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    # Cylinder  compressed by enforced displacement, axially symmetric model

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

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] = fens.xyz[:, 1] .+ rin
    bdryfes = meshboundary(fes);

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

    # the symmetry plane
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    # The other end
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, ua)

    applyebc!(u)
    numberdofs!(u)
    # println("Number of degrees of freedom = $(u.nfreedofs)")
    @test u.nfreedofs == 240

    # Property and material
    material = MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)
    # display(material)
    # println("$(material.D)")
    # @show MR

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

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

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - 0.0) < 1.0e-5
    @test abs(maximum(fld.values) - 0.0) < 1.0e-5
    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 2)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - 0.0) < 1.0e-5
    @test abs(maximum(fld.values) - 0.0) < 1.0e-5
    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 3)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - -0.050) < 1.0e-5
    @test abs(maximum(fld.values) - -0.04999999999999919) < 1.0e-5

    # File =  "mocylpull14.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    #               vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
end
end
using .m111ocylpull14nnn
m111ocylpull14nnn.test()

module mophysun13 # From miscellaneous
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    E1=1.0;
    nu23=0.19;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    tolerance=rin/1000.

    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] = fens.xyz[:, 1] .+ rin
    bdryfes = meshboundary(fes);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2))
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)
    @test u.nfreedofs == 240

    material=MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)
    # println("success? ")
    # @code_llvm FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material, true)
    # femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material, true)
    # println("failure? ")
    # @code_llvm FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    true
end
end
using .mophysun13
mophysun13.test()

module m111ocylpull14n1 # From miscellaneous
using FinEtools
using FinEtoolsDeforLinear
using Test

function test()
    E1=1.0;
    nu23=0.19;
    rin=1.;
    rex =1.2;
    Length = 1*rex

    MR = DeforModelRed2DAxisymm
    fens,fes = Q4block(rex-rin,Length,5,20);
    material = MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)

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

    true
end
end
using .m111ocylpull14n1
m111ocylpull14n1.test()

module mocylpull14a
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    MR = DeforModelRed2DAxisymm
    material = MatDeforElastIso(MR, 0.0, 1.0, 0.0, 0.0)
    @test MR == material.mr
    femm = FEMMDeforLinear(MR, IntegDomain(FESetP1(reshape([1],1,1)), GaussRule(2, 2), true), material)

end
end
using .mocylpull14a
mocylpull14a.test()

module mstressconversionm
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm
function test()
  symmtens(N) = begin t=rand(N, N); t = (t+t')/2.0; end
  t = symmtens(2)
  v = zeros(3)
  strainttov!(DeforModelRed2DStrain, v, t)
  to = zeros(2, 2)
  strainvtot!(DeforModelRed2DStrain, to, v)
  @test norm(t-to) < eps(1.0)

  t = symmtens(3)
  v = zeros(6)
  strainttov!(DeforModelRed3D, v, t)
  to = zeros(3, 3)
  strainvtot!(DeforModelRed3D, to, v)
  @test norm(t-to) < eps(1.0)

  t = symmtens(2)
  v = zeros(3)
  stressttov!(DeforModelRed2DStress, v, t)
  to = zeros(2, 2)
  stressvtot!(DeforModelRed2DStress, to, v)
  @test norm(t-to) < eps(1.0)

  v = vec([1. 2. 3.])
  t = zeros(3, 3)
  stressvtot!(DeforModelRed2DStrain, t, v)
  to = [1. 3. 0; 3. 2. 0; 0 0 0]
  @test norm(t-to) < eps(1.0)

  v = vec([1. 2 3 4])
  t = zeros(3, 3)
  stressvtot!(DeforModelRed2DStrain, t, v)
  to = [1. 3 0; 3 2 0; 0 0 4]
  @test norm(t-to) < eps(1.0)

  v = rand(6)
  t = zeros(3, 3)
  stressvtot!(DeforModelRed3D, t, v)
  vo = zeros(6)
  stressttov!(DeforModelRed3D, vo, t)
  @test norm(v-vo) < eps(1.0)

  # v = rand(9)
  # t = zeros(3, 3)
  # strain9vto3x3t!(t, v)
  # t = (t + t')/2.0 # symmetrize
  # strain3x3tto9v!(v, t)
  # v6 = zeros(6)
  # strain9vto6v!(v6, v)
  # v9 = zeros(9)
  # strain6vto9v!(v9, v6)
  # @test norm(v-v9) < eps(1.0)

  # v = vec([1. 2 3 4 4 5 5 6 6])
  # v6 = zeros(6)
  # stress9vto6v!(v6, v)
  # v9 = zeros(9)
  # stress6vto9v!(v9, v6)
  # @test norm(v-v9) < eps(1.0)

end
end
using .mstressconversionm
mstressconversionm.test()

module mmmtdt1a
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
function test()
	MR = DeforModelRed3D
	symmet(a) = a + transpose(a)

	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(3, 3)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6

	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6

	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(3, 3))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
end
end
using .mmmtdt1a
mmmtdt1a.test()

module mmmtdt2
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Test
function test()
	MR = DeforModelRed2DStrain
	symmet(a) = a + transpose(a)

	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = rand(2, 2)
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	@test abs(det(a) - dett(MR, a)) / abs(det(a)) <= 1.0e-6

	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 6)
	strainttov!(MR, av, a)
	@test abs(det(a) - strainvdet(MR, av)) / abs(det(a)) <= 1.0e-6

	MR = DeforModelRed2DStrain
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 3)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 3)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
	a = symmet(rand(2, 2))
	av = fill(zero(eltype(a)), 3)
	strainttov!(MR, av, a)
	@test abs(tr(a) - strainvtr(MR, av)) / abs(tr(a)) <= 1.0e-6
end
end
using .mmmtdt2
mmmtdt2.test()

module mmtwistedeexportmm
using FinEtools
using FinEtoolsDeforLinear
using Test
using FinEtools.MeshExportModule
function test()
  E = 0.29e8;
  nu = 0.22;
  W = 1.1;
  L = 12.;
  t =  0.32;
  nl = 2; nt = 1; nw = 1; ref = 3;
  p =   1/W/t;
  #  Loading in the Z direction
  loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
  #   Loading in the Y direction
  #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
  tolerance  = t/1000;

  fens,fes  = H8block(L,W,t, nl*ref,nw*ref,nt*ref)

  # Reshape into a twisted beam shape
  for i = 1:count(fens)
    a = fens.xyz[i,1]/L*(pi/2); y = fens.xyz[i,2]-(W/2); z = fens.xyz[i,3]-(t/2);
    fens.xyz[i,:] = [fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
  end

  # Clamped end of the beam
  l1  = selectnode(fens; box = [0 0 -100*W 100*W -100*W 100*W], inflate  =  tolerance)
  e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
  e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
  e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

  # Traction on the opposite edge
  boundaryfes  =   meshboundary(fes);
  Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
  el1femm  = FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
  flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>loadv)

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

  # Make model data
  modeldata =  FDataDict(
  "fens"=> fens, "regions"=>  [region1],
  "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

  AE = AbaqusExporter("twisted_beam");
  HEADING(AE, "Twisted beam example");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(region1["femm"].integdomain.fes))
  ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdomain.fes), connasarray(flux1["femm"].integdomain.fes))
  NSET_NSET(AE, "l1", l1)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglass");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  SECTION_CONTROLS(AE, "section1", "HOURGLASS=ENHANCED")
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3)
  DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open("twisted_beam.inp") do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 223
  rm("twisted_beam.inp")

  true
end
end
using .mmtwistedeexportmm
mmtwistedeexportmm.test()

module mmtwistedeexport2mm
using FinEtools
using FinEtoolsDeforLinear
using Test
using FinEtools.MeshExportModule
function test()
  E = 0.29e8;
  nu = 0.22;
  W = 1.1;
  L = 12.;
  t =  0.32;
  nl = 2; nt = 1; nw = 1; ref = 3;
  p =   1/W/t;
  #  Loading in the Z direction
  loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
  #   Loading in the Y direction
  #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
  tolerance  = t/1000;

  fens,fes  = H8block(L,W,t, nl*ref,nw*ref,nt*ref)

  # Reshape into a twisted beam shape
  for i = 1:count(fens)
    a = fens.xyz[i,1]/L*(pi/2); y = fens.xyz[i,2]-(W/2); z = fens.xyz[i,3]-(t/2);
    fens.xyz[i,:] = [fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
  end

  # Clamped end of the beam
  l1  = selectnode(fens; box = [0 0 -100*W 100*W -100*W 100*W], inflate  =  tolerance)
  e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
  e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
  e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

  # Traction on the opposite edge
  boundaryfes  =   meshboundary(fes);
  Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
  el1femm  = FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
  flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>loadv)

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

  # Make model data
  modeldata =  FDataDict(
  "fens"=> fens, "regions"=>  [region1],
  "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

  AE = AbaqusExporter("twisted_beam");
  HEADING(AE, "Twisted beam example");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(region1["femm"].integdomain.fes))
  ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdomain.fes), connasarray(flux1["femm"].integdomain.fes))
  NSET_NSET(AE, "l1", l1)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglass");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  SECTION_CONTROLS(AE, "section1", "HOURGLASS=ENHANCED")
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1, 0.0)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2, 0.0)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3, 0.0)
  DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open("twisted_beam.inp") do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 223
  rm("twisted_beam.inp")

  true
end
end
using .mmtwistedeexport2mm
mmtwistedeexport2mm.test()

module mstresscomponentmap
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    MR = DeforModelRed1D
    @test stresscomponentmap(MR)[:x] == 1
    MR = DeforModelRed2DAxisymm
    @test stresscomponentmap(MR)[:x] == 1
    @test stresscomponentmap(MR)[:zz] == 3
end
end
using .mstresscomponentmap
mstresscomponentmap.test()

module mmMeasurement_3a
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  3.32;
    nl, nt, nw = 5, 3, 4;
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;

    # println("New segmentation fault?")
    for orientation in [:a :b :ca :cb]
        fens,fes  = T4block(L,W,t, nl,nw,nt, orientation)
        geom  =  NodalField(fens.xyz)
        MR  =  DeforModelRed3D
        material = MatDeforElastIso(MR, 0.0, Ea, nua, 0.0)

        femm  =  FEMMDeforLinearNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
        V = integratefunction(femm, geom, (x) ->  1.0)
        @test abs(V - W*L*t)/V < 1.0e-5
    end

end
end
using .mmMeasurement_3a
mmMeasurement_3a.test()

module mmvgradmat1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: infsup_sh
using Test
import LinearAlgebra: norm, cholesky
function test()
	Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = ( 1.0, 1.0, 1.0, 1, 1, 1, :a)
	Ea, nua, alphaa = ( 1.0, 0.3, 0.0)
	fens = FENodeSet(Float64[0     0     0
     0     3     3
     0     0     3
     3     0     3]);
	fes = FESetT4(reshape([1 2 3 4], 1, 4));

	MR  =  DeforModelRed3D

	# Property and material
	material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

	femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)

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

	Sh = infsup_sh(femm, geom, u);

	Sh1 = [0.5   0   0   0   0   0   -0.5   0   0   -0.0000   0   0
		   0    0.5   0   0   0   0   0   -0.5   0   0   -0.0000   0
		   0   0    0.5   0   0   0   0   0   -0.5   0   0   -0.0000
		   0   0   0    0.5   0   0   -0.5   0   0   0   0   0
		   0   0   0   0    0.5   0   0   -0.5   0   0   0   0
		   0   0   0   0   0    0.5   0   0   -0.5   0   0   0
		   -0.5   0   0   -0.5   0   0    1.5000   0   0   -0.5   0   0
		   0   -0.5   0   0   -0.5   0   0    1.5000   0   0   -0.5   0
		   0   0   -0.5   0   0   -0.5   0   0    1.5000   0   0   -0.5
		   -0.0000   0   0   0   0   0   -0.5   0   0    0.5   0   0
		   0   -0.0000   0   0   0   0   0   -0.5   0   0    0.5   0
		   0   0   -0.0000   0   0   0   0   0   -0.5   0   0    0.5]
	@test norm(Sh - Sh1)  <= 1.0e-9
end
end
using .mmvgradmat1
mmvgradmat1.test()

module mmdivmat2
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: infsup_gh
using Test
import LinearAlgebra: norm, cholesky
function test()
	Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = ( 1.0, 1.0, 1.0, 1, 1, 1, :a)
	Ea, nua, alphaa = ( 1.0, 0.3, 0.0)
	fens = FENodeSet(Float64[           0            0            0
							  -1.6200e-01   3.0000e+00   2.9791e+00
							            0            0   3.0000e+00
							   3.0000e+00   1.6200e-01   3.0000e+00]);
	fes = FESetT4(reshape([1 2 3 4], 1, 4));

	MR  =  DeforModelRed3D

	# Property and material
	material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

	femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)

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

	Gh = infsup_gh(femm, geom, u);

	Gh1 = [  7.0319e-08  -1.3022e-06  -1.8778e-04  -1.0111e-05   1.8724e-04            0  -1.7720e-04  -1.9604e-04   1.8778e-04   1.8724e-04   1.0111e-05            0
  -1.3022e-06   2.4115e-05   3.4774e-03   1.8724e-04  -3.4673e-03            0   3.2814e-03   3.6305e-03  -3.4774e-03  -3.4673e-03  -1.8724e-04            0
  -1.8778e-04   3.4774e-03   5.0146e-01   2.7000e-02  -5.0000e-01            0   4.7319e-01   5.2352e-01  -5.0146e-01  -5.0000e-01  -2.7000e-02            0
  -1.0111e-05   1.8724e-04   2.7000e-02   1.4538e-03  -2.6921e-02            0   2.5478e-02   2.8188e-02  -2.7000e-02  -2.6921e-02  -1.4538e-03            0
   1.8724e-04  -3.4673e-03  -5.0000e-01  -2.6921e-02   4.9855e-01            0  -4.7181e-01  -5.2200e-01   5.0000e-01   4.9855e-01   2.6921e-02            0
            0            0            0            0            0            0            0            0            0            0            0            0
  -1.7720e-04   3.2814e-03   4.7319e-01   2.5478e-02  -4.7181e-01            0   4.4651e-01   4.9401e-01  -4.7319e-01  -4.7181e-01  -2.5478e-02            0
  -1.9604e-04   3.6305e-03   5.2352e-01   2.8188e-02  -5.2200e-01            0   4.9401e-01   5.4656e-01  -5.2352e-01  -5.2200e-01  -2.8188e-02            0
   1.8778e-04  -3.4774e-03  -5.0146e-01  -2.7000e-02   5.0000e-01            0  -4.7319e-01  -5.2352e-01   5.0146e-01   5.0000e-01   2.7000e-02            0
   1.8724e-04  -3.4673e-03  -5.0000e-01  -2.6921e-02   4.9855e-01            0  -4.7181e-01  -5.2200e-01   5.0000e-01   4.9855e-01   2.6921e-02            0
   1.0111e-05  -1.8724e-04  -2.7000e-02  -1.4538e-03   2.6921e-02            0  -2.5478e-02  -2.8188e-02   2.7000e-02   2.6921e-02   1.4538e-03            0
            0            0            0            0            0            0            0            0            0            0            0            0]
    # @show Matrix(Gh)
	@test norm(Gh - Gh1) / norm(Gh1)  <= 2.0e-5
end
end
using .mmdivmat2
mmdivmat2.test()

module mmvgradmat2
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: infsup_sh
using Test
import LinearAlgebra: norm, cholesky
function test()
	Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = ( 1.0, 1.0, 1.0, 1, 1, 1, :a)
	Ea, nua, alphaa = ( 1.0, 0.3, 0.0)
	fens = FENodeSet(Float64[        0         0         0
								   -0.1620    3.0000    2.9791
								         0         0    3.0000
								    3.0000    0.1620    3.0000]);
	fes = FESetT4(reshape([1 2 3 4], 1, 4));

	MR  =  DeforModelRed3D

	# Property and material
	material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

	femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)

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

	Sh = infsup_sh(femm, geom, u);

	Sh1 = [5.0148e-01            0            0  -3.4774e-03            0            0  -4.9800e-01            0            0  -2.6321e-17            0            0
            0   5.0148e-01            0            0  -3.4774e-03            0            0  -4.9800e-01            0            0  -2.6321e-17            0
            0            0   5.0148e-01            0            0  -3.4774e-03            0            0  -4.9800e-01            0            0  -2.6321e-17
  -3.4774e-03            0            0   5.0000e-01            0            0  -4.9652e-01            0            0  -4.7374e-18            0            0
            0  -3.4774e-03            0            0   5.0000e-01            0            0  -4.9652e-01            0            0  -4.7374e-18            0
            0            0  -3.4774e-03            0            0   5.0000e-01            0            0  -4.9652e-01            0            0  -4.7374e-18
  -4.9800e-01            0            0  -4.9652e-01            0            0   1.4945e+00            0            0  -5.0000e-01            0            0
            0  -4.9800e-01            0            0  -4.9652e-01            0            0   1.4945e+00            0            0  -5.0000e-01            0
            0            0  -4.9800e-01            0            0  -4.9652e-01            0            0   1.4945e+00            0            0  -5.0000e-01
  -2.6321e-17            0            0  -4.7374e-18            0            0  -5.0000e-01            0            0   5.0000e-01            0            0
            0  -2.6321e-17            0            0  -4.7374e-18            0            0  -5.0000e-01            0            0   5.0000e-01            0
            0            0  -2.6321e-17            0            0  -4.7374e-18            0            0  -5.0000e-01            0            0   5.0000e-01]
	@test norm(Sh - Sh1) / norm(Sh1)  <= 1.0e-4
end
end
using .mmvgradmat2
mmvgradmat2.test()

module mmdivmat1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: infsup_gh
using Test
import LinearAlgebra: norm, cholesky
function test()
	Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = ( 1.0, 1.0, 1.0, 1, 1, 1, :a)
	Ea, nua, alphaa = ( 1.0, 0.3, 0.0)
	fens = FENodeSet(Float64[0     0     0
     0     3     3
     0     0     3
     3     0     3]);
	fes = FESetT4(reshape([1 2 3 4], 1, 4));

	MR  =  DeforModelRed3D

	# Property and material
	material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

	femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)

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

	Gh = infsup_gh(femm, geom, u);

	Gh1 = [0.0  0    0.0  0   -0.0  0    0.0    0.0   -0.0   -0.0  0  0
			0  0  0  0  0  0  0  0  0  0  0  0
			0.0  0    0.5  0   -0.5  0    0.5    0.5   -0.5   -0.5  0  0
			0  0  0  0  0  0  0  0  0  0  0  0
			-0.0  0   -0.5  0    0.5  0   -0.5   -0.5    0.5    0.5  0  0
			0  0  0  0  0  0  0  0  0  0  0  0
			0.0  0    0.5  0   -0.5  0    0.5    0.5   -0.5   -0.5  0  0
			0.0  0    0.5  0   -0.5  0    0.5    0.5   -0.5   -0.5  0  0
			-0.0  0   -0.5  0    0.5  0   -0.5   -0.5    0.5    0.5  0  0
			-0.0  0   -0.5  0    0.5  0   -0.5   -0.5    0.5    0.5  0  0
			0  0  0  0  0  0  0  0  0  0  0  0
			0  0  0  0  0  0  0  0  0  0  0  0]
	@test norm(Gh - Gh1)  <= 1.0e-9
end
end
using .mmdivmat1
mmdivmat1.test()

module mmmMergem1
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 12; nl  = 40; nc = 120;
    nh = 6; nl  = 20; nc = 60;
    # nh = 3; nl  = 8; nc = 30;
    tolerance = h/nh/100;

    t0 = time()
    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2.0*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
    # println("  before merging  = $(count(fens))")
    @test count(fens) == 8967
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    candidates = selectnode(fens, box = boundingbox([Rmed-h -Inf 0.0; Rmed+h +Inf 0.0]), inflate = tolerance)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);
    @test count(fens) == 8820

    # fens,fes = mergenodes(fens, fes,  tolerance);
    # println("  after merging  = $(count(fens))")
    #
    # println("Mesh generation ($(time() - t0) sec)")
    #
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
end
end
using .mmmMergem1
 mmmMergem1.test()

module mmmMergem2
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 12; nl  = 40; nc = 120;
    nh = 6; nl  = 20; nc = 60;
    # nh = 3; nl  = 8; nc = 30;
    tolerance = h/nh/100;

    t0 = time()
    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2.0*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
    # println("  before merging  = $(count(fens))")
    @test count(fens) == 8967
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], inflate = tolerance)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);
    @test count(fens) == 8820

    # fens,fes = mergenodes(fens, fes,  tolerance);
    # println("  after merging  = $(count(fens))")
    #
    # println("Mesh generation ($(time() - t0) sec)")
    #
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
end
end
using .mmmMergem2
 mmmMergem2.test()

module mmmMergem3
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 12; nl  = 40; nc = 120;
    nh = 6; nl  = 20; nc = 60;
    # nh = 3; nl  = 8; nc = 30;
    tolerance = h/nh/100;

    t0 = time()
    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2.0*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
    # println("  before merging  = $(count(fens))")
    @test count(fens) == 8967
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);
    @test count(fens) == 8820

    # fens,fes = mergenodes(fens, fes,  tolerance);
    # println("  after merging  = $(count(fens))")
    #
    # println("Mesh generation ($(time() - t0) sec)")
    #
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
end
end
using .mmmMergem3
 mmmMergem3.test()

 module mtwistedalgofin # From miscellaneous
 using FinEtools
 using FinEtoolsDeforLinear
 using FinEtoolsDeforLinear.AlgoDeforLinearModule
 using Statistics
 using LinearAlgebra: norm
 using Test
 function test()
     # println("""
     # The initially twisted cantilever beam is one of the standard test    problems for verifying the finite-element accuracy [1]. The beam is        clamped at one end and loaded either with unit in-plane or        unit out-of-plane force at the other. The centroidal axis of the beam s
     # straight at the undeformed  configuration, while its cross-sections are
     # twisted about the centroidal axis from 0 at the clamped end to pi/2 at
     # the free end.
     #
     # Reference:
     # Zupan D, Saje M (2004) On "A proposed standard set of problems to test
     # finite element accuracy": the twisted beam. Finite Elements in Analysis
     # and Design 40: 1445-1451.
     # """)
     E = 0.29e8;
     nu = 0.22;
     W = 1.1;
     L = 12.;
     t =  0.32;
     nl = 2; nt = 1; nw = 1; ref = 4;
     p =   1/W/t;
     #  Loading in the Z direction
     loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
     #   Loading in the Y direction
     #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
     tolerance  = t/1000;

     fens,fes  = H8block(L,W,t, nl*ref,nw*ref,nt*ref)

     # Reshape into a twisted beam shape
     for i = 1:count(fens)
         a = fens.xyz[i,1]/L*(pi/2); y = fens.xyz[i,2]-(W/2); z = fens.xyz[i,3]-(t/2);
         fens.xyz[i,:] = [fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
     end

     # Clamped end of the beam
     l1  = selectnode(fens; box = [0 0 -100*W 100*W -100*W 100*W], inflate  =  tolerance)
     e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
     e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
     e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

     # Traction on the opposite edge
     boundaryfes  =   meshboundary(fes);
     Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
     el1femm  = FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
     flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>ForceIntensity(loadv))

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

     # Make model data
     modeldata =  FDataDict(
     "fens"=> fens, "regions"=>  [region1],
     "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

     # Call the solver
     modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
     geom = modeldata["geom"]
     u = modeldata["u"]

     # Extract the solution
     nl = selectnode(fens, box = [L L -100*W 100*W -100*W 100*W],inflate = tolerance);
     theutip = mean(u.values[nl,:], dims = 1)
     # println("displacement  = $(theutip[dir]) as compared to converged $uex")
     # println("normalized displacement  = $(theutip[dir]/uex*100) %")
     @test (theutip[dir] - 0.005443006890338948) / 0.005443006890338948 <= 1.0e-6

     # Write out mesh with displacements
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam")
     modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

     # Write out mesh with stresses
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
     "quantity"=> :Cauchy, "component"=> :xy)
     modeldata = AlgoDeforLinearModule.exportstress(modeldata)

     # Write out mesh with stresses
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
     "quantity"=> :Cauchy, "component"=> :xz)
     modeldata = AlgoDeforLinearModule.exportstress(modeldata)

     # Write out mesh with von Mises stresses
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
     "quantity"=> :vm)
     modeldata = AlgoDeforLinearModule.exportstress(modeldata)

     # Write out mesh with von Mises stresses, elementwise
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-ew",
     "quantity"=> :vm)
     modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
     vm  = modeldata["postprocessing"]["exported"][1]["field"]
     # println("extremes of von Mises: $([minimum(vm.values),   maximum(vm.values)])")
     @test norm([minimum(vm.values),   maximum(vm.values)] - [5.21266051216496, 452.7873821727328]) <= 1.0e-3

     # Write out mesh with von Mises stresses, elementwise
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-ew",
     "quantity"=> :Cauchy, "component"=> :xz)
     modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

     # Write out mesh with principal stresses, elementwise
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-principal-1-ew",
     "quantity"=> :princCauchy, "component"=> 1)
     modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
     ps  = modeldata["postprocessing"]["exported"][1]["field"]
     # println("extremes of first principal stress: $([minimum(ps.values),   maximum(ps.values)])")
     @test norm([minimum(ps.values),   maximum(ps.values)] - [-5.9919354276389285, 461.8866079275564]) <= 1.0e-3

     # Write out mesh with principal stresses, elementwise
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-principal-3-ew",
     "quantity"=> :princCauchy, "component"=> 3)
     modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
     ps  = modeldata["postprocessing"]["exported"][1]["field"]
     # println("extremes of third principal stress: $([minimum(ps.values),   maximum(ps.values)])")
     @test norm([minimum(ps.values),   maximum(ps.values)] - [-461.88660792935354, 5.9919354278587855]) <= 1.0e-3

     # Write out mesh with principal stresses, elementwise
     modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-press-ew",
     "quantity"=> :pressure, "component"=> 1)
     modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
     ps  = modeldata["postprocessing"]["exported"][1]["field"]
     # println("extremes of pressure: $([minimum(ps.values),   maximum(ps.values)])")
     @test norm([minimum(ps.values),   maximum(ps.values)] - [-167.8998735447784, 167.89987354544647]) <= 1.0e-3
     try
         for f    in  ["twisted_beam-Cauchyxy-region 1.vtk",  "twisted_beam-press-ew-pressure1-region 1.vtk",
             "twisted_beam-Cauchyxz-region 1.vtk",     "twisted_beam-principal-1-ew-princCauchy1-region 1.vtk",
             "twisted_beam-ew-Cauchyxz-region 1.vtk",  "twisted_beam-principal-3-ew-princCauchy3-region 1.vtk",
             "twisted_beam1.vtk",  "twisted_beam-ew-vm1-region 1.vtk",       "twisted_beam-vm1-region 1.vtk"]
             rm(f)
         end
     catch end

     # println("Done")
     true
 end
 end
 using .mtwistedalgofin
 mtwistedalgofin.test()


text:
"""
    transpose_to_axes(coord, b, n)

Convert from Hilbert index to a coordinatefor `b` bits
in `n` dimensions.
"""
function transpose_to_axes!(X::Vector{T}, b, n) where {T <: Integer}
    N = T(2) << (b - 1)
    # Gray decode by H^(H/2)
    t = X[n] >> 1
    for io = n:-1:2
        X[io] ⊻= X[io - 1]
    end
    X[1] ⊻= t
    # Undo excess work
    Q = T(2)
    while Q != N
        P = Q - one(T)
        for jo = n:-1:1
            if (X[jo] & Q) != zero(T)
                X[1] ⊻= P  # invert
            else  # exchange
                t = (X[1] ⊻ X[jo]) & P
                X[1] ⊻= t
                X[jo] ⊻= t
            end
        end
        Q <<= 1
    end
end

function axes_to_transpose!(X::Vector{T}, b, n) where {T <: Integer}
    M = one(T) << (b - 1)
    # Inverse undo
    Q = M
    while Q > one(T)
        P = Q - one(T)
        for io = 1:n
            if (X[io] & Q) != zero(T)
                X[1] ⊻= P
            else
                t = (X[1] ⊻ X[io]) & P
                X[1] ⊻= t
                X[io] ⊻= t
            end
        end
        Q >>= 1
    end
    # Gray encode
    for jo = 2:n
        X[jo] ⊻= X[jo - 1]
    end
    t2 = zero(T)
    Q = M
    while Q > one(T)
        if (X[n] & Q) != 0
            t2 ⊻= (Q - one(T))
        end
        Q >>= one(T)
    end
    for ko = 1:n
        X[ko] ⊻= t2
    end
end

"""
Takes a vector of length `n` and places the bits of all `n` integers
into a single integer. The vector's 1st component is the most significant bit.
"""
function interleave_transpose(::Type{T}, X::Vector, b, n) where {T}
    h = zero(T)
    for i in 0:(b - 1)
        for d in 1:n
            ith_bit = (X[d] & (1<<i)) >> i
            h |= ith_bit << (i * n + n - d)
        end
    end
    h
end

"""
Takes a single integer and places its values into components of a vector,
bit-by-bit.
"""
function outerleave_transpose!(X::Vector{T}, h, b, n) where {T <: Integer}
    X .= zero(T)
    for i in 0:(b-1)
        for d in 1:n
            ith_bit = (h & (one(h) << (i * n + n - d))) >> (i * n + n - d)
            X[d] |= ith_bit << i
        end
    end
end

"""
Takes a vector of length `n` and places the bits of all `n` integers
into a single integer. The vector's 1st component is the least significant bit.
"""
function interleave_transpose_low(::Type{T}, X::Vector{T}, b, n) where {T}
    h = zero(T)
    for i in 0:(b - 1)
        for d in 1:n
            h |= ((X[d] & (one(T)<<i))) << (i*(n - 1) + d - 1)
        end
    end
    h
end

function outerleave_transpose_low!(X::Vector{T}, h, b, n) where {T <: Integer}
    X .= zero(T)
    for i in 0:(b-1)
        for d in 1:n
            X[d] |= (h & (one(T) << (i * n + d - 1))) >> (i * (n - 1) + d - 1)
        end
    end
end

"""
    GlobalGray(b, n)
    GlobalGray(T, b, n)

`T` is a data type for the Hilbert index. It can be signed or unsigned,
as long as it has at least `n * b` bits. `n` is the number of dimensions,
and `b` is the bits per dimension, so each axis value should be between
0 and ``2^b - 1``, inclusive, for the zero-based interface. They should be
between 1 and ``2^b``, inclusive, for the one-based interface.

The GlobalGray algorithm is an n-dimensional Hilbert curve with a simplified
implementation. It follows an article, "Programming the Hilbert Curve," by
John Skilling, 707 (2004), http://dx.doi.org/10.1063/1.1751381.
I call it "Global Gray" because the insight of the article
is that a single, global Gray code can be applied to all
np bits of a Hilbert length.
"""
struct GlobalGray{T} <: HilbertAlgorithm{T}
    b::Int
    n::Int
end

axis_type(gg::GlobalGray) = large_enough_unsigned(gg.b)

function GlobalGray(b, n)
    ttype = large_enough_unsigned(b * n)
    GlobalGray{ttype}(b, n)
end

function GlobalGray(::Type{T}, b, n) where {T}
    GlobalGray{T}(b, n)
end

function encode_hilbert_zero!(g::GlobalGray{T}, X::Vector)::T where {T}
    axes_to_transpose!(X, g.b, g.n)
    interleave_transpose(T, X, g.b, g.n)
end

"""
    encode_hilbert_zero(ha::HilbertAlgorithm{T}, X::Vector{A})

Takes an n-dimensional vector `X` and returns a single integer of type
`T` which orders `X` to improve spatial locality. The input `X` has multiple
axes and the output is called a Hilbert index. This version is zero-based,
so each axis counts from 0, and the smallest Hilbert index is 0.
"""
function encode_hilbert_zero(g::GlobalGray{T}, X::Vector)::T where {T}
    Y = copy(X)
    encode_hilbert_zero!(g, Y)
end

"""
    decode_hilbert_zero!(ha::HilbertAlgorithm{T}}, X::Vector{A}, h::T)

Given a Hilbert index, `h`, computes an n-dimensional coordinate `X`. The type of
the Hilbert index is large enought to contain the bits of all dimensions of the
axis vector, `X`.
"""
function decode_hilbert_zero!(g::GlobalGray{T}, X::Vector, h::T) where {T}
    outerleave_transpose!(X, h, g.b, g.n)
    transpose_to_axes!(X, g.b, g.n)
end

ext: jl


text:
module KuboOscillator

    using GeometricEquations
    using Parameters

    export kubo_oscillator_sde_1, kubo_oscillator_psde_1, kubo_oscillator_spsde_1
    export kubo_oscillator_sde_2, kubo_oscillator_psde_2, kubo_oscillator_spsde_2
    export kubo_oscillator_sde_3, kubo_oscillator_psde_3, kubo_oscillator_spsde_3
    export kubo_oscillator_ode

    q_init_A=[0.5, 0.0]
    q_init_B=[[ 0.5, 0.0],
              [ 0.0, 0.5],
              [-0.5, 0.0]]

    const noise_intensity = 0.1
    const p = (ν = noise_intensity,)

    const Δt = 0.01
    const nt = 10

    function kubo_oscillator_sde_v(t, q, v_out, params)
        v_out[1]=  q[2]
        v_out[2]= -q[1]
    end

    function kubo_oscillator_sde_B(t, q, B_out::AbstractVector, params)
        @unpack ν = params
        B_out[1] = +ν*q[2]
        B_out[2] = -ν*q[1]
    end

    function kubo_oscillator_sde_B(t, q, B_out::AbstractMatrix, params, col=0)
        @unpack ν = params
        if col == 0
            B_out[1,1] = +ν*q[2]
            B_out[2,1] = -ν*q[1]
        else
            B_out[1,col] = +ν*q[2]
            B_out[2,col] = -ν*q[1]
        end
    end

    function kubo_oscillator_sde_1(q₀=q_init_A, params=p)
        # q_init_A - interpreted as one random initial conditions with one sample path
        # 1-dimensional noise
        SDE(1, 1, kubo_oscillator_sde_v, kubo_oscillator_sde_B, q₀; parameters=params)
    end

    function kubo_oscillator_sde_2(q₀=q_init_A, params=p)
        # q_init_A - single deterministic initial condition
        # Generating 3 sample paths
        # 1-dimensional noise
        SDE(1, 3, kubo_oscillator_sde_v, kubo_oscillator_sde_B, q₀; parameters=params)
    end

    function kubo_oscillator_sde_3(q₀=q_init_B, params=p)
        # q_init_B - interpreted as three random initial conditions
        # The 3 columns correspond to 3 sample paths
        # 1-dimensional noise
        SDE(1, 1, kubo_oscillator_sde_v, kubo_oscillator_sde_B, q₀; parameters=params)
    end

    # ODE

    function kubo_oscillator_ode(q₀=q_init_A, params=p)
        ODE(kubo_oscillator_sde_v, q₀; parameters=params)
    end

    # PSDE

    q_init_C=[0.5]
    p_init_C=[0.0]

    q_init_D=[[0.5], [0.0], [-0.5]]
    p_init_D=[[0.0], [0.5], [ 0.0]]

    function kubo_oscillator_psde_v(t, q, p, v_out, params)
        v_out[1] =  p[1]
    end

    function kubo_oscillator_psde_f(t, q, p, f_out, params)
        f_out[1] = -q[1]
    end

    function kubo_oscillator_psde_B(t, q, p, B_out, params)
        @unpack ν = params
        B_out[1,1] = +ν*p[1]
    end

    function kubo_oscillator_psde_G(t, q, p, G_out, params)
        @unpack ν = params
        G_out[1,1] = -ν*q[1]
    end

    function kubo_oscillator_psde_1(q₀=q_init_C, p₀=p_init_C, params=p)
        # q_init_C - interpreted as a single random initial condition with one sample path
        # 1-dimensional noise
        PSDE(1, 1, kubo_oscillator_psde_v, kubo_oscillator_psde_f,
                   kubo_oscillator_psde_B, kubo_oscillator_psde_G,
                   q₀, p₀; parameters=params)
    end

    function kubo_oscillator_psde_2(q₀=q_init_C, p₀=p_init_C, params=p)
        # q_init_C - single deterministic initial condition
        # Generating 3 sample paths
        # 1-dimensional noise
        PSDE(1, 3, kubo_oscillator_psde_v, kubo_oscillator_psde_f,
                   kubo_oscillator_psde_B, kubo_oscillator_psde_G,
                   q₀, p₀; parameters=params)
    end

    function kubo_oscillator_psde_3(q₀=q_init_D, p₀=p_init_D, params=p)
        # q_init_D - interpreted as a single random initial condition
        # The 3 columns correspond to 3 sample paths
        # 1-dimensional noise
        PSDE(1, 1, kubo_oscillator_psde_v, kubo_oscillator_psde_f,
                   kubo_oscillator_psde_B, kubo_oscillator_psde_G,
                   q₀, p₀; parameters=params)
    end

    # SPSDE

    function kubo_oscillator_spsde_v(t, q, p, v_out, params)
        v_out[1] =  p[1]
    end

    function kubo_oscillator_spsde_f1(t, q, p, f_out, params)
        f_out[1] = -q[1]
    end

    function kubo_oscillator_spsde_f2(t, q, p, f_out, params)
        f_out[1] = 0
    end

    function kubo_oscillator_spsde_B(t, q, p, B_out, params)
        @unpack ν = params
        B_out[1,1] = +ν*p[1]
    end

    function kubo_oscillator_spsde_G1(t, q, p, G_out, params)
        @unpack ν = params
        G_out[1,1] = -ν*q[1]
    end

    function kubo_oscillator_spsde_G2(t, q, p, G_out, params)
        G_out[1,1] = 0
    end

    function kubo_oscillator_spsde_1(q₀=q_init_C, p₀=p_init_C, params=p)
        # q_init_C - interpreted as a single random initial condition with one sample path
        # 1-dimensional noise
        SPSDE(1, 1, kubo_oscillator_spsde_v, kubo_oscillator_spsde_f1, kubo_oscillator_spsde_f2,
                    kubo_oscillator_spsde_B, kubo_oscillator_spsde_G1, kubo_oscillator_spsde_G2,
                    q₀, p₀; parameters=params)
    end

    function kubo_oscillator_spsde_2(q₀=q_init_C, p₀=p_init_C, params=p)
        # q_init_C - single deterministic initial condition
        # Generating 3 sample paths
        # 1-dimensional noise
        SPSDE(1, 3, kubo_oscillator_spsde_v, kubo_oscillator_spsde_f1, kubo_oscillator_spsde_f2,
                    kubo_oscillator_spsde_B, kubo_oscillator_spsde_G1, kubo_oscillator_spsde_G2,
                    q₀, p₀; parameters=params)
    end

    function kubo_oscillator_spsde_3(q₀=q_init_D, p₀=p_init_D, params=p)
        # q_init_D - interpreted as a single random initial condition
        # The 3 columns correspond to 3 sample paths
        # 1-dimensional noise
        SPSDE(1, 1, kubo_oscillator_spsde_v, kubo_oscillator_spsde_f1, kubo_oscillator_spsde_f2,
                    kubo_oscillator_spsde_B, kubo_oscillator_spsde_G1, kubo_oscillator_spsde_G2,
                    q₀, p₀; parameters=params)
    end

end


text:
#!/usr/bin/env julia

nums = [1,2,3,4]

for i in nums
	println(i)
end

function myfun(x)
	println("myfun:",x)
end

function anotherfun(y)
	println("another fun: ",y)
end

function thirdfun(z)
	println("this is z: ",z)
end

# Can put functions in an array
analyses = [myfun,anotherfun,thirdfun]
for f in analyses
	f(99)
end

push!(nums,22)

println(nums)

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function createdataranges(start::Point{Int}, data, incols::Bool=true)
    tbldata = Matrix[]
    tblrngs = CellRange{Int}[]
    nelements = incols ? size(data[1],2) : size(data[1],1)
    kd = length(data)

    if incols
        numRows = size(data[1],1)
        for col=1:nelements
            for k=1:kd
                el = data[k]
                j = (col-1)*kd + k
                statrow = CellRange(
                    Box{Int}(Point(start.x + j - 1, start.y),
                             Point(start.x + j, start.y + numRows)),
                j+5,k+5)
                push!(tblrngs, statrow)
                push!(tbldata, reshape(el[:, col], numRows, 1))
            end
        end
    else
        numCols = size(data[1],2)
        for row=1:nelements
            for k=1:kd
                el = data[k]
                i = (row-1)*kd + k
                statrow = CellRange(
                    Box{Int}(Point(start.x, start.y + i - 1),
                             Point(start.x + numCols, start.y + i)),
                    i+5,k+5)
                push!(tblrngs, statrow)
                push!(tbldata, el[row,:])
            end
        end
    end
    return tbldata, tblrngs
end

const EMPTY = 1
const COLNAME = 2
const ROWNAME = 3
const COLLEVELS = 4
const ROWLEVELS = 5

checkstylelen(styles, p::Pivot) = (length(p.data) + 5 == length(styles)) ||
  error("length of data + 3 != length of styles")

function checkdatasize(p::Pivot)
    res = [size(el) for el in p.data]
    n = length(res)
    if n > 0
        el = res[1]
        for i=2:n
            if el != res[i]
                error("data sizes are not all the same\n$res")
            end
        end
    end
    return true
end

function totable(p::Pivot)
    styles = pivotstyles(p)
    checkstylelen(styles, p)
    checkdatasize(p)

    pm = p.metadata
    collevels = CartesianDataElement(colLevels(pm), true)
    rowlevels = CartesianDataElement(rowLevels(pm), false)
    rownames = rowNames(pm)
    rownames = reshape(rownames, 1, length(rownames))
    colnames = colNames(pm)

    # figure out where to place the tables
    rowRibbonHeight, rowRibbonWidth  = size(rowlevels)
    colRibbonHeight, colRibbonWidth = size(collevels)
    numCols = rowRibbonWidth + colRibbonWidth
    numRows = rowRibbonHeight + colRibbonHeight + 1

    emptydata = ExpandDataElement("", colRibbonHeight, rowRibbonWidth)
    colnames =  ExpandDataElement(join(colnames, " / "), 1, colRibbonWidth)

    emptyrng = CellRange(
        Box{Int}(Point(0,0), Point(rowRibbonWidth, colRibbonHeight)),
        EMPTY,EMPTY)
    colnamerng = CellRange(
        Box{Int}(Point(rowRibbonWidth, 0), Point(numCols, 1)),
        COLNAME,COLNAME)
    rownamerng = CellRange(
        Box{Int}(Point(0, colRibbonHeight), Point(rowRibbonWidth, colRibbonHeight + 1)),
        ROWNAME,ROWNAME)
    collevelrng = CellRange(
        Box{Int}(Point(rowRibbonWidth, 1), Point(numCols, 1+colRibbonHeight)),
        COLLEVELS,COLLEVELS)
    rowlevelrng = CellRange(
        Box{Int}(Point(0, 1+colRibbonHeight), Point(rowRibbonWidth, numRows)),
        ROWLEVELS,ROWLEVELS)

    data = Any[emptydata, colnames, rownames, collevels, rowlevels]
    rngs = CellRange{Int}[emptyrng; colnamerng; rownamerng; collevelrng; rowlevelrng]

    ndata, nrngs = createdataranges(Point(rowRibbonWidth, colRibbonHeight + 1),
                                    p.data, incols(p.metadata))
    append!(data, ndata)
    append!(rngs, nrngs)

    return Table(data, Layout(rngs), styles)
end

function strings(xs)
    n = length(xs)
    res = fill(String[], n)
    for i=1:n
        res[i] = map(string, xs[i])
    end
    return res
end

function levellength(xs)
    if length(xs) > 0
        return prod([length(x) for x in xs])
    else
        return 1
    end
end

@doc """
Args:
  df:           DataFrame to summarize
  rowvars:      row variables in table
  colvars:      column variables in table
  stats:        aggregate statistics to decribe each combination
                of row and column variables
  statnames:    statistic names used for table
  varorders:    makes the levels of the row and column variables
                follow the order given
  stat_in_cols: if true then the statistics are considered column
                levels of the "stats" variable. Otherwise they are
                row levels.

Usage:

```
using RDatasets
using Pivots
using Tables.HTMLFormat
const tb = Tables.HTMLFormat

house = RDatasets.dataset("Ecdat", "Housing")
p = summarize(
    house,
    rowvars      = [:Bedrooms],
    colvars      = Symbol[],
    statnames    = String["mean(price)", "std(price)", "count"],
    stats        = [df->mean(df[:Price]), df->std(df[:Price]), df->length(df[:Price])],
    varorders    = [:Bedrooms => Int32[1:6], :Driveway => ["no", "yes"]],
    stat_in_cols = true)

tbl = totable(p)
tb.HTML(tbl)
```
""" ->
function summarize{T <: String}(
    df::AbstractDataFrame;
    rowvars::Vector{Symbol} = Symbol[],
    colvars::Vector{Symbol} = Symbol[],
    stats::Vector{Function} = Function[],
    statnames::Vector{T} = fill("", length(stats)),
    varorders=Dict{Symbol,Any}(),
    stat_in_cols=true)

    @assert length(statnames) == length(stats)

    levels_of_rowvars = [haskey(varorders, rowvar) ? varorders[rowvar] :
                           unique(df[rowvar]) for rowvar in rowvars]
    levels_of_colvars = [haskey(varorders, colvar) ? varorders[colvar] :
                           unique(df[colvar]) for colvar in colvars]

    rownames = String[string(rowvar) for rowvar in rowvars]
    colnames = String[string(colvar) for colvar in colvars]

    pm = PivotMetadata(rownames, colnames,
                       strings(levels_of_rowvars),
                       strings(levels_of_colvars), statnames, stat_in_cols)

    statlen = length(stats)

    # determine how big the statistic matrices are
    n = levellength(levels_of_rowvars)
    m = levellength(levels_of_colvars)

    # create the global statistics
    frsts = [stat(df[1,:]) for stat in stats]
    types = [typeof(frst) for frst in frsts]

    # make the statistic matrices
    results = [fill(Nullable{tp}(), n, m) for tp in types]

    allvars = [colvars; rowvars]
    if length(allvars) > 0
        for (i,grp) in groupby(df, allvars, varorders)
            for k=1:statlen
                results[k][i] = Nullable{types[k]}(stats[k](grp))
            end
        end
    else
        for k=1:length(types)
            results[k][1,1] = Nullable{types[k]}(statistics[k](df))
        end
    end

    p = Pivot(results, pm)

    return p
end

max_stars_repo_name: cpritcha/Pivots

max_issues_repo_name: cpritcha/Pivots

max_forks_repo_name: cpritcha/Pivots


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

"""
    Some{T}

A wrapper type used in `Union{Some{T}, Nothing}` to distinguish between the absence
of a value ([`nothing`](@ref)) and the presence of a `nothing` value (i.e. `Some(nothing)`).

Use [`something`](@ref) to access the value wrapped by a `Some` object.
"""
struct Some{T}
    value::T
end

promote_rule(::Type{Some{T}}, ::Type{Some{S}}) where {T, S<:T} = Some{T}

nonnothingtype(::Type{T}) where {T} = Core.Compiler.typesubtract(T, Nothing)
promote_rule(T::Type{Nothing}, S::Type) = Union{S, Nothing}
function promote_rule(T::Type{>:Nothing}, S::Type)
    R = nonnothingtype(T)
    R >: T && return Any
    T = R
    R = promote_type(T, S)
    return Union{R, Nothing}
end

function nonnothingtype_checked(T::Type)
    R = nonnothingtype(T)
    R >: T && error("could not compute non-nothing type")
    return R
end

convert(::Type{T}, x::T) where {T>:Nothing} = x
convert(::Type{T}, x) where {T>:Nothing} = convert(nonnothingtype_checked(T), x)
convert(::Type{Some{T}}, x::Some{T}) where {T} = x
convert(::Type{Some{T}}, x::Some) where {T} = Some{T}(convert(T, x.value))

function show(io::IO, x::Some)
    if get(io, :typeinfo, Any) == typeof(x)
        show(io, x.value)
    else
        print(io, "Some(")
        show(io, x.value)
        print(io, ')')
    end
end

"""
    notnothing(x)

Throw an error if `x === nothing`, and return `x` if not.
"""
notnothing(x::Any) = x
notnothing(::Nothing) = throw(ArgumentError("nothing passed to notnothing"))

"""
    isnothing(x)

Return `true` if `x === nothing`, and return `false` if not.

!!! compat "Julia 1.1"
    This function requires at least Julia 1.1.
"""
isnothing(::Any) = false
isnothing(::Nothing) = true

"""
    something(x, y...)

Return the first value in the arguments which is not equal to [`nothing`](@ref),
if any. Otherwise throw an error.
Arguments of type [`Some`](@ref) are unwrapped.

See also [`coalesce`](@ref).

# Examples
```jldoctest
julia> something(nothing, 1)
1

julia> something(Some(1), nothing)
1

julia> something(missing, nothing)
missing

julia> something(nothing, nothing)
ERROR: ArgumentError: No value arguments present
```
"""
function something end

something() = throw(ArgumentError("No value arguments present"))
something(x::Nothing, y...) = something(y...)
something(x::Some, y...) = x.value
something(x::Any, y...) = x

max_stars_repo_name: sloede/julia

max_issues_repo_name: sloede/julia

max_forks_repo_name: sloede/julia


text:
using Documenter
using ImplicitDifferentiation
using Literate

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

EXAMPLES_DIR_JL = joinpath(dirname(@__DIR__), "test")
EXAMPLES_DIR_MD = joinpath(@__DIR__, "src", "examples")

for file in readdir(EXAMPLES_DIR_MD)
    if endswith(file, ".md")
        rm(joinpath(EXAMPLES_DIR_MD, file))
    end
end

for file in readdir(EXAMPLES_DIR_JL)
    if (!endswith(file, ".jl")) ||
        startswith(file, "runtests") ||
        startswith(file, "profiling")
        continue
    else
        Literate.markdown(
            joinpath(EXAMPLES_DIR_JL, file),
            EXAMPLES_DIR_MD;
            documenter=true,
            flavor=Literate.DocumenterFlavor(),
        )
    end
end

function markdown_title(path)
    title = "?"
    open(path, "r") do file
        for line in eachline(file)
            if startswith(line, '#')
                title = strip(line, [' ', '#'])
                break
            end
        end
    end
    return String(title)
end

pages = [
    "Home" => "index.md",
    "Mathematical background" => "background.md",
    "API reference" => "api.md",
    "Examples" => [
        markdown_title(joinpath(EXAMPLES_DIR_MD, file)) => joinpath("examples", file)
        for file in sort(readdir(EXAMPLES_DIR_MD)) if endswith(file, ".md")
    ],
]

makedocs(;
    modules=[ImplicitDifferentiation],
    authors="Guillaume Dalle, Mohamed Tarek and contributors",
    repo="https://github.com/gdalle/ImplicitDifferentiation.jl/blob/{commit}{path}#{line}",
    sitename="ImplicitDifferentiation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gdalle.github.io/ImplicitDifferentiation.jl",
        assets=String[],
        edit_branch="main",
    ),
    pages=pages,
)

deploydocs(; repo="github.com/gdalle/ImplicitDifferentiation.jl", devbranch="main")


text:
using Documenter

try
    using Nash
catch
    if !("../src/" in LOAD_PATH)
       push!(LOAD_PATH,"../src/")
       @info "Added \"../src/\"to the path: $LOAD_PATH "
       using Nash
    end
end

makedocs(
    sitename = "Nash.jl",
    format = format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [Nash],
    pages = ["index.md", "reference.md"],
    doctest = true
)

deploydocs(
    repo ="github.com/KrainskiL/Nash.jl.git",
    target="build"
)


text:
using Gridap
using GridapGmsh

model = GmshDiscreteModel("demo/demo.msh")

order = 1
dirichlet_tags = ["boundary1","boundary2"]
u_boundary1(x) = 0.0
u_boundary2(x) = 1.0

reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe,dirichlet_tags=dirichlet_tags)
U = TrialFESpace(V,[u_boundary1,u_boundary2])

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

a(u,v) = ∫( ∇(u)⋅∇(v) )dΩ
l(v) = 0

op = AffineFEOperator(a,l,U,V)
uh = solve(op)

writevtk(Ω,"demo",cellfields=["uh"=>uh])


text:
network1 = @min_reaction_network rnm begin
    2.0, X + Y --> XY
    1.0, XY ← X + Y
    (1,X), Z ↔ Z1 + Z2
end
addsdes!(network1)
addjumps!(network1)

@test length(network1.f_func) == 6
@test typeof(network1.f) <: Function
@test length(network1.g_func) == 24
@test typeof(network1.g) <: Function
@test size(network1.p_matrix) == (6,4)
@test length(network1.params) == 0
@test length(network1.syms) == 6
@test typeof(network1.jumps[1]) <:ConstantRateJump
@test length(network1.jump_rate_expr) == 4
@test length(network1.jump_affect_expr) == 4

network2 = @min_reaction_network rnm begin
    1.0, ∅ > X1
    X, 0 → X1
    (1.0,X1,X2), (X1,X2,X3) --> (X1,X2,X3)
    (1.0,X2,X1), (X2,X1,X3) --> X1
    ((1.0,2,3.), (P1,X2,X1)), X2 ↔ (X2+X4,X1,X2)
    hill(P2,2.,1,1), asdasdrqr < X1
    mm(X1,2.,1), X1 ⟾ X5
end
addsdes!(network2)
addjumps!(network2)

@test length(network2.f_func) == 6
@test typeof(network2.f) <: Function
@test length(network2.g_func) == 96
@test typeof(network2.g) <: Function
@test size(network2.p_matrix) == (6,16)
@test length(network2.params) == 0
@test length(network2.syms) == 6
@test typeof(network2.jumps[1]) <:ConstantRateJump
@test length(network2.jump_rate_expr) == 16
@test length(network2.jump_affect_expr) == 16

network3 = @min_reaction_network rnm2 begin
    p3, X + Y --> XY
    p2, XY ← X + Y
    (p1,X), Z ↔ Z1 + Z2
end p1 p2 p3

@test length(network3.params) == 3
@test typeof(network3) == rnm2
@test typeof(network1) == typeof(network2)

network4 = @min_reaction_network rnm begin
    1, X + Y > XY
    p1, X + Y < XY
    2, Z --> 0
end p1
addsdes!(network4)
addjumps!(network4)
network5 = @min_reaction_network rnm begin
    (1,p1), X + Y ↔ XY
    2Z, ∅ ⟽ Z
end p1
addsdes!(network5)
addjumps!(network5)

@test network4.f_func == network5.f_func
@test network4.g_func == network5.g_func
@test network4.p_matrix == network5.p_matrix
@test network4.params == network5.params
@test network4.jump_rate_expr == network5.jump_rate_expr
@test network4.jump_affect_expr == network5.jump_affect_expr

for i = 1:100
    u = 5*rand(4)
    du4 = 3*rand(4); du5 = du4;
    du4g = 2.5*rand(4,3); du5g = du4g;
    t = 9*rand(1)[1]
    p12 = 2*rand(1)[1]; p =[p12]

    @test network4.f(du4,u,p,t) == network5.f(du5,u,p,t)
    @test network4.g(du4g,u,p,t) == network5.g(du5g,u,p,t)
    @test network4.jumps[1].rate(u,p,t) == network5.jumps[1].rate(u,p,t)
    @test network4.jumps[2].rate(u,p,t) == network5.jumps[2].rate(u,p,t)
    @test network4.jumps[3].rate(u,p,t) == network5.jumps[3].rate(u,p,t)
end

network6 = @min_reaction_network begin
    p2, X + Y --> XY
    p1, XY ← X + Y
    (1,X), Z ↔ Z1 + Z2
end p1 p2
addsdes!(network6)

@test length(network6.f_func) == 6
@test length(network6.g_func) == 24
@test size(network6.p_matrix) == (6,4)
@test length(network6.params) == 2
@test length(network6.syms) == 6
@test typeof(network6) <: DiffEqBase.AbstractReactionNetwork
@test typeof(network6) == min_reaction_network


text:
## Example for using LSTNet

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Flux, Plots, BSON
include("../shared/Sequentialize.jl")
include("../data/dataloader.jl")
include("DARNN.jl")

# Load some sample data
poollength = 10
horizon = 6
datalength = 500
input, target = get_data(:solar, poollength, datalength, horizon)

# Define the network architecture
inputsize = size(input,1)
encodersize = 10
decodersize = 10

# Define the neural net
model = DARNN(inputsize, encodersize, decodersize, poollength, 1)

# MSE loss
function loss(x,y)
  #Flux.reset!(model)
  return sum(abs2,model(x) - y')
end

# Callback for plotting the training
cb = function()
  Flux.reset!(model)
  pred = model(input)'
  Flux.reset!(model)
  p1=plot(pred, label="Predict")
  p1=plot!(target, label="Data")
  display(plot(p1))
end

# Training loop
Flux.train!(loss, Flux.params(model),Iterators.repeated((input, target),10),
            ADAM(0.007), cb=cb)

max_stars_repo_name: cbdavis/FluxArchitectures

max_issues_repo_name: cbdavis/FluxArchitectures

max_forks_repo_name: cbdavis/FluxArchitectures


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Generic function to compose rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export compose_rotation

################################################################################
#                              Compose Rotations
################################################################################

"""
    compose_rotation(R1, [, R2, R3, R4, R5, ...])

Compute a composed rotation using the rotations `R1`, `R2`, `R3`, `R4`, ..., in
the following order:

     First rotation
     |
     |
    R1 => R2 => R3 => R4 => ...
           |
           |
           Second rotation

The rotations can be described by:

- A direction cosine matrix ([`DCM`](@ref));
- An Euler angle and axis ([`EulerAngleAxis`](@ref));
- A set of Euler angles ([`EulerAngles`](@ref)); or
- A quaternion ([`Quaternion`](@ref)).

Notice, however, that all rotations **must be** of the same type (DCM or
quaternion).

The output will have the same type as the inputs.

# Example

```jldoctest
julia> D1 = angle_to_dcm(pi / 3, pi / 4, pi / 5, :ZYX);

julia> D2 = angle_to_dcm(-pi / 5, -pi / 4, -pi / 3, :XYZ);

julia> compose_rotation(D1, D2)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.0          1.08801e-17  3.54837e-17
 1.08801e-17  1.0          2.88714e-17
 3.54837e-17  2.88714e-17  1.0

julia> ea1 = EulerAngleAxis(30 * pi / 180, [0, 1, 0]);

julia> ea2 = EulerAngleAxis(45 * pi / 180, [0, 1, 0]);

julia> compose_rotation(ea1, ea2)
EulerAngleAxis{Float64}:
  Euler angle : 1.309 rad  (75.0°)
  Euler axis  : [0.0, 1.0, 0.0]

julia> Θ1 = EulerAngles(1, 2, 3, :ZYX);

julia> Θ2 = EulerAngles(-3, -2, -1, :XYZ);

julia> compose_rotation(Θ1, Θ2)
EulerAngles{Float64}:
  R(X) : -1.66533e-16 rad  (-9.54166e-15°)
  R(Y) :  9.24446e-33 rad  ( 5.29669e-31°)
  R(Z) : -1.11022e-16 rad  (-6.36111e-15°)

julia> q1 = angle_to_quat(pi / 3, pi / 4, pi / 5, :ZYX);

julia> q2 = angle_to_quat(-pi / 5, -pi / 4, -pi / 3, :XYZ);

julia> compose_rotation(q1, q2)
Quaternion{Float64}:
  + 1.0 + 0.0⋅i + 2.08167e-17⋅j + 5.55112e-17⋅k
```
"""
@inline compose_rotation(D::DCM) = D
@inline compose_rotation(D::DCM, Ds::DCM...) = compose_rotation(Ds...)*D

@inline compose_rotation(ea::EulerAngleAxis) = ea
@inline function compose_rotation(ea::EulerAngleAxis, eas::EulerAngleAxis...)
    return compose_rotation(eas...)*ea
end

@inline compose_rotation(Θ::EulerAngles) = Θ
@inline compose_rotation(Θ::EulerAngles, Θs::EulerAngles...) = compose_rotation(Θs...)*Θ

@inline compose_rotation(q::Quaternion) = q
@inline compose_rotation(q::Quaternion, qs::Quaternion...) = q*compose_rotation(qs...)

# This algorithm was proposed by @Per in
#
#   https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/24


text:
# Autogenerated wrapper script for protoc_jll for i686-linux-musl-cxx11
export libprotoc, protoc

JLLWrappers.@generate_wrapper_header("protoc")
JLLWrappers.@declare_library_product(libprotoc, "libprotoc.so.30")
JLLWrappers.@declare_executable_product(protoc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libprotoc,
        "lib/libprotoc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        protoc,
        "bin/protoc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
# Method p51

Plane or axisymmetric strain analysis of an elastic solid (Plane
structural element) using 3-, 6-, 10- or 15-node right-angled triangles
(Triangle finite elements) or 4-, 8- or 9-node rectangular quadrilaterals
(Quadrilateral finite elements). Mesh numbered in x(r)- or y(z)- direction.

### Constructors
```julia
p51(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}` : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Structural element
* support::Array{Tuple{Int,Array{Int,1}},1}            : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1}   : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Floalt64}                       : x-coordinate vector
* y_coords::FloatRange{Floalt64}                       : y-coordinate vector
* thickness:: Float64                                  : Thickness of plate
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20             : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}         : Element material vector if np_types > 1
```

### Return values
```julia
* (fem, fm_df, sigma_df)     : Tuple of jFem, dis_df and fm_df
                               where:
                                 fm_df         : Forces and moments data table
                                 sigma_df      : Stresses data table
```

### Related help
```julia
?StructuralElement           : List of available structural element types
?Plane                       : Help on a Plane structural element
?FiniteElement               : List finite element types
?Quadrilateral               : Help on Quadrilateral finite element
```
"""
function p51(data::Dict{Symbol, Any})

  # Setup basic dimensions of arrays

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Add radial stress
  struc_el.axisymmetric && (nst = 4)

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = fin_el.nodof           # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  g_num = zeros(Int, fin_el.nod, nels)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  formnf!(nodof, nn, nf)
  neq = maximum(nf)

  # Find global array sizes

  for iel in 1:nels
    geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, struc_el.direction)
    num_to_g!(num, nf, g)
    g_num[:, iel] = num
    g_coord[:, num] = coord'
    g_g[:, iel] = g
  end

  println("There are $(neq) equations.")

  sample!(fin_el, points, weights)

  gsm = spzeros(neq, neq)
  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'              # Transpose
    g = g_g[:, iel]
    km = zeros(ndof, ndof)
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      jac = der*coord
      detm = det(jac)
      jac = inv(jac)
      deriv = jac*der
      beemat!(bee, deriv)
      if struc_el.axisymmetric
        gc = fun'*coord
        bee[4, 1:2:(ndof-1)] = fun[:]/gc[1]
      end
      km += (bee')*dee*bee*detm*weights[i]*gc[1]
    end
    PtFEM.fsparm!(gsm, g, km)
  end

  loads = OffsetArray(zeros(neq+1), 0:neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      no[i] = nf[data[:fixed_freedoms][i][2], data[:fixed_freedoms][i][1]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i],no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] * value[i]
    end
  end

  cfgsm = cholesky(gsm)
  loads[1:neq] = cfgsm \ loads[1:neq]

  displacements = zeros(size(nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]]
      end
    end
  end
  displacements = displacements'

  local dis_df
  local fm_df
  if struc_el.axisymmetric
    dis_df = DataFrame(
      x_disp = displacements[:, 1],
      z_disp = displacements[:, 2]
    )
  else
    dis_df = DataFrame(
      x_disp = displacements[:, 1],
      y_disp = displacements[:, 2]
    )
  end

  loads[0] = 0.0
  struc_el.nip = 1
  points = zeros(struc_el.nip, ndim)
  weights = zeros(struc_el.nip)
  sample!(fin_el, points, weights)

  gc1 = Vector{Float64}()
  gc2 = Vector{Float64}()
  s1 = Vector{Float64}()
  s2 = Vector{Float64}()
  s3 = Vector{Float64}()
  struc_el.axisymmetric && (s4 = Vector{Float64}())

  for iel in 1:nels
    deemat!(dee, prop[etype[iel], 1], prop[etype[iel], 2])
    num = g_num[:, iel]
    coord = g_coord[:, num]'
    g = g_g[:, iel]
    eld = loads[g]
    for i in 1:struc_el.nip
      shape_fun!(fun, points, i)
      shape_der!(der, points, i)
      gc = fun'*coord
      jac = inv(der*coord)
      deriv = jac*der
      beemat!(bee, deriv)
      if struc_el.axisymmetric
        gc = fun'*coord
        bee[4, 1:2:(ndof-1)] = fun[:]/gc[1]
      end
      sigma = dee*(bee*eld)
      append!(gc1, gc[1])
      append!(gc2, gc[2])
      append!(s1, sigma[1])
      append!(s2, sigma[2])
      append!(s3, sigma[3])
      struc_el.axisymmetric && append!(s4, sigma[4])
    end
  end
  if struc_el.axisymmetric
    fm_df = DataFrame(
      r_coord = gc1,
      z_coord = gc2,
      sig_r = s1,
      sig_z = s2,
      tau_rz = s3,
      sig_t = s4
    )
  else
    fm_df = DataFrame(
      x_coord = gc1,
      y_coord = gc2,
      sig_x = s1,
      sig_y = s2,
      tau_xy = s3
    )
  end

  fem = PtFEM.jFEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof,
    neq, penalty, etype, g, g_g, g_num, nf, no,
    node, num, sense, actions, bee, coord, gamma, dee,
    der, deriv, displacements, eld, fun, gc, g_coord, jac,
    km, mm, kg, cfgsm, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

  (fem, dis_df, fm_df)
end

max_stars_repo_path: src/5 Elastic Solids/p51.jl

max_issues_repo_path: src/5 Elastic Solids/p51.jl

max_forks_repo_path: src/5 Elastic Solids/p51.jl


text:
using Test
using Printf
using Random
using Statistics
using LinearAlgebra
using Logging

using CUDA
using JLD2
using FFTW
using OffsetArrays
using SeawaterPolynomials

using Oceananigans
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.Advection
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.Coriolis
using Oceananigans.Buoyancy
using Oceananigans.Forcings
using Oceananigans.Solvers
using Oceananigans.Models
using Oceananigans.Simulations
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TurbulenceClosures
using Oceananigans.AbstractOperations
using Oceananigans.Logger
using Oceananigans.Utils
using Oceananigans.Architectures: device # to resolve conflict with CUDA.device

using Dates: DateTime, Nanosecond
using TimesDates: TimeDate
using Statistics: mean
using LinearAlgebra: norm
using NCDatasets: Dataset
using KernelAbstractions: @kernel, @index, Event

import Oceananigans.Fields: interior
import Oceananigans.Utils: launch!, datatuple

Logging.global_logger(OceananigansLogger())

#####
##### Testing parameters
#####

float_types = (Float32, Float64)

         archs = (CPU(),)
@hascuda archs = (GPU(),)

closures = (
    :IsotropicDiffusivity,
    :AnisotropicDiffusivity,
    :AnisotropicBiharmonicDiffusivity,
    :TwoDimensionalLeith,
    :SmagorinskyLilly,
    :BlasiusSmagorinsky,
    :RozemaAnisotropicMinimumDissipation,
    :VerstappenAnisotropicMinimumDissipation
)

#####
##### Run tests!
#####

CUDA.allowscalar(true)

include("utils_for_runtests.jl")

group = get(ENV, "TEST_GROUP", :all) |> Symbol

@testset "Oceananigans" begin
    if group == :unit || group == :all
        @testset "Unit tests" begin
            include("test_grids.jl")
            include("test_operators.jl")
            include("test_boundary_conditions.jl")
            include("test_fields.jl")
            include("test_averaged_field.jl")
            include("test_kernel_computed_field.jl")
            include("test_halo_regions.jl")
            include("test_solvers.jl")
            include("test_pressure_solvers.jl")
            include("test_coriolis.jl")
            include("test_buoyancy.jl")
            include("test_surface_waves.jl")
            include("test_utils.jl")
        end
    end

    if group == :time_stepping_1 || group == :all
        @testset "Model and time stepping tests (part 1)" begin
            include("test_incompressible_models.jl")
            include("test_time_stepping.jl")
        end
    end

    if group == :time_stepping_2 || group == :all
        @testset "Model and time stepping tests (part 2)" begin
            include("test_boundary_conditions_integration.jl")
            include("test_forcings.jl")
            include("test_turbulence_closures.jl")
            include("test_dynamics.jl")
        end
    end

    if group == :simulation || group == :all
        @testset "Simulation tests" begin
            include("test_simulations.jl")
            include("test_diagnostics.jl")
            include("test_output_writers.jl")
            include("test_abstract_operations_computed_field.jl")
        end
    end

    if group == :regression || group == :all
        include("test_regression.jl")
    end

    if group == :scripts || group == :all
        @testset "Scripts" begin
            include("test_validation.jl")
        end
    end

    if group == :convergence
        include("test_convergence.jl")
    end

    if group == :shallow_water || group == :all
        include("test_shallow_water_models.jl")
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_jobs_data_plane
using AWS.Compat
using AWS.UUIDs

"""
    DescribeJobExecution()

Gets details of a job execution.

# Required Parameters
- `jobId`: The unique identifier assigned to this job when it was created.
- `thingName`: The thing name associated with the device the job execution is running on.

# Optional Parameters
- `executionNumber`: Optional. A number that identifies a particular job execution on a particular device. If not specified, the latest job execution is returned.
- `includeJobDocument`: Optional. When set to true, the response contains the job document. The default is false.
"""
describe_job_execution(jobId, thingName; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs/$(jobId)"; aws_config=aws_config)
describe_job_execution(jobId, thingName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs/$(jobId)", args; aws_config=aws_config)

"""
    GetPendingJobExecutions()

Gets the list of all jobs for a thing that are not in a terminal status.

# Required Parameters
- `thingName`: The name of the thing that is executing the job.

"""
get_pending_job_executions(thingName; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs"; aws_config=aws_config)
get_pending_job_executions(thingName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("GET", "/things/$(thingName)/jobs", args; aws_config=aws_config)

"""
    StartNextPendingJobExecution()

Gets and starts the next pending (status IN_PROGRESS or QUEUED) job execution for a thing.

# Required Parameters
- `thingName`: The name of the thing associated with the device.

# Optional Parameters
- `statusDetails`: A collection of name/value pairs that describe the status of the job execution. If not specified, the statusDetails are unchanged.
- `stepTimeoutInMinutes`: Specifies the amount of time this device has to finish execution of this job. If the job execution status is not set to a terminal state before this timer expires, or before the timer is reset (by calling UpdateJobExecution, setting the status to IN_PROGRESS and specifying a new timeout value in field stepTimeoutInMinutes) the job execution status will be automatically set to TIMED_OUT. Note that setting this timeout has no effect on that job execution timeout which may have been specified when the job was created (CreateJob using field timeoutConfig).
"""
start_next_pending_job_execution(thingName; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("PUT", "/things/$(thingName)/jobs/$next"; aws_config=aws_config)
start_next_pending_job_execution(thingName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("PUT", "/things/$(thingName)/jobs/$next", args; aws_config=aws_config)

"""
    UpdateJobExecution()

Updates the status of a job execution.

# Required Parameters
- `jobId`: The unique identifier assigned to this job when it was created.
- `status`: The new status for the job execution (IN_PROGRESS, FAILED, SUCCESS, or REJECTED). This must be specified on every update.
- `thingName`: The name of the thing associated with the device.

# Optional Parameters
- `executionNumber`: Optional. A number that identifies a particular job execution on a particular device.
- `expectedVersion`: Optional. The expected current version of the job execution. Each time you update the job execution, its version is incremented. If the version of the job execution stored in Jobs does not match, the update is rejected with a VersionMismatch error, and an ErrorResponse that contains the current job execution status data is returned. (This makes it unnecessary to perform a separate DescribeJobExecution request in order to obtain the job execution status data.)
- `includeJobDocument`: Optional. When set to true, the response contains the job document. The default is false.
- `includeJobExecutionState`: Optional. When included and set to true, the response contains the JobExecutionState data. The default is false.
- `statusDetails`:  Optional. A collection of name/value pairs that describe the status of the job execution. If not specified, the statusDetails are unchanged.
- `stepTimeoutInMinutes`: Specifies the amount of time this device has to finish execution of this job. If the job execution status is not set to a terminal state before this timer expires, or before the timer is reset (by again calling UpdateJobExecution, setting the status to IN_PROGRESS and specifying a new timeout value in this field) the job execution status will be automatically set to TIMED_OUT. Note that setting or resetting this timeout has no effect on that job execution timeout which may have been specified when the job was created (CreateJob using field timeoutConfig).
"""
update_job_execution(jobId, status, thingName; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("POST", "/things/$(thingName)/jobs/$(jobId)", Dict{String, Any}("status"=>status); aws_config=aws_config)
update_job_execution(jobId, status, thingName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_jobs_data_plane("POST", "/things/$(thingName)/jobs/$(jobId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("status"=>status), args)); aws_config=aws_config)


text:
using DSGE, Test, CSV, DataFrames
using ModelConstructors, Nullables, Dates, OrderedCollections

prob_vecs = [[1., 0.], [.5, .5], [0., 1.]]
Hbar_vec = [14, 10, 6]

m = Model1002("ss59"; custom_settings = Dict{Symbol, Setting}(:flexible_ait_policy_change =>
                                                              Setting(:flexible_ait_policy_change,
                                                                      false),
                                                              :add_pgap => Setting(:add_pgap, false),
                                                              :add_ygap => Setting(:add_ygap, false),
                                                              :add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true),
                                                              :add_altpolicy_ygap => Setting(:add_altpolicy_ygap, true))) # Set to false unless you want to re-generate any saved output
m <= Setting(:flexible_ait_policy_change, false) # Set to false unless you want to re-generate any saved output
m <= Setting(:date_forecast_start, Date(2020, 6, 30))
m <= Setting(:date_conditional_end, Date(2020, 6, 30))
m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, # Have to add anticipated rates to conditional data
                                   :obs_nominalrate, :obs_longrate,
                                   :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                   :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])
m <= Setting(:cond_semi_names, [:obs_spread,
                                   :obs_nominalrate, :obs_longrate,
                                   :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                   :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])

output_vars = [:histobs, :forecastobs, :histpseudo, :forecastpseudo]

df_full = CSV.read(joinpath(dirname(@__FILE__), "../reference/uncertain_altpolicy_data.csv"), DataFrame)
modal_params = map(x -> x.value, m.parameters)

hist_rule = DSGE.taylor_rule()
setup_regime_switching_inds!(m; cond_type = :full) # Reset the regime dates b/c using full conditional forecast

# Set up parameters
for j in 5:30 # Just add a bunch
    ModelConstructors.set_regime_val!(m[:σ_g], j, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_b], j, m[:σ_b].value)
    ModelConstructors.set_regime_val!(m[:σ_μ], j, m[:σ_μ].value)
    ModelConstructors.set_regime_val!(m[:σ_ztil], j, m[:σ_ztil].value)
    ModelConstructors.set_regime_val!(m[:σ_λ_f], j, m[:σ_λ_f].value)
    ModelConstructors.set_regime_val!(m[:σ_λ_w], j, m[:σ_λ_w].value)
    ModelConstructors.set_regime_val!(m[:σ_r_m], j, m[:σ_r_m].value)
    ModelConstructors.set_regime_val!(m[:σ_σ_ω], j, m[:σ_σ_ω].value)
    ModelConstructors.set_regime_val!(m[:σ_μ_e], j, m[:σ_μ_e].value)
    ModelConstructors.set_regime_val!(m[:σ_γ], j, m[:σ_γ].value)
    ModelConstructors.set_regime_val!(m[:σ_π_star], j, m[:σ_π_star].value)
    ModelConstructors.set_regime_val!(m[:σ_lr], j, m[:σ_lr].value)
    ModelConstructors.set_regime_val!(m[:σ_z_p], j, m[:σ_z_p].value)
    ModelConstructors.set_regime_val!(m[:σ_tfp], j, m[:σ_tfp].value)
    ModelConstructors.set_regime_val!(m[:σ_gdpdef], j, m[:σ_gdpdef].value)
    ModelConstructors.set_regime_val!(m[:σ_corepce], j, m[:σ_corepce].value)
    ModelConstructors.set_regime_val!(m[:σ_gdp], j, m[:σ_gdp].value)
    ModelConstructors.set_regime_val!(m[:σ_gdi], j, m[:σ_gdi].value)

    for l = 1:DSGE.n_mon_anticipated_shocks(m)
        ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(l)")], j, m[Symbol("σ_r_m$(l)")].value)
    end

    ModelConstructors.set_regime_val!(m[:σ_φ], j, 0.)
    ModelConstructors.set_regime_val!(m[:σ_biidc], j, 0.)
    ModelConstructors.set_regime_val!(m[:σ_ziid], j, 0.)
    if haskey(m.keys, :σ_φ_prop)
        ModelConstructors.set_regime_val!(m[:σ_φ_prop], j, 0.)
    end
    if haskey(m.keys, :σ_biidc_prop)
        ModelConstructors.set_regime_val!(m[:σ_biidc_prop], j, 0.)
    end
    if haskey(m.keys, :σ_ziid_prop)
        ModelConstructors.set_regime_val!(m[:σ_ziid_prop], j, 0.)
    end
end

# Impose the historical policy as the permanent and "alternative" rule
T_unc = Dict()
R_unc = Dict()
C_unc = Dict()
for j in 1:length(prob_vecs)
    T_unc[j], R_unc[j], C_unc[j] = DSGE.gensys_uncertain_altpol(m, prob_vecs[j], DSGE.AltPolicy[hist_rule])
    T_unc[j], R_unc[j], C_unc[j] = DSGE.augment_states(m, T_unc[j], R_unc[j], C_unc[j])
end

# Some setup
tempZLB_regimes = Dict()
m <= Setting(:alternative_policies, AltPolicy[hist_rule])

autoTs = Dict() # automatically calculated, these are what we are testing
autoRs = Dict()
autoCs = Dict()
auto_fcasts = Dict()

manTs = Dict()  # manually calculated, these are the "truth"
manRs = Dict()
manCs = Dict()
fcasts = Dict()

# Calculate baseline forecast to calculate the states in 2020-Q3
baseline_fcasts = DSGE.forecast_one_draw(m, :mode, :full, vcat(:forecaststates, output_vars), modal_params,
                                         df_full; regime_switching = true, n_regimes =
                                         get_setting(m, :n_regimes))

# Now calculate everything else
for j in 1:length(prob_vecs)
    # Set up initial state
    s₀ = baseline_fcasts[:forecaststates][:, 1]

    # Set up temporary ZLB regimes
    tempZLB_regimes[j] = 4:(4 + Hbar_vec[j] - 1)
    new_regime_dates = Dict{Int, Date}()
    new_regime_dates[1] = date_presample_start(m)
    for (regind, date) in zip(2:(maximum(tempZLB_regimes[j]) + 1), # This won't actually go to 3000, but this makes sure we have
                              Date(2020, 3, 31):Dates.Month(3):Date(3000,3,31)) # enough regimes to fill up hist_regime_dates properly
        new_regime_dates[regind] = date
    end

    m <= Setting(:regime_dates, new_regime_dates)
    setup_regime_switching_inds!(m; cond_type = :full)

    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys2_separate_cond_regimes, true) # required b/c not using an alternative policy technically
    replace_eqcond = Dict{Int, DSGE.EqcondEntry}()     # Which rule to replace with in which periods
    for regind in tempZLB_regimes[j]
        replace_eqcond[regind] = DSGE.EqcondEntry(DSGE.zlb_rule(), prob_vecs[j]) # Temp ZLB rule in this regimes
    end
    replace_eqcond[tempZLB_regimes[j][end]+1] = DSGE.EqcondEntry(AltPolicy(:historical, eqcond, solve), prob_vecs[j]) # Will set the "lift-off" regime directly
    m <= Setting(:regime_eqcond_info, replace_eqcond)
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])
    fcasts[j] = deepcopy(baseline_fcasts)

    # Set up matrices for regime-switching. Note that we have
    # plus 1 for lift-off b/c historical and conditional forecast period are accounted for by fact
    # tempZLB_regimes calculated from forecastobs (which includes conditional forecast periods)
    Γ0s = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)
    Γ1s = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)
    Cs = Vector{Vector{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)
    Ψs = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)
    Πs = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)

    TTTs_uzlb = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1) # uncertain ZLB
    RRRs_uzlb = Vector{Matrix{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)
    CCCs_uzlb = Vector{Vector{Float64}}(undef, maximum(tempZLB_regimes[j]) + 1)

    # Automatic calculation
    m <= Setting(:uncertain_altpolicy, true)
    m <= Setting(:uncertain_temporary_altpolicy, true)
    autoTs[j], autoRs[j], autoCs[j] = solve(m; regime_switching = true,
                                            gensys_regimes = [1:get_setting(m, :n_hist_regimes) + 1],
                                            gensys2_regimes =
                                            [(get_setting(m, :n_hist_regimes) + 1):get_setting(m, :n_regimes)],
                                            regimes = collect(1:get_setting(m, :n_regimes)))
    auto_fcasts[j] = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params,
                                            df_full; regime_switching = true, n_regimes =
                                            get_setting(m, :n_regimes))
    m <= Setting(:uncertain_altpolicy, false) # Need to turn these off afterward
    m <= Setting(:uncertain_temporary_altpolicy, false)

    # Onto the manual calculation

    # Only need to fill up the matrices for the forecast regimes for this part
    for fcast_reg in 3:length(TTTs_uzlb)
        Γ0s[fcast_reg], Γ1s[fcast_reg], Cs[fcast_reg], Ψs[fcast_reg], Πs[fcast_reg] =
        DSGE.zlb_rule().eqcond(m, fcast_reg)
    end
    TTT_gensys_final, RRR_gensys_final, CCC_gensys_final = solve(m)
    TTT_gensys_final = TTT_gensys_final[1:n_states(m), 1:n_states(m)]
    RRR_gensys_final = RRR_gensys_final[1:n_states(m), :]
    CCC_gensys_final = CCC_gensys_final[1:n_states(m)]

    # Calculate temp ZLB regimes under assumption smooth AIT-GDP occurs forever after ZLB ends
    gensys2_regimes = (tempZLB_regimes[j][1] - 1):(tempZLB_regimes[j][end] + 1)
    Tcal, Rcal, Ccal = DSGE.gensys2(m, Γ0s[gensys2_regimes], Γ1s[gensys2_regimes],
                                    Cs[gensys2_regimes], Ψs[gensys2_regimes], Πs[gensys2_regimes],
                                    TTT_gensys_final, RRR_gensys_final, CCC_gensys_final, Hbar_vec[j] + 1)
    Tcal[end] = TTT_gensys_final
    Rcal[end] = RRR_gensys_final
    Ccal[end] = CCC_gensys_final

    Th, Rh, Ch = hist_rule.solve(m)
    TTTs_uzlb[1] = Th
    TTTs_uzlb[2] = Th
    TTTs_uzlb[3] = Th
    CCCs_uzlb[1] = Ch
    CCCs_uzlb[2] = Ch
    CCCs_uzlb[3] = Ch

    # Prep for calculation of new transition matrices
    Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = DSGE.gensys_to_predictable_form(Γ0s[4], Γ1s[4], Cs[4], Ψs[4], Πs[4])

    # Get new transition matrices for "uncertain" ZLB
    # Note we pass in Tcal[2:end] since we want the period t + 1 ZLB matrix for period t.
    # TTTs_uzlb[inreg] is the same size as Tcal though b/c TTTs_uzbl[inreg[end]] is just Tcal[end], etc.
    inreg = gensys2_regimes[2:end]
    TTTs_uzlb[inreg], RRRs_uzlb[inreg], CCCs_uzlb[inreg] =
        DSGE.gensys2_uncertain_altpol(prob_vecs[j], Th[1:n_states(m), 1:n_states(m)],
                                      Ch[1:n_states(m)], Tcal[2:end], Rcal[2:end], Ccal[2:end],
                                      Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til)
    for li in inreg
        TTTs_uzlb[li], RRRs_uzlb[li], CCCs_uzlb[li] = DSGE.augment_states(m, TTTs_uzlb[li], RRRs_uzlb[li], CCCs_uzlb[li])
    end

    # Now update the last regime matrices to feature "uncertainty" about the alt policy, too, post ZLB
    TTTs_uzlb[end] = T_unc[j]
    RRRs_uzlb[end] = R_unc[j]
    CCCs_uzlb[end] = C_unc[j]

    manTs[j] = TTTs_uzlb
    manRs[j] = RRRs_uzlb
    manCs[j] = CCCs_uzlb
end

@testset "Check automatic calculation of uncertain ZLB with an uncertain alternative policy." begin
    obs_inds = vcat(1:9, 12:13) # Only consider non-anticipating variables b/c still handling how to calculate measurement equations for them
    for j in 1:length(prob_vecs)
        for i in 4:length(manTs[j])
            @test @test_matrix_approx_eq autoTs[j][i] manTs[j][i]
        end
    end
end

nothing