text:
using CompScienceMeshes
using SauterSchwabQuadrature
#create a subdivision surfaces for a sphere exaple
G = readmesh("/Users/Benjamin/Documents/sphere.in")
subdG=GSubdMesh(G)

chart = chart(subdG,1)
u1= [0.5,0.25]
u2= [0.7,0.1]
d1 = neighborhood(chart,u1)
d2 = neighborhood(chart,u2)
global p1 = cartesian(d1)
global p2 = cartesian(d2)
function kernel(x,y)
			return(((x-p1)'*(y-p2))*exp(-im*1*norm(x-y))/(4pi*norm(x-y)))	#k=1
end
result = sauterschwabintegral(chart, kernel)


text:
"""
$(SIGNATURES)

Solve Schrodinger equation defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total annealing time.
- `tspan` = (0, tf): time interval to solve the dynamics.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_schrodinger(A::Annealing, tf::Real; tspan = (0, tf), kwargs...)
    u0 = build_u0(A.u0, :v)
    p = ODEParams(A.H, float(tf), A.annealing_parameter)
    update_func = function (C, u, p, t)
        update_cache!(C, p.L, p, p(t))
    end
    cache = get_cache(A.H)
    diff_op = DiffEqArrayOperator(cache, update_func = update_func)
    jac_cache = similar(cache)
    jac_op = DiffEqArrayOperator(jac_cache, update_func = update_func)
    ff = ODEFunction(diff_op, jac_prototype = jac_op)

    prob = ODEProblem{true}(ff, u0, float.(tspan), p)
    alg_keyword_warning(;kwargs...)
    solve(prob; kwargs...)
end

"""
$(SIGNATURES)

Solve the unitary defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total annealing time.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `tspan` = (0, tf): time interval to solve the dynamics.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_unitary(
    A::Annealing,
    tf::Real;
    vectorize::Bool = false,
    tspan = (0, tf),
    kwargs...,
)

    H = A.H
    u0 = build_u0(Matrix{ComplexF64}(I, size(H)), :m, vectorize = vectorize)
    p = ODEParams(A.H, float(tf), A.annealing_parameter)

    uni_jac = function (J, u, p, t)
        hmat = p.L(t)
        J .= -1.0im * one(hmat) ⊗ hmat
    end

    cache = get_cache(H)
    if vectorize == false
        diff_op_update = (C, u, p, t) -> update_cache!(C, p.L, p, p(t))
    else
        cache = vectorize_cache(cache)
        diff_op_update = uni_jac
    end
    j_cache = similar(cache)
    diff_op = DiffEqArrayOperator(cache, update_func = diff_op_update)
    jac_op = DiffEqArrayOperator(j_cache, update_func = uni_jac)
    ff = ODEFunction(diff_op, jac_prototype = jac_op)

    prob = ODEProblem{true}(ff, u0, float.(tspan), p)
    alg_keyword_warning(;kwargs...)
    solve(prob; kwargs...)
end

"""
$(SIGNATURES)

Solve the von Neumann equation defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total annealing time.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `tspan` = (0, tf): time interval to solve the dynamics.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_von_neumann(
    A::Annealing,
    tf::Real;
    tspan = (0, tf),
    vectorize::Bool = false,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize = vectorize)
    von_f = function (du, u, p, t)
        p.L(du, u, p, p(t))
    end

    von_jac = function (J, u, p, t)
        hmat = p.L(p(t))
        iden = one(hmat)
        J .= -1.0im * (iden ⊗ hmat - transpose(hmat) ⊗ iden)
    end

    if vectorize == false
        ff = ODEFunction{true}(von_f, jac = von_jac)
    else
        cache = vectorize_cache(get_cache(A.H))
        diff_op = DiffEqArrayOperator(
            cache,
            update_func = (A, u, p, t) ->
                update_vectorized_cache!(A, p.L, p, p(t)),
        )
        jac_cache = similar(cache)
        jac_op = DiffEqArrayOperator(
            jac_cache,
            update_func = (A, u, p, t) ->
                update_vectorized_cache!(A, p.L, p, p(t)),
        )
        ff = ODEFunction(diff_op, jac_prototype = jac_op)
    end

    p = ODEParams(A.H, float(tf), A.annealing_parameter)
    prob = ODEProblem{true}(ff, u0, float.(tspan), p)
    solve(prob; kwargs...)
end

ext: jl


text:
y=Array(Array{Int,1},2)
y[1]=[5,7,3]
y[2]=[10]

sa=Array(Array{Int,2},2)
sa[1]=[5 1; 2 2; 3 1]
sa[2]=[5 2]

dims=(5,2)

@test RustModels.sa2y(dims,sa[1])==y[1]
@test RustModels.sa2y(dims,sa)==y
@test RustModels.y2sa(dims,y[1])==sa[1]
@test RustModels.y2sa(dims,y)==sa


text:
module Tokenizer

export tokenize,
       to_lines,
       to_sentences,
       to_letters,
       cleanup,
       words,
       letters,
       lines

"""
    tokenize(text[, on=letters])

Split `text` into SupTokens (array of arrays of tokens). An optional function
of general type `func(::Any) -> Vector{Vector{Any}}` can be provided to be used
for the tokenization.

For possible *combinators* which can be composed to obtain `func`, see:
[`to_lines`](@ref), [`to_sentences`](@ref), [`to_letters`](@ref), [`to_words`](@ref),
[`cleanup`](@ref).
"""
function tokenize(text; on=letters)
    return on(text)
end

"""
    to_lines(text::AbstractString)

Return an array of lines in `text`.
"""
function to_lines(text::AbstractString)
    return split(text, "\n")
end

"""
    to_sentences(text::AbstractString)

Return an array of sentences in `text`. The text is split along dots; the dots
remain in the strings, only the spaces after the dots are stripped.

The function tries to be as smart as possible. For example, the string
`"Channel No. 5 is a perfume."` will be treated as one sentence,
although it has two dots.
"""
function to_sentences(text::AbstractString)
    # Split on (optional) whitespace if it's preceeded by a dot
    # and if it's followed by a capital letter
    rule = r"((?<=[.])\s*(?=[A-Z]))|((?<=[?!])\s*)"
    split(text, rule; keepempty=false)
end

"""
    to_letters(tokens::Vector{<:AbstractString})

Split all of the tokens in `tokens` into individual characters.
"""
function to_letters(tokens::Vector{<:AbstractString})
    return [split(token, "") for token in tokens]
end

"""
    to_words(tokens::Vector{<:AbstractString}; keeppunctuation=true)

Split all of the tokens in `tokens` into individual words by whitespace.
If `keeppunctuation` is true, all of the special characters are preserved
(and thus "glued" to the preceding/following word).
"""
function to_words(tokens::Vector{<:AbstractString}; keeppunctuation=true)
    rule = if keeppunctuation r"\s+" else r"\W+" end
    return [split(token, rule; keepempty=false) for token in tokens]
end

"""
    cleanup(suptokens::Vector{<:Vector{<:AbstractString}}; badchars="»«\\n-_()[]{}<>–—\$=\'\"„“\r\t")

Remove all characters that are in `badchars` from all tokens in `suptokens`.
"""
function cleanup(suptokens::Vector{<:Vector{<:AbstractString}}; badchars="»«\n-_()[]{}<>–—\$=\'\"„“\r\t")
    cleanup_token(token) = filter(c -> !(c in badchars), token)

    return [
        # A list of non-empty cleaned-up tokens
        [cleanup_token(token) for token in suptoken if cleanup_token(token) != ""]
        # The list is built for every item in suptokens
        for suptoken in suptokens
    ]
end

"""
    letters = cleanup ∘ to_letters ∘ to_sentences

Composite function which splits its input into sentences, then the sentences
into letters, and then removes special characters.
"""
letters = cleanup ∘ to_letters ∘ to_sentences

"""
    lines = cleanup ∘ to_letters ∘ to_sentences

Composite function which splits its input into lines, then the line
into letters, and then removes special characters.
"""
lines = cleanup ∘ to_letters ∘ to_lines

"""
    words = cleanup ∘ to_letters ∘ to_sentences

Composite function which splits its input into sentences, then the sentences
into words, and then removes special characters. Please note that dots and
commas are not removed.
"""
words = cleanup ∘ to_words ∘ to_sentences

end


text:
module scratch1_06092017_ortho
using FinEtools
using FinEtoolsDeforLinear
using Test

mutable struct MyIData
  c::FInt
  r::FFltVec
  s::FFltVec
end

function test()

  # println("Thick pipe with internal pressure: axially symmetric model")
  #=
  This is a simple modification of the full three-dimensional simulation of
  the tutorial pub_thick_pipe that implements the axially-symmetric model
  reduction procedure.

  An infinitely long thick walled cylindrical pipe
  with inner boundary radius of 3 mm and outer boundary radius of 9 mm is
  subjected to an internal pressure of 1.0 MPa. A wedge   with thickness of
  2 mm and a 90-degree angle sector is considered for the finite element
  analysis. The material properties are taken as  isotropic linear elastic
  with $E=1000$ MPa and $\nu=0.4999$ to represent nearly incompressible
  behavior. This problem has been proposed to by MacNeal and Harder as a
  test of an element's ability to represent the  response of a nearly
  incompressible material. The plane-strain condition is assumed in the
  axial direction of the pipe which together with the radial symmetry
  confines the material in all but the radial direction and therefore
  amplifies the numerical difficulties associated with the confinement of
  the nearly incompressible material.

  There is an analytical solution to this problem. Timoshenko and Goodier
  presented the original solution of Lame in their textbook. We are going
  to compare with  both the stress distribution (radial and hoop stresses)
  and the displacement of the inner  cylindrical surface.

  References:
  - Macneal RH, Harder RL (1985) A proposed standard set of problems to test
  finite element accuracy. Finite Elements in Analysis and Design 1: 3-20.
  - Timoshenko S. and Goodier J. N., Theory of Elasticity, McGraw-Hill, 2nd ed., 1951.

  =#

  # Internal radius of the pipe.
  a=3*phun("MM");
  ##
  # External radius of the pipe.
  b = 9*phun("MM");
  ##
  # Thickness of the slice.
  t = 2*phun("MM");

  ##
  # Geometrical tolerance.
  tolerance   = a/10000.;
  ##
  # Young's modulus and Poisson's ratio.
  E = 1000*phun("MEGA*PA");
  nu = 0.499;
  ##
  # Applied pressure on the internal surface.
  press =   1.0*phun("MEGA*PA");

  ##
  # Analytical solutions.   Radial stress:
  radial_stress(r)  = press*a.^2/(b^2-a^2).*(1-(b^2)./r.^2);
  ##
  # Circumferential (hoop) stress:
  hoop_stress(r) = press*a.^2/(b^2-a^2).*(1+(b^2)./r.^2);

  ##
  # Radial displacement:
  radial_displacement(r) = press*a^2*(1+nu)*(b^2+r.^2*(1-2*nu))/(E*(b^2-a^2).*r);;

  ##
  # Therefore the radial displacement of the loaded surface will be:
  urex  =  radial_displacement(a);

  ##
  # The mesh parameters: The numbers of element edges axially,
  # and through the thickness of the pipe wall (radially).

  na = 1; nt = 10;

  ##
  # Note that the material object needs to be created with the proper
  # model-dimension reduction in effect.  In this case that is the axial symmetry
  # assumption.
  MR = DeforModelRed2DAxisymm
  axisymmetric = true

  # Create the mesh and initialize the geometry.  First we are going
  # to construct the block of elements with the first coordinate
  # corresponding to the thickness in the radial direction, and the second
  # coordinate is the thickness in the axial direction.
  fens,fes  =   Q8block(b-a, t, nt, na);

  # Extract the boundary  and mark the finite elements on the
  # interior surface.
  bdryfes = meshboundary(fes);

  bcl = selectelem(fens, bdryfes, box=[0.,0.,-Inf,Inf], inflate=tolerance);
  internal_fenids= connectednodes(subset(bdryfes,bcl));
  # Now  shape the block  into  the actual wedge piece of the pipe.
  for i=1:count(fens)
    fens.xyz[i,:] = fens.xyz[i,:] + [a; 0.0];
  end

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

  # The plane-strain condition in the axial direction  is specified by selecting nodes
  # on the plane y=0 and y=t.
  l1 = selectnode(fens; box=[-Inf Inf 0.0 0.0], inflate = tolerance)
  setebc!(u, l1, true, 2, 0.0)
  l1 = selectnode(fens; box=[-Inf Inf t t], inflate = tolerance)
  setebc!(u, l1, true, 2, 0.0)

  applyebc!(u)
  numberdofs!(u)

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

  el1femm =  FEMMBase(IntegDomain(subset(bdryfes,bcl), GaussRule(1, 3), axisymmetric))
  fi = ForceIntensity([press; 0.0]);
  F2= distribloads(el1femm, geom, u, fi, 2);

  # Property and material
  material = MatDeforElastOrtho(MR,  E, nu)

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

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

  # Transfer the solution of the displacement to the nodes on the
  # internal cylindrical surface and convert to
  # cylindrical-coordinate displacements there.
  uv=u.values[internal_fenids,:]
  # Report the  relative displacement on the internal surface:
  # println("(Approximate/true displacement) at the internal surface: $( mean(uv[:,1])/urex*100  ) %")

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

  fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)

  # File =  "thick_pipe_sigmax.vtk"
  # vtkexportmesh(File, fens, fes; scalars=[("sigmax", fld.values)])

  # Produce a plot of the solution components in the cylindrical
  # coordinate system.

  function inspector(idat::MyIData, elnum, conn, xe,  out,  xq)
    push!(idat.r, xq[1])
    push!(idat.s, out[idat.c])
    return idat
  end

  idat = MyIData(1, FInt[], FInt[])
  idat = inspectintegpoints(femm, geom, u, collect(1:count(fes)),
  inspector, idat, :Cauchy)

  # using Plots
  # plotly()
  #
  # # Plot the analytical solution.
  # r = linearspace(a,b,100);
  # plot(r, radial_stress(r))
  # # Plot the computed  integration-point data
  # plot!(idat.r, idat.s, m=:circle, color=:red)
  # gui()

  @test  abs(idat.r[1]-0.003126794919243112)<1.0e-9
  @test  abs(idat.s[1]- -910911.9777008593)<1.0e-2

  ## Discussion
  #
  ##
  # The axially symmetric model is clearly very effective
  # computationally, as the size is much reduced compared to the 3-D
  # model.  In conjunction with  uniform or selective reduced integration
  # it can be very accurate as well.
end
end
using .scratch1_06092017_ortho
scratch1_06092017_ortho.test()

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

  # NAFEMS LE11 benchmark with Q8 elements.
  # # This is a test recommended by the National Agency for Finite Element
  # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
  # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
  # #
  # # Target solution: Direct stress,   =  –105 MPa at point A.
  #function  LE11NAFEMS()
  # Parameters:
  Ea =  210000*phun("MEGA*Pa")
  nua =  0.3;
  alphaa = 2.3e-4;              # thermal expansion coefficient
  sigmaA = -105*phun("MEGA*Pa")
  nref =  1;                        # how many times should we refine the mesh?
  X = [1.     0.;#A
  1.4    0.;#B
  0.995184726672197   0.098017140329561;
  1.393258617341076 0.137223996461385;
  0.980785  0.195090;#
  1.37309939 0.27312645;
  0.956940335732209   0.290284677254462
  1.339716470025092 0.406398548156247
  0.9238795  0.38268;#C
  1.2124  0.7;#D
  0.7071  0.7071;#E
  1.1062  1.045;#F
  0.7071  (0.7071+1.79)/2;#(E+H)/2
  1.      1.39;#G
  0.7071  1.79;#H
  1.      1.79;#I
  ]*phun("M")
  tolerance  = 1.e-6*phun("M")
  ##
  # Note that the material object needs to be created with the proper
  # model-dimension reduction in mind.  In this case that is the axial symmetry
  # assumption.
  MR  =  DeforModelRed2DAxisymm

  fens = FENodeSet(X);
  fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
  for ref = 1:nref
    fens,fes = Q4refine(fens,fes);
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
  end
  fens,fes = Q4toQ8(fens,fes)
  list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
  fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

  #     File  =   "mesh.vtk"
  # vtkexportmesh(File, fens, fes)

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

  # Apply EBC's
  l1 = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
  setebc!(u, l1, true, 2, 00.0)
  l1 = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
  setebc!(u, l1, true, 2, 00.0)
  applyebc!(u)
  numberdofs!(u)

  # Temperature field
  dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

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

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

  K  = stiffness(femm, geom, u)
  F  =  thermalstrainloads(femm, geom, u, dT)
  #K = cholesky(K)
  U =   K\F
  scattersysvec!(u, U[:])

  nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

  File  =   "LE11NAFEMS_Q8_sigmay.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.443052182185006e8, -1.4106181545272605e7]) < 1.0e-2
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  sA  =  fld.values[nA]/phun("MEGA*Pa")
  sAn  =  fld.values[nA]/sigmaA
  # println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")
  @test abs(sA[1] - (-93.8569)) < 1.0e-3

  fen2fe  = FENodeToFEMap(connasarray(fes), nnodes(geom))
  function inspector(idat, elnum, conn, xe,  out,  xq)
    println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
    return idat
  end

  # inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
  # inspector, []; quantity = :Cauchy)

end
end
using .mmLE11Q8mm
mmLE11Q8mm.test()

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

  # NAFEMS LE11 benchmark with Q8 elements.
  # # This is a test recommended by the National Agency for Finite Element
  # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
  # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
  # #
  # # Target solution: Direct stress,   =  –105 MPa at point A.
  #function  LE11NAFEMS()
  # Parameters:
  Ea =  210000*phun("MEGA*Pa")
  nua =  0.3;
  alphaa = 2.3e-4;              # thermal expansion coefficient
  sigmaA = -105*phun("MEGA*Pa")
  nref =  1;                        # how many times should we refine the mesh?
  X = [1.     0.;#A
  1.4    0.;#B
  0.995184726672197   0.098017140329561;
  1.393258617341076 0.137223996461385;
  0.980785  0.195090;#
  1.37309939 0.27312645;
  0.956940335732209   0.290284677254462
  1.339716470025092 0.406398548156247
  0.9238795  0.38268;#C
  1.2124  0.7;#D
  0.7071  0.7071;#E
  1.1062  1.045;#F
  0.7071  (0.7071+1.79)/2;#(E+H)/2
  1.      1.39;#G
  0.7071  1.79;#H
  1.      1.79;#I
  ]*phun("M")
  tolerance  = 1.e-6*phun("M")
  ##
  # Note that the material object needs to be created with the proper
  # model-dimension reduction in mind.  In this case that is the axial symmetry
  # assumption.
  MR  =  DeforModelRed2DAxisymm

  fens = FENodeSet(X);
  fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
  for ref = 1:nref
    fens,fes = Q4refine(fens,fes);
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
  end
  fens,fes = Q4toQ8(fens,fes)
  list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
  fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

  #     File  =   "mesh.vtk"
  # vtkexportmesh(File, fens, fes)

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

  # Apply EBC's
  l1 = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
  setebc!(u, l1, true, 2, 00.0)
  l1 = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
  setebc!(u, l1, true, 2, 00.0)
  applyebc!(u)
  numberdofs!(u)

  # Temperature field
  dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

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

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

  K  = stiffness(femm, geom, u)
  F  =  thermalstrainloads(femm, geom, u, dT)
  #K = cholesky(K)
  U =   K\F
  scattersysvec!(u, U[:])

  nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

  File  =   "LE11NAFEMS_Q8_sigmay.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.443052182185006e8, -1.4106181545272605e7]) < 1.0e-2
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  sA  =  fld.values[nA]/phun("MEGA*Pa")
  sAn  =  fld.values[nA]/sigmaA
  # println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")
  @test abs(sA[1] - (-93.8569)) < 1.0e-3

  fen2fe  = FENodeToFEMap(connasarray(fes), nnodes(geom))
  function inspector(idat, elnum, conn, xe,  out,  xq)
    println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
    return idat
  end

  # inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
  # inspector, []; quantity = :Cauchy)

end
end
using .mmLE11Q8mmortho
mmLE11Q8mmortho.test()

module mLE11Q8aximmm
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
  # NAFEMS LE11 benchmark with Q8 elements.
  # # This is a test recommended by the National Agency for Finite Element
  # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
  # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
  # #
  # # Target solution: Direct stress,   =  –105 MPa at point A.
  #function  LE11NAFEMS()
  # Parameters:
  Ea =  210000*phun("MEGA*Pa")
  nua =  0.3;
  alphaa = 2.3e-4;              # thermal expansion coefficient
  sigmaA = -105*phun("MEGA*Pa")
  nref = 2;                        # how many times should we refine the mesh?
  X = [1.     0.;#A
  1.4    0.;#B
  0.995184726672197   0.098017140329561;
  1.393258617341076 0.137223996461385;
  0.980785  0.195090;#
  1.37309939 0.27312645;
  0.956940335732209   0.290284677254462
  1.339716470025092 0.406398548156247
  0.9238795  0.38268;#C
  1.2124  0.7;#D
  0.7071  0.7071;#E
  1.1062  1.045;#F
  0.7071  (0.7071+1.79)/2;#(E+H)/2
  1.      1.39;#G
  0.7071  1.79;#H
  1.      1.79;#I
  ]*phun("M")
  tolerance  = 1.e-6*phun("M")
  ##
  # Note that the material object needs to be created with the proper
  # model-dimension reduction in mind.  In this case that is the axial symmetry
  # assumption.
  MR  =  DeforModelRed2DAxisymm

  fens = FENodeSet(X);
  fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
  for ref = 1:nref
    fens,fes = Q4refine(fens,fes);
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
  end
  fens,fes = Q4toQ8(fens,fes)
  list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
  fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

  #     File  =   "mesh.vtk"
  # vtkexportmesh(File, fens, fes)

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

  # Apply EBC's
  lbottom = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
  setebc!(u, lbottom, true, 2, 00.0)
  ltop = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
  setebc!(u, ltop, true, 2, 00.0)
  applyebc!(u)
  numberdofs!(u)

  # Temperature field
  dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

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

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

  K  = stiffness(femm, geom, u)
  F  =  thermalstrainloads(femm, geom, u, dT)
  #K = cholesky(K)
  U =   K\F
  scattersysvec!(u, U[:])

  nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

  File  =   "LE11NAFEMS_Q8_sigmay.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.6338426447540134e8, -4.961956343464769e6]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  sA  =  fld.values[nA]/phun("MEGA*Pa")
  sAn  =  fld.values[nA]/sigmaA
  # println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

  fen2fe  = FENodeToFEMap(connasarray(fes), nnodes(geom))
  function inspector(idat, elnum, conn, xe,  out,  xq)
    # println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
    return idat
  end

  inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
  inspector, []; quantity = :Cauchy)

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Pressure, 1)
  File  =   "LE11NAFEMS_Q8_pressure.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  pressure = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.1881819144904878e7, 7.555030948761216e7]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  fld =  fieldfromintegpoints(femm, geom, u, dT, :vm, 1)
  File  =   "LE11NAFEMS_Q8_vm.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
  vectors = [("u", u.values)])
  # println("range of von Mises = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [3.221370699152578e7, 1.4437590830351183e8]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  AE = AbaqusExporter("LE11NAFEMS_Q8_export_stress");
  HEADING(AE, "NAFEMS LE11 benchmark with Q8 elements.");
  COMMENT(AE, "sigmaA = -105 MPa ");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "cax8", "AllElements", 1, connasarray(fes))
  NSET_NSET(AE, "ltop", ltop)
  NSET_NSET(AE, "lbottom", lbottom)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, Ea, nua)
  EXPANSION(AE, alphaa)
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.ltop", 2)
  BOUNDARY(AE, "ASSEM1.INSTNC1.lbottom", 2)
  TEMPERATURE(AE, "ASSEM1.INSTNC1.", collect(1:count(fens)), vec(dT.values))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open(AE.filename) do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 963
  try rm(AE.filename); catch end

end
end
using .mLE11Q8aximmm
mLE11Q8aximmm.test()

module mLE11Q8aximorthom
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
  # NAFEMS LE11 benchmark with Q8 elements.
  # # This is a test recommended by the National Agency for Finite Element
  # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
  # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
  # #
  # # Target solution: Direct stress,   =  –105 MPa at point A.
  #function  LE11NAFEMS()
  # Parameters:
  Ea =  210000*phun("MEGA*Pa")
  nua =  0.3;
  alphaa = 2.3e-4;              # thermal expansion coefficient
  sigmaA = -105*phun("MEGA*Pa")
  nref = 2;                        # how many times should we refine the mesh?
  X = [1.     0.;#A
  1.4    0.;#B
  0.995184726672197   0.098017140329561;
  1.393258617341076 0.137223996461385;
  0.980785  0.195090;#
  1.37309939 0.27312645;
  0.956940335732209   0.290284677254462
  1.339716470025092 0.406398548156247
  0.9238795  0.38268;#C
  1.2124  0.7;#D
  0.7071  0.7071;#E
  1.1062  1.045;#F
  0.7071  (0.7071+1.79)/2;#(E+H)/2
  1.      1.39;#G
  0.7071  1.79;#H
  1.      1.79;#I
  ]*phun("M")
  tolerance  = 1.e-6*phun("M")
  ##
  # Note that the material object needs to be created with the proper
  # model-dimension reduction in mind.  In this case that is the axial symmetry
  # assumption.
  MR  =  DeforModelRed2DAxisymm

  fens = FENodeSet(X);
  fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
  for ref = 1:nref
    fens,fes = Q4refine(fens,fes);
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
  end
  fens,fes = Q4toQ8(fens,fes)
  list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
  fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

  #     File  =   "mesh.vtk"
  # vtkexportmesh(File, fens, fes)

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

  # Apply EBC's
  lbottom = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
  setebc!(u, lbottom, true, 2, 00.0)
  ltop = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
  setebc!(u, ltop, true, 2, 00.0)
  applyebc!(u)
  numberdofs!(u)

  # Temperature field
  dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

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

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

  K  = stiffness(femm, geom, u)
  F  =  thermalstrainloads(femm, geom, u, dT)
  #K = cholesky(K)
  U =   K\F
  scattersysvec!(u, U[:])

  nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

  File  =   "LE11NAFEMS_Q8_sigmay.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.6338426447540134e8, -4.961956343464769e6]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  sA  =  fld.values[nA]/phun("MEGA*Pa")
  sAn  =  fld.values[nA]/sigmaA
  # println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

  fen2fe  = FENodeToFEMap(connasarray(fes), nnodes(geom))
  function inspector(idat, elnum, conn, xe,  out,  xq)
    # println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
    return idat
  end

  inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
  inspector, []; quantity = :Cauchy)

  fld =  fieldfromintegpoints(femm, geom, u, dT, :Pressure, 1)
  File  =   "LE11NAFEMS_Q8_pressure.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
  vectors = [("u", u.values)])
  # println("range of  pressure = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [-1.1881819144904878e7, 7.555030948761216e7]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  fld =  fieldfromintegpoints(femm, geom, u, dT, :vm, 1)
  File  =   "LE11NAFEMS_Q8_vm.vtk"
  vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
  vectors = [("u", u.values)])
  # println("range of von Mises = $((minimum(fld.values), maximum(fld.values)))")
  @test norm([minimum(fld.values), maximum(fld.values)] - [3.221370699152578e7, 1.4437590830351183e8]) < 1.e-1
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  AE = AbaqusExporter("LE11NAFEMS_Q8_export_stress");
  HEADING(AE, "NAFEMS LE11 benchmark with Q8 elements.");
  COMMENT(AE, "sigmaA = -105 MPa ");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "cax8", "AllElements", 1, connasarray(fes))
  NSET_NSET(AE, "ltop", ltop)
  NSET_NSET(AE, "lbottom", lbottom)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, Ea, nua)
  EXPANSION(AE, alphaa)
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.ltop", 2)
  BOUNDARY(AE, "ASSEM1.INSTNC1.lbottom", 2)
  TEMPERATURE(AE, "ASSEM1.INSTNC1.", collect(1:count(fens)), vec(dT.values))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open(AE.filename) do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 963
  try rm(AE.filename); catch end

end
end
using .mLE11Q8aximorthom
mLE11Q8aximorthom.test()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
abstract type Forecast <: Component end

struct Deterministic <: Forecast
    device::Device
    horizon::Int
    resolution::Dates.Period
    interval::Dates.Period
    initialtime::Dates.DateTime
    data::Dict{Any,Dict{Int,TimeSeries.TimeArray}}
end

struct Scenarios <: Forecast
    horizon::Int
    resolution::Dates.Period
    interval::Dates.Period
    initialtime::Dates.DateTime
    scenarioquantity::Int
    data::Dict{Any,Dict{Int,TimeSeries.TimeArray}}
end

struct Probabilistic <: Forecast
    horizon::Int
    resolution::Dates.Period
    interval::Dates.Period
    initialtime::Dates.DateTime
    percentilequantity::Int
    data::Dict{Any,Dict{Int,TimeSeries.TimeArray}}
end


text:
"""
    abstract type TangentBundle{N, B}; end

This type represents the `N`-th order (iterated) tangent bundle [1] `TⁿB` over some base
(Riemannian) manifold `B`. Note that `TⁿB` is itself another manifold and thus
in particular a vector space (over ℝ). As such, subtypes of this abstract type
are expected to support the usual vector space operations.

However, beyond that, this abstract type makes no guarantee about the
representation. That said, to gain intution for what this object is,
it makes sense to pick some explicit bases and write down examples.

To that end, suppose that B=ℝ. Then `T¹B=T¹ℝ` is just our usual notion of a dual
number, i.e. for some element `η ∈ T¹ℝ`, we may consider `η = a + bϵ` for
real numbers `(a, b)` and `ϵ` an infinitessimal differential such that `ϵ^2 = 0`.

Equivalently, we may think of `η` as being identified with the vector
`(a, b) ∈ ℝ²` with some additional structure.
The same form essentially holds for general `B`, where we may write (as sets):

    T¹B = {(a, b) | a ∈ B, b ∈ Tₐ B }

Note that these vectors are orthogonal to those in the underlying base space.
For example, if `B=ℝ²`, then we have:

    T¹ℝ² = {([aₓ, a_y], [bₓ, b_y]) | [aₓ, a_y] ∈ ℝ², [bₓ, b_y] ∈ Tₐ ℝ² }

For convenience, we will sometimes writes these in one as:

    η ∈ T ℝ²  = aₓ x̂ + a_y ŷ + bₓ ∂/∂x|_aₓ x + b_y ∂/∂y|_{a_y}
             := aₓ x̂ + a_y ŷ + bₓ ∂/∂¹x + b_y ∂/∂¹y
             := [aₓ, a_y] + [bₓ, b_y] ∂/∂¹
             := a + b ∂/∂¹
             := a + b ∂₁

These are all definitional equivalences and we will mostly work with the final
form. An important thing to keep in mind though is that the subscript on ∂₁ does
not refer to a dimension of the underlying base manifold (for which we will
rarely pick an explicit basis here), but rather tags the basis of the tangent
bundle.

Let us iterate this construction to second order. We have:

    T²B = T¹(T¹B) = { (α, β) | α ∈ T¹B, β ∈ T_α T¹B }
                  = { ((a, b), (c, d)) | a ∈ B, b ∈ Tₐ B, c ∈ Tₐ B, d ∈ T²ₐ B}

(where in the last equality we used the linearity of the tangent vector).

Following our above notation, we will canonically write such an element as:

      a + b ∂₁ + c ∂₂ + d ∂₂ ∂₁
    = a + b ∂₁ + c ∂₂ + d ∂₁ ∂₂

It is worth noting that there still only is one base point `a` of the
underlying manifold and thus `TⁿB` is a vector bundle over `B` for all `N`.

# Further Reading

[1] https://en.wikipedia.org/wiki/Tangent_bundle
"""
abstract type AbstractTangentBundle{N, B}; end
basespace(::Type{<:AbstractTangentBundle{N, B} where N}) where {B} = B

const ATB = AbstractTangentBundle

abstract type TangentIndex; end

struct CanonicalTangentIndex <: TangentIndex
    i::Int
end

struct TaylorTangentIndex <: TangentIndex
    i::Int
end

function Base.getindex(a::AbstractTangentBundle, b::TaylorTangentIndex)
    error("$(typeof(a)) is not taylor-like. Taylor indexing is ambiguous")
end

"""
    struct TangentBundle{N, B, P}

A fully explicit coordinate representation of the tangent bundle.
Represented by a primal value in `B` and a vector of `2^(N-1)` partials.
"""
struct TangentBundle{N, B, P} <: AbstractTangentBundle{N, B}
    primal::B
    partials::P
end

check_tangent_invariant(lp, N) = @assert lp == 2^N - 1
@ChainRulesCore.non_differentiable check_tangent_invariant(lp, N)

function TangentBundle{N}(primal::B, partials::P) where {N, B, P}
    check_tangent_invariant(length(partials), N)
    TangentBundle{N, Core.Typeof(primal), P}(primal, partials)
end

function TangentBundle{N,B}(primal::B, partials::P) where {N, B, P}
    check_tangent_invariant(length(partials), N)
    TangentBundle{N, B, P}(primal, partials)
end

function Base.show(io::IO, x::TangentBundle)
    print(io, x.primal)
    print(io, " + ")
    print(io, x.partials[1], " ∂₁")
    length(x.partials) >= 2 && print(io, " + ", x.partials[2], " ∂₂")
    length(x.partials) >= 3 && print(io, " + ", x.partials[3], " ∂₁ ∂₂")
    length(x.partials) >= 4 && print(io, " + ", x.partials[4], " ∂₃")
    length(x.partials) >= 5 && print(io, " + ", x.partials[5], " ∂₁ ∂₃")
    length(x.partials) >= 6 && print(io, " + ", x.partials[6], " ∂₂ ∂₃")
    length(x.partials) >= 7 && print(io, " + ", x.partials[7], " ∂₁ ∂₂ ∂₃")
end

function Base.getindex(a::TangentBundle{N}, b::TaylorTangentIndex) where {N}
    if b.i === N
        return a.partials[end]
    end
    error("$(typeof(a)) is not taylor-like. Taylor indexing is ambiguous")
end

"""
    struct TaylorBundle{N, B, P}

The taylor bundle construction mods out the full N-th order tangent bundle
by the equivalence relation that coefficients of like-order basis elements be
equal, i.e. rather than a generic element

    a + b ∂₁ + c ∂₂ + d ∂₃ + e ∂₂ ∂₁ + f ∂₃ ∂₁ + g ∂₃ ∂₂ + h ∂₃ ∂₂ ∂₁

we have a tuple (c₀, c₁, c₂, c₃) corresponding to the full element

    c₀ + c₁ ∂₁ + c₁ ∂₂ + c₁ ∂₃ + c₂ ∂₂ ∂₁ + c₂ ∂₃ ∂₁ + c₂ ∂₃ ∂₂ + c₃ ∂₃ ∂₂ ∂₁

i.e.

    c₀ + c₁ (∂₁ + ∂₂ + ∂₃) + c₂ (∂₂ ∂₁ + ∂₃ ∂₁ + ∂₃ ∂₂) + c₃ ∂₃ ∂₂ ∂₁

This restriction forms a submanifold of the original manifold. The naming is
by analogy with the (truncated) Taylor series

    c₀ + c₁ x + 1/2 c₂ x² + 1/3! c₃ x³ + O(x⁴)

"""
struct TaylorBundle{N, B, P} <: AbstractTangentBundle{N, B}
    primal::B
    coeffs::P

    function TaylorBundle{N, B}(primal::B, coeffs::P) where {N, B, P}
        check_taylor_invariants(coeffs, primal, N)
        new{N, B, P}(primal, coeffs)
    end
end

function check_taylor_invariants(coeffs, primal, N)
    @assert length(coeffs) == N
    if isa(primal, TangentBundle)
        @assert isa(coeffs[1], TangentBundle)
    end
end
@ChainRulesCore.non_differentiable check_taylor_invariants(coeffs, primal, N)

function TaylorBundle{N}(primal, coeffs) where {N}
    TaylorBundle{N, Core.Typeof(primal)}(primal, coeffs)
end

Base.getindex(tb::TaylorBundle, tti::TaylorTangentIndex) = tb.coeffs[tti.i]
function Base.getindex(tb::TaylorBundle, tti::CanonicalTangentIndex)
    tb.coeffs[count_ones(tti.i)]
end

"""
    struct UniformBundle{N, B}

Represents an N-th order tangent bundle with all unform partials. Particularly
useful for representing singleton values.
"""
struct UniformBundle{N, B, U} <: AbstractTangentBundle{N, B}
    primal::B
    partial::U
end
UniformBundle{N, B, U}(primal::B) where {N,B,U} = UniformBundle{N, B, U}(primal, U.instance)
UniformBundle{N, B}(primal::B, partial::U) where {N,B,U} = UniformBundle{N, Core.Typeof(primal), U}(primal, partial)
UniformBundle{N}(primal, partial::U) where {N,U} = UniformBundle{N, Core.Typeof(primal), U}(primal, partial)
UniformBundle{N, <:Any, U}(primal, partial::U) where {N, U} = UniformBundle{N, Core.Typeof(primal), U}(primal, U.instance)
UniformBundle{N, <:Any, U}(primal) where {N, U} = UniformBundle{N, Core.Typeof(primal), U}(primal, U.instance)

const ZeroBundle{N, B} = UniformBundle{N, B, ZeroTangent}
const DNEBundle{N, B} = UniformBundle{N, B, NoTangent}

Base.getindex(u::UniformBundle, ::TaylorTangentIndex) = u.partial

"""
    TupleTangentBundle{N, B <: Tuple}

Represents the tagent bundle where the base space is some tuple type.
Mathematically, this tangent bundle is the product bundle of the individual
element bundles.
"""
struct CompositeBundle{N, B, T<:Tuple{Vararg{AbstractTangentBundle{N}}}} <: AbstractTangentBundle{N, B}
    tup::T
end
CompositeBundle{N, B}(tup::T) where {N, B, T} = CompositeBundle{N, B, T}(tup)

function Base.getindex(tb::CompositeBundle{N, B} where N, tti::TaylorTangentIndex) where {B}
    B <: SArray && error()
    Tangent{B}(map(tb.tup) do el
        el[tti]
    end...)
end

primal(b::CompositeBundle{N, <:Tuple} where N) = map(primal, b.tup)
function primal(b::CompositeBundle{N, T} where N) where T<:CompositeBundle
    T(map(primal, b.tup)...)
end
@generated primal(b::CompositeBundle{N, B} where N) where {B} =
    quote
        x = map(primal, b.tup)
        $(Expr(:splatnew, B, :x))
    end

expand_singleton_to_array(asize, a::AbstractZero) = fill(a, asize...)
expand_singleton_to_array(asize, a::AbstractArray) = a

function unbundle(atb::TangentBundle{Order, A}) where {Order, Dim, T, A<:AbstractArray{T, Dim}}
    asize = size(atb.primal)
    StructArray{TangentBundle{Order, T}}((atb.primal, map(a->expand_singleton_to_array(asize, a), atb.partials)...))
end

function StructArrays.staticschema(::Type{<:TangentBundle{N, B, T}}) where {N, B, T}
    Tuple{B, T.parameters...}
end

function StructArrays.component(m::TangentBundle{N, B, T}, i::Int) where {N, B, T}
    i == 1 && return m.primal
    return m.partials[i - 1]
end

function StructArrays.createinstance(T::Type{<:TangentBundle}, args...)
    T(first(args), Base.tail(args))
end

function unbundle(atb::TaylorBundle{Order, A}) where {Order, Dim, T, A<:AbstractArray{T, Dim}}
    StructArray{TaylorBundle{Order, T}}((atb.primal, atb.coeffs...))
end

function ChainRulesCore.rrule(::typeof(unbundle), atb::TaylorBundle)
    unbundle(atb), Δ->throw(Δ)
end

function StructArrays.staticschema(::Type{<:TaylorBundle{N, B, T}}) where {N, B, T}
    Tuple{B, T.parameters...}
end

function StructArrays.staticschema(::Type{<:TaylorBundle{N, B}}) where {N, B}
    Tuple{B, Vararg{Any, N}}
end

function StructArrays.component(m::TaylorBundle{N, B}, i::Int) where {N, B, T}
    i == 1 && return m.primal
    return m.coeffs[i - 1]
end

function StructArrays.createinstance(T::Type{<:TaylorBundle}, args...)
    T(first(args), Base.tail(args))
end

function unbundle(zb::ZeroBundle{N, A}) where {N,T,Dim,A<:AbstractArray{T, Dim}}
    StructArray{ZeroBundle{N, T}}((zb.primal, fill(zb.partial, size(zb.primal)...)))
end

function ChainRulesCore.rrule(::typeof(unbundle), atb::ZeroBundle)
    unbundle(atb), Δ->throw(Δ)
end

function StructArrays.createinstance(T::Type{<:ZeroBundle}, args...)
    T(args[1], args[2])
end

function rebundle(A::AbstractArray{<:TangentBundle{N}}) where {N}
    TangentBundle{N}(
        map(x->x.primal, A),
        ntuple(2^N-1) do i
            map(x->x.partials[i], A)
        end)
end

function rebundle(A::AbstractArray{<:TaylorBundle{N}}) where {N}
    TaylorBundle{N}(
        map(x->x.primal, A),
        ntuple(N) do i
            map(x->x.coeffs[i], A)
        end)
end

function ChainRulesCore.rrule(::typeof(rebundle), atb)
    rebundle(atb), Δ->throw(Δ)
end


text:
module Markets

using Genie, Logging, LoggingExtras

function main()
  Core.eval(Main, :(const UserApp = $(@__MODULE__)))

  Genie.genie(; context = @__MODULE__)

  Core.eval(Main, :(const Genie = UserApp.Genie))
  Core.eval(Main, :(using Genie))
end

end


text:
module Epi

using NamedArrays, Distributions

include("Utils.jl")
include("crossTable.jl")
include("ORPRRR.jl")
include("chisquare.jl")

export crossTable
export removeMissing, removeNaN, unlist
export oddsratio, riskratio
export chisquaretest

end


text:
MOIU.@model(DualizableModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan,),
            (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
             MOI.ExponentialCone, MOI.DualExponentialCone,
             MOI.PositiveSemidefiniteConeTriangle,),
            (MOI.PowerCone, MOI.DualPowerCone),
            (),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

mutable struct PrimalDualMap{T}
    primal_var_dual_con::Dict{VI, CI}
    primal_con_dual_var::Dict{CI, Vector{VI}}
    primal_con_dual_con::Dict{CI, CI}
    primal_con_constants::Dict{CI, Vector{T}}

    function PrimalDualMap{T}() where T
        return new(Dict{VI, CI}(),
                   Dict{CI, Vector{VI}}(),
                   Dict{CI, CI}(),
                   Dict{CI, Vector{T}}())
    end
end

function is_empty(primal_dual_map::PrimalDualMap{T}) where T
    if isempty(primal_dual_map.primal_var_dual_con) &&
       isempty(primal_dual_map.primal_con_dual_var) &&
       isempty(primal_dual_map.primal_con_dual_con) &&
       isempty(primal_dual_map.primal_con_constants)
       return true
    end
    return false
end

function empty!(primal_dual_map::PrimalDualMap{T}) where T
    primal_dual_map.primal_var_dual_con = Dict{VI, CI}()
    primal_dual_map.primal_con_dual_var = Dict{CI, Vector{VI}}()
    primal_dual_map.primal_con_dual_con = Dict{CI, CI}()
    primal_dual_map.primal_con_constants = Dict{CI, Vector{T}}()
end

struct DualProblem{T, OT <: MOI.ModelLike}
    dual_model::OT #It can be a model from an optimizer or a DualizableModel{T}
    primal_dual_map::PrimalDualMap{T}

    function DualProblem{T}(dual_optimizer::OT, pdmap::PrimalDualMap{T}) where {T, OT <: MOI.ModelLike}
        return new{T, OT}(dual_optimizer, pdmap)
    end
end

function DualProblem{T}(dual_optimizer::OT) where {T, OT <: MOI.ModelLike}
    return DualProblem{T}(dual_optimizer, PrimalDualMap{T}())
end

function DualProblem(dual_optimizer::OT) where {OT <: MOI.ModelLike}
    return DualProblem{Float64}(dual_optimizer)
end

# Empty DualProblem cosntructor
function DualProblem{T}() where {T}
    return DualProblem{T}(DualizableModel{T}(), PrimalDualMap{T}())
end


text:
using CategoricalArrays

function crammerv(x::CategoricalArray, y::CategoricalArray)
	@assert length(x) == length(y)
	lx = levelcode.(x)
	ly = levelcode.(y)
	cx = countmap(lx)
	cy = countmap(ly)
	ux = unique(lx)
	uy = unique(ly)
	k = length(unique(x))
	r = length(unique(y))
	cxy = countmap(zip(lx, ly))
	n = length(x)
	chisq = sum(begin
			nume = (get(cxy, (i, j), 0) - cx[i] * cy[j] / n)^2
			deno = (cx[i] * cy[j] / n)
			nume / deno
		end for (i, j) in Iterators.product(ux, uy))
	sqrt(chisq / n / min(k - 1, r - 1))
end
crammerv(x::CategoricalArray, y) = crammerv(x, bin_number_vector(y, 10))
crammerv(x, y::CategoricalArray) = crammerv(bin_number_vector(x, 10), y)

function bin_number_vector(X, k=10)
	minx = minimum(X)
	maxx = maximum(X)
	bins = range(start=minx, stop=maxx, length=k)
	map(x -> findfirst(x .< bins), X)
end

function compute_cordf(df)
	result = DataFrame(col=names(df))
	function compute_cor(name_1, name_2)
		col_1 = df[!, name_1]
		col_2 = df[!, name_2]
		if eltype(col_1) <: Number && eltype(col_2) <: Number
			StatsBase.cor(col_1, col_2)
		else
			crammerv(categorical(col_1), categorical(col_2))
		end
	end
	transformations = map(names(df)) do name_1
		:col => ByRow(name_2 -> compute_cor(name_1, name_2)) => Symbol(name_1)
	end
	transformations
	transform(result, transformations...)
end

function find_multicolinearity(cordf, ycol; threshold=0.7)
	to_be_removed = String[]
	sizehint!(to_be_removed, size(cordf, 1)÷3)
	considers = setdiff(cordf.col, ["SalePrice"])
	name_index = String.(cordf.col)
	df = select(cordf, Not("col"))
	for n1 in considers, n2 in considers
		if n1 in to_be_removed || n2 in to_be_removed || n1 == n2
			continue
		end
		row_ind = findfirst(n1 .== name_index)
		col_ind = findfirst(n2 .== name_index)
		if abs(df[row_ind, col_ind]) > threshold
			append!(to_be_removed, [abs(df[row_ind, ycol]) > abs(df[col_ind, ycol]) ? n2 : n1])
		end
	end
	to_be_removed
end


text:
@testset "gc2gd" begin
    xyz = [2.e6, 3.e6, 5.244e6]
    e, p, h = gc2gd(WGS84, xyz)
    @test e ≈ 0.98279372324732907 atol=1e-14
    @test p ≈ 0.97160184819075459 atol=1e-14
    @test h ≈ 331.41724614260599 atol=1e-8
    e, p, h = gc2gd(GRS80, xyz)
    @test e ≈ 0.98279372324732907 atol=1e-14
    @test p ≈ 0.97160184820607853 atol=1e-14
    @test h ≈ 331.41731754844348 atol=1e-8
    e, p, h = gc2gd(WGS72, xyz)
    @test e ≈ 0.98279372324732907 atol=1e-14
    @test p ≈ 0.97160181811015119 atol=1e-14
    @test h ≈ 333.27707261303181 atol=1e-8
    @test_throws ArgumentError gc2gd(WGS84, xyz[1:2])
end

@testset "gc2gde" begin
    a = 6378136.0
    f = 0.0033528
    xyz = [2e6, 3e6, 5.244e6]
    e, p, h = gc2gde(a, f, xyz)
    @test e ≈ 0.98279372324732907 atol=1e-14
    @test p ≈ 0.97160183775704115 atol=1e-14
    @test h ≈ 332.36862495764397 atol=1e-8
    @test_throws ArgumentError gc2gde(a, f, xyz[1:2])
    @test_throws ERFAException gc2gde(-100, f, xyz)
    @test_throws ERFAException gc2gde(a, -100, xyz)
end

@testset "gd2gc" begin
    e = 3.1
    p = -0.5
    h = 2500.0
    xyz = gd2gc(WGS84, e, p, h)
    @test xyz[1] ≈ -5599000.5577049947 atol=1e-7
    @test xyz[2] ≈ 233011.67223479203 atol=1e-7
    @test xyz[3] ≈ -3040909.4706983363 atol=1e-7
    xyz = gd2gc(GRS80, e, p, h)
    @test xyz[1] ≈ -5599000.5577260984 atol=1e-7
    @test xyz[2] ≈ 233011.6722356703 atol=1e-7
    @test xyz[3] ≈ -3040909.4706095476 atol=1e-7
    xyz = gd2gc(WGS72, e, p, h)
    @test xyz[1] ≈ -5598998.7626301490 atol=1e-7
    @test xyz[2] ≈ 233011.5975297822 atol=1e-7
    @test xyz[3] ≈ -3040908.6861467111 atol=1e-7
end

@testset "gd2gce" begin
    a = 6378136.0
    f = 0.0033528
    e = 3.1
    p = -0.5
    h = 2500.0
    xyz = gd2gce(a, f, e, p, h)
    @test xyz[1] ≈ -5598999.6665116328 atol=1e-7
    @test xyz[2] ≈ 233011.63514630572 atol=1e-7
    @test xyz[3] ≈ -3040909.0517314132 atol=1e-7
end

@testset "gmst00" begin
    g = gmst00(2400000.5, 53736.0, 2400000.5, 53736.0)
    @test g ≈ 1.754174972210740592 atol=1e-14
end

@testset "gmst06" begin
    g = gmst06(2400000.5, 53736.0, 2400000.5, 53736.0)
    @test g ≈ 1.754174971870091203 atol=1e-14
end

@testset "gmst82" begin
    g = gmst82(2400000.5, 53736.0)
    @test g ≈ 1.754174981860675096 atol=1e-14
end

@testset "gst00a" begin
    g = gst00a(2400000.5, 53736.0, 2400000.5, 53736.0)
    @test g ≈ 1.754166138018281369 atol=1e-14
end

@testset "gst00b" begin
    g = gst00b(2400000.5, 53736.0)
    @test g ≈ 1.754166136510680589 atol=1e-14
end

@testset "gst06" begin
    rnpb = [0.9999989440476103608 -0.1332881761240011518e-2 -0.5790767434730085097e-3;
            0.1332858254308954453e-2 0.9999991109044505944 -0.4097782710401555759e-4;
            0.5791308472168153320e-3 0.4020595661593994396e-4 0.9999998314954572365]
    g = gst06(2400000.5, 53736.0, 2400000.5, 53736.0, rnpb)
    @test g ≈ 1.754166138018167568 atol=1e-14
    @test_throws ArgumentError gst06(2400000.5, 53736.0, 2400000.5, 53736.0, rnpb[1:2,:])
end

@testset "gst06a" begin
    g = gst06a(2400000.5, 53736.0, 2400000.5, 53736.0)
    @test g ≈ 1.754166137675019159 atol=1e-14
end

@testset "gst94" begin
    g = gst94(2400000.5, 53736.0)
    @test g ≈ 1.754166136020645203 atol=1e-14
end

@testset "g2icrs" begin
    dr, dd = g2icrs(5.5850536063818546461558105, -0.7853981633974483096156608)
    @test dr ≈ 5.9338074302227188048671 atol=1e-14
    @test dd ≈ -1.1784870613579944551541 atol=1e14
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Two Body orbit propagator.
#
#       This algorithm considers a perfect Keplerian orbit. In other words, no
#       perturbation is considered during the propagation and the Earth is
#       modeled as a perfect sphere.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export TwoBody_Structure
export twobody_init, twobody!

################################################################################
#                                  Overloads
################################################################################

# Copy for TwoBody_Structure.
Base.copy(m::TwoBody_Structure) =
    TwoBody_Structure(m.epoch, m.a_0, m.e_0, m.i_0, m.Ω_0, m.ω_0, m.M_0, m.Δt,
                      m.a, m.M_k, m.f_k, m.μ)

################################################################################
#                                  Functions
################################################################################

"""
    function twobody_init(epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number, μ::T) where T

Initialize the data structure of Two Body orbit propagator algorithm.

# Args

* `epoch`: Epoch of the orbital elements [s].
* `a_0`: Initial semi-major axis [m].
* `e_0`: Initial eccentricity.
* `i_0`: Initial inclination [rad].
* `Ω_0`: Initial right ascension of the ascending node [rad].
* `ω_0`: Initial argument of perigee [rad].
* `f_0`: Initial true anomaly.
* `μ`: Standard gravitational parameter of the central body [m^3/s^2].

# Returns

The structure `TwoBody_Structure` with the initialized parameters.

"""
function twobody_init(epoch::Number, a_0::Number, e_0::Number, i_0::Number,
                      Ω_0::Number, ω_0::Number, f_0::Number, μ::T) where T

    # The propagator is only defined for 0 <= e < 1.
    if !(0 <= e_0 < 1)
        throw(ArgumentError("The Two Body propagator only supports eccentricities in the interval [0,1)"))
    end

    # Compute the mean motion using the semi-major axis.
    n_0 = sqrt(μ/a_0^3)

    # Compute the initial mean anomaly.
    M_0 = f_to_M(e_0, f_0)

    # Create and return the Two Body orbit propagator structure.
    TwoBody_Structure{T}(epoch, a_0, e_0, i_0, Ω_0, ω_0, M_0, f_0, 0, n_0, M_0, f_0, μ)
end

"""
    function twobody!(tbd::TwoBody_Structure, t::Number)

Propagate the orbit defined in `tbd` (see `TwoBody_Structure`) until the time
`t` [s]. Notice that the values in `tbd` will be modified.

# Returns

* The position vector represented in the inertial frame at time `t` [m].
* The velocity vector represented in the inertial frame at time `t` [m/s]

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. If the orbit parameters are obtained
from a TLE, then the inertial frame will be TEME.

"""
function twobody!(tbd::TwoBody_Structure, t::Number)
    # Time elapsed since epoch.
    Δt     = t
    tbd.Δt = Δt

    # Update the mean anomaly.
    tbd.M_k = tbd.M_0 + tbd.n_0*Δt

    # Convert the mean anomaly to true anomaly.
    tbd.f_k = M_to_f(tbd.e_0, tbd.M_k)

    # Compute the position and velocity vectors given the orbital elements.
    (r_i_k, v_i_k) =
        kepler_to_rv(tbd.a_0, tbd.e_0, tbd.i_0, tbd.Ω_0, tbd.ω_0, tbd.f_k)

    # Return the position and velocity vector represented in the inertial
    # reference frame.
    (r_i_k, v_i_k)
end


text:
type MoRMax<:FiniteHorizonAgent
    policy::Matrix{UInt64}

    ns::Array{UInt64, 4}
    n::Array{UInt64, 3}

    phat::Array{Float64, 4}
    rhat::Array{Float64, 3}

    pprime::Array{Float64, 4}
    rprime::Array{Float64, 3}

    rewards::Array{Float64, 3}

    τ::Array{Int64, 3}
    τstar::Int64

    Qhat::Array{Float64, 3}
    Qprime::Array{Float64, 3}
    policyprime::Matrix{UInt64}
    Qprev::Array{Float64, 3}

    δ::Float64
    ϵ::Float64
    rewards_known::Bool
    maxRet::Float64
    mfactor::Float64
end

function MoRMax(S, A, H, rewards, δ, ϵ, mfactor=1.)
    policy = zeros(UInt64, S, H)
    policyprime = copy(policy)

    ns = zeros(UInt64, S, S, A, H)
    n = zeros(UInt64, S, A, H)

    phat = zeros(S, S, A, H)
    pprime = copy(phat)

    maxrs = maximum(reshape(rewards, S*A, H), 1)
    rhat = copy(rewards)
    for t=1:H
        rhat[:, :, t] = sum(maxrs[t:end])
    end
    rprime = copy(rhat)
    τ = zeros(Int64, S, A, H)

    Qhat = zeros(S, A, H)
    Qprime = copy(Qhat)
    Qprev = zeros(S, A, H)
    compute_optQ!(Qhat, policy, phat, rhat)
    rand!(policy, 1:A) # Initialize with uniformly random policy

    MoRMax(policy, ns,n,  phat, rhat, pprime, rprime, rewards,
           τ, 1, Qhat, Qprime, policyprime, Qprev, δ, ϵ, true, sum(maxrs), mfactor)
end

maxV(m::MoRMax) = m.maxRet
function mval(m::MoRMax)
    S = nS(m)
    A = nA(m)
    H = horizon(m)
    mval = 1800*maxV(m)^2 / m.ϵ^2 * H^2
    mval *= log(144 * H * S^2 * A^2 * maxV(m) / m.δ / m.ϵ)
    mval
end

function update_policy!(m::MoRMax)
    S = nS(m)
    A = nA(m)
    H = horizon(m)
    mvalue = mval(m)
    if !m.rewards_known
        error("Not Implemented")
    end

    for t ∈ 1:H, a ∈ 1:A, s ∈ 1:S
        if m.n[s,a,t] >= mvalue * m.mfactor

            if (m.τ[s, a, t] == 0 ||
                (m.τstar > m.τ[s,a,t] && m.Qprev[s, a, t] > m.ϵ / H + m.Qhat[s, a, t]))
                    copy!(m.pprime, m.phat)
                    copy!(m.rprime, m.rhat)
                    m.pprime[:, s, a, t] .= m.ns[:, s, a, t] / m.n[s, a, t]
                    m.rprime[s, a, t] = m.rewards[s, a, t]

                    #println("Propose model update for s=$s, a=$a, t=$t")
                    compute_optQ!(m.Qprime, m.policyprime, m.pprime, m.rprime)
                    if m.Qprime[s, a, t] <= m.Qhat[s, a, t]
                        m.phat[:, s, a, t] = m.pprime[:, s, a, t]
                        m.rhat[s, a, t] = m.rprime[s, a, t]
                        copy!(m.Qhat, m.Qprime)
                        copy!(m.policy, m.policyprime)
                        #println("Accepted model update for s=$s, a=$a, t=$t")
                    end
                    m.τstar += 1
                    m.τ[s,a,t] = m.τstar
                    m.Qprev[s,a,t] = m.Qhat[s,a,t]
            end

            m.n[s,a,t] = 0
            m.ns[:, s, a, t] = 0
        end
    end
end

function compute_optQ!(Q, pol, P, R)
    S, A, H = size(R)
    V = zeros(S)
    for t ∈ H:-1:1
        Q[:, :, t] .= R[:, :, t] + squeeze(sum(V .* P[:, :, :, t], 1), 1)
        for s ∈ 1:S
            a = indmax(Q[s, :, t])
            pol[s, t] = a
            V[s] = Q[s, a, t]
        end
    end
end


text:
"""
    module Mocking

This module defines the @mock_test macro, which deepcopies a new module in which you can
execute arbitrary setup before running tests. This is useful because it allows you to
redefine behavior without worrying about corrupting the remaining tests.
"""
module Mocking
import DeepcopyModules

export @mock_test

"""
    @mock_test(m::Module, begin ... end, @testset begin ... end)
This macro deepcopies the provided Module, executes the provided setup code in that Module,
and then executes the `test` code in that new module.
This can be used to override behavior of certain parts of your code without worrying about
corrupting behavior in future tests.
# Examples
```julia
@mock_test(
    MyModule,
    begin
        # Override behavior of expensive_code
        struct TestSucceeded end
        expensive_code() = throw(TestSucceeded())
    end,
    @testset "expensive_code is correctly triggered" begin
        @test_throws TestSucceeded MyPackage.parse_args("...")
    end)
)
```
"""
macro mock_test(_module, setup, test)
    quote
        m2 = $DeepcopyModules.deepcopy_module($(esc(_module)))

        @eval m2 eval($(Expr(:quote, setup)))
        @eval m2 eval($(Expr(:quote, test)))
    end
end

end