text:
#module CenterOfMass

using PyPlot
using Convex
using SCS
using GLPKMathProgInterface
include("Histogram.jl")
include("plotting_tools.jl")

# all wrenches computed around the origin
point_to_wrench(p::Array{Float64,1}) = [1.0;-p[2]; p[1]]

# normal(x, var) = exp(-0.5*x*x/var)/sqrt(2*pi*var)
const one_over_sqrt_2pi = 1.0 / sqrt(2*pi)
function normal(x, var::Float64)
  ov = 1.0/var
  sqrt2pivar = sqrt(ov) * one_over_sqrt_2pi
  normal(x, -0.5*ov, sqrt2pivar)
end

normal(x::Float64, nhalfov::Float64, ov2sqrt2pi::Float64) = exp(x*x*nhalfov) * ov2sqrt2pi

function critical_force(wrench_applied, wrench_offset, W_boundary)
  fa = Variable(1)
  Fr = Variable(size(W_boundary, 2))

  stability = wrench_applied * fa +
              W_boundary * Fr +
              wrench_offset == 0
  feasibility = [fa >= 0; Fr >= 0]

  problem = maximize(fa, [stability; feasibility])

  solver = SCSSolver(verbose = 0)

  solve!(problem, solver)

  problem.optval, evaluate(Fr)
end

function solve_minimal(a, w1, w2, b)
  A = [a w1 w2]

  #if abs(det(A)) <= 1e-3
    #return -1, [-1;-1]
  #end

  x = A \ (-b)

  for ii = 1:3
    if x[ii] < -1e-3
      return -1.0
    end
  end

  return x[1]
end

# assume points form the convex hull and are arranged clockwise or
# counter-clockwise
function critical_force_iterative(wrench_applied, wrench_offset, W_boundary,
  actuator_limit::Float64 = 1e6)

  max_val::Float64 = -1

  n = size(W_boundary, 2)

  # TODO: implement binary search
  for ii = 1:n
    i2 = mod(ii, n) + 1
    val = solve_minimal(wrench_applied, W_boundary[:,ii], W_boundary[:,i2], wrench_offset)

    if val > max_val
      if val > actuator_limit
        max_val = actuator_limit
      else
        max_val = val
      end

      break
    end
  end

  max_val
end

flatten(v::Array) = hcat(v...)
flatten(v) = v

function rand_in_circle()
  r = sqrt(rand())
  theta = 2*pi*rand()
  r * [cos(theta);sin(theta);0]
end

function to_wrench_matrix(points)
  n = length(points)
  ws = zeros(3, n)
  for ii = 1:n
    ws[:,ii] = point_to_wrench(points[ii])
  end

  ws
end

function critical_force_from_points(boundary_ps, com_p, applied_p, mass,
    actuator_limit = 1e6, solver = critical_force_iterative)

  g = -9.8

  boundary_ws = to_wrench_matrix(boundary_ps)
  gravity_w = mass*g*point_to_wrench(com_p)
  applied_w = point_to_wrench(applied_p)

  fa = solver(applied_w, gravity_w, boundary_ws, actuator_limit)

  #stability = applied_w * fa +
              #boundary_ws * f_boundary +
              #gravity_w

  fa
end

function get_reaction_points(boundary_ps, boundary_fs)
  inds = find(boundary_fs .>= 1e-3)
  force_points = boundary_ps[inds]

  force_points
end

function plot_solution(applied_p)
  scatter3D(applied_p[1,:]', applied_p[2,:]', [0], color="k", s=800, alpha=0.8)
end

function plot_attachment_points(attachment_ps)
  for p = attachment_ps
    #scatter3D(p[1,:]', p[2,:]', [0], s=400, facecolors="",edgecolors = "k",
    #alpha = 0.8)
    plot_circle(p, radius = 0.2)
  end
end

function plot_new_point(point)
  # nothing to see here
end

function plot_occupied_points(occupied_points)
  for p = occupied_points
    plot_quadrotor(p+[0;0;0.025], color="k", scale=0.15)
  end
end

function plot_measurement_points(measurement_points)
  for p = measurement_points
    plot_quadrotor(p+[0;0;0.025], color="r", scale=0.15)
  end
end

function initialize_prior(boundary_ps, resolution, interior_q, masses)
  boundary_vec = flatten(boundary_ps)

  min_p = minimum(boundary_vec, 2)[1:2]
  max_p = maximum(boundary_vec, 2)[1:2]
  position_ranges = map(x->x[1]:resolution:x[2], zip(min_p, max_p))
  ranges = (position_ranges..., masses)

  prior = Histogram(ranges)

  num_inside = 0
  for kk = 1:length(1:length(get_range(prior, 3)))
    for jj = 1:length(1:length(get_range(prior, 2)))
      for ii = 1:length(1:length(get_range(prior, 1)))
        p = from_indices(prior, [ii,jj,kk])
        if interior_q(p[1:2])::Bool
          num_inside += 1
        end
      end
    end
  end
  p_interior = (1 / (num_inside))

  for kk = 1:length(1:length(get_range(prior, 3)))
    for jj = 1:length(1:length(get_range(prior, 2)))
      for ii = 1:length(1:length(get_range(prior, 1)))
        p = from_indices(prior, [ii,jj,kk])
        value = 0
        if interior_q(p[1:2])::Bool
          value = p_interior
        end
        get_data(prior)[ii,jj,kk] = value
      end
    end
  end

  prior
end

function update_prior!(prior, interior_q, boundary_ws, applied_w, f_hat, sigma, actuator_limit = 1e6)
  critical_fs = get_critical_values(boundary_ws, applied_w, prior, interior_q, actuator_limit)
  update_prior!(prior, critical_fs, f_hat, sigma)
end

shift_dims(field) = reshape(field[1,:,:], (size(field, 2), size(field, 3)))

function update_prior!(prior, critical_fs, f_hat, sigma)
  p_prior = get_data(prior)

  p_f_given_com = map(y->normal(y, sigma^2), f_hat - critical_fs)

  p_f_and_com = p_f_given_com .* p_prior

  p_f = sum(p_f_and_com)

  p_com_given_f = p_f_and_com / p_f

  get_data(prior)[:] = p_com_given_f[:]

  assert(abs(sum(get_data(prior)) - 1.0) < 1e-6)

  prior
end

function get_critical_values(boundary_ws, applied_w, prior, interior_q, actuator_limit = 1e6)
  values = zeros(size(prior))

  for kk = 1:length(get_range(prior, 3))
    for jj = 1:length(get_range(prior, 2))
      for ii = 1:length(get_range(prior, 1))
        p = from_indices(prior, [ii, jj, kk])
        if interior_q(p[1:2])::Bool
          mass = p[3]
          gravity_w = mass*g*point_to_wrench(p)
          applied_f = critical_force_iterative(applied_w, gravity_w, boundary_ws, actuator_limit)
          values[ii,jj,kk] = applied_f
          if isnan(applied_f)
            values[ii,jj,kk] = 0.0
          end
        end
      end
    end
  end

  values
end

function plot_attachment_csqmis(boundary_ps, attachment_ps, csqmis)
  attachment_vec = hcat(attachment_ps...)

  csqmi_scaled = 800 * csqmis / maximum(csqmis)
  scatter3D(attachment_vec[1,:]', attachment_vec[2,:]',
    zeros(size(attachment_vec, 2)), s = csqmi_scaled,
    color="b", alpha = 0.5, edgecolor="k", linewidth=2)
end

function to_cloud(prior, interior_q)
  out = Any[]
  ind = 1
  for kk = 1:length(get_range(prior, 3))
    for jj = 1:length(get_range(prior, 2))
      for ii = 1:length(get_range(prior, 1))
        p = from_indices(prior, [ii, jj, kk])
        if(interior_q(p[1:2]))::Bool
          mass = p[3]
          probability = get_data(prior)[ii, jj, kk]
          push!(out, [p; probability])
        end
      end
    end
  end
  out = hcat(out...)

  out
end

function plot_field(field, title_string="")
  fig, ax = plt[:subplots](1)
  imshow(field' , cmap= "BuPu",
         interpolation="nearest", origin="lower")
  fig[:axes][1][:get_yaxis]()[:set_visible](false)
  fig[:axes][1][:get_xaxis]()[:set_visible](false)
  title(title_string)
  axis("off")
end

include("csqmi.jl")

typealias CombinationMap Dict{Array{Int64,1},Array{Float64,2}}
function generate_over_combinations(action_ws, num_choose, prior,
  boundary_ws, actuator_limit = 1e6, interior_q = x->true)

  total_limit = actuator_limit * num_choose

  d = CombinationMap()
  for c = combinations(1:length(action_ws), num_choose)
    action_w = mean(action_ws[c])
    field = get_critical_values(boundary_ws, action_w, prior, interior_q, total_limit)
    normals = normal_matrix(field[:], sigma^2)
    d[c] = normals
  end
  d
end

function generate_combinations_maps(action_ws, max_choose, prior,
  boundary_ws, actuator_limit = 1e6, interior_q = x->true)
  out = Array{CombinationMap,1}()
  for ii = 1:max_choose
    println("Generating normals for $(ii)")
    @time combination_map = generate_over_combinations(action_ws, ii, prior,
      boundary_ws, actuator_limit, interior_q)
    push!(out, combination_map)
  end
  out
end

function maximize_csqmi_combinations(robot_indices, num_select, prior,
  boundary_ws, action_ws, sigma, interior_q, actuator_limit)

  max_csqmi = 0.0
  best_combination = Array{Float64,1}()
  data = get_data(prior)[:]

  total_limit = actuator_limit * num_select

  for c = combinations(robot_indices, num_select)
    action_w = mean(action_ws[c])

    field = get_critical_values(boundary_ws, action_w, prior, interior_q, total_limit)
    normals = normal_matrix(field[:], sigma^2)

    csqmi = compute_mutual_information(data, normals)

    if csqmi > max_csqmi
      max_csqmi = csqmi
      best_combination = c
    end
  end

  max_csqmi, best_combination
end

function maximize_csqmi_available_robots(robot_indices, prior, boundary_ws,
  action_ws, sigma, interior_q, actuator_limit)

  max_csqmi = 0.0
  best_combination = Array{Float64,1}()

  for ii = 1:length(robot_indices)
    csqmi, c = maximize_csqmi_combinations(robot_indices, ii, prior,
      boundary_ws, action_ws, sigma, interior_q, actuator_limit)

    if csqmi > max_csqmi
      max_csqmi = csqmi
      best_combination = c
    end
  end

  max_csqmi, best_combination
end

function maximize_csqmi_additional_robot(robot_indices, prior, boundary_ws,
  action_ws, sigma, interior_q, actuator_limit)

  max_csqmi = 0.0
  best_combination = Array{Float64,1}()
  additional_robot = 0

  all_robots = collect(1:length(action_ws))

  remaining_robots = setdiff(all_robots, robot_indices)

  # iterate over all combinations of robots including the additional robot
  for remaining_robot = remaining_robots
    for ii = 0:length(robot_indices)
      for c = combinations(robot_indices, ii)
        robots = [remaining_robot; c]
        total_limit = actuator_limit * length(robots)

        action_w = mean(action_ws[robots])

        field = get_critical_values(circle_ws, action_w, prior, interior_q, total_limit)
        normals = normal_matrix(field[:], sigma^2)

        csqmi = compute_mutual_information(get_data(prior)[:], normals)

        if csqmi > max_csqmi
          max_csqmi = csqmi
          best_combination = robots
          additional_robot = remaining_robot
        end
      end
    end
  end

  max_csqmi, best_combination, additional_robot
end

function maximize_csqmi_additional_robot_feasibility(robot_indices, prior, boundary_ws,
  action_ws, sigma, interior_q, actuator_limit, feasibility_constraint)

  max_csqmi = 0.0
  best_combination = Array{Float64,1}()
  additional_robot = 0

  all_robots = collect(1:length(action_ws))

  remaining_robots = setdiff(all_robots, robot_indices)

  # iterate over all combinations of robots including the additional robot
  for remaining_robot = remaining_robots
    robot_indices_p = [robot_indices; remaining_robot]
    robots = action_ws[robot_indices_p]
    doubly_remaining = action_ws[setdiff(all_robots, robot_indices_p)]

    fp = feasibility_probability(belief, robots, doubly_remaining, actuator_limit, max_robots)

    if fp >= feasibility_constraint
      println("adding robot $(remaining_robot) is feasible: $(fp) >= $(feasibility_constraint)")
      for ii = 0:length(robot_indices)
        for c = combinations(robot_indices, ii)
          robots = [remaining_robot; c]
          total_limit = actuator_limit * length(robots)

          action_w = mean(action_ws[robots])

          field = get_critical_values(circle_ws, action_w, prior, interior_q, total_limit)
          normals = normal_matrix(field[:], sigma^2)

          csqmi = compute_mutual_information(get_data(prior)[:], normals)

          if csqmi > max_csqmi
            max_csqmi = csqmi
            best_combination = robots
            additional_robot = remaining_robot
          end
        end
      end
    else
      println("adding robot $(remaining_robot) isn't feasible: $(fp) < $(feasibility_constraint)")
    end
  end

  max_csqmi, best_combination, additional_robot
end

# feasibility
function lifting_feasibility(applied_wrenches, wrench_offset, actuator_limit)
  W = hcat(applied_wrenches...)

  fs = Variable(length(applied_wrenches))

  lifting_condition = W*fs + wrench_offset == 0

  feasibility = [fs >= 0.0; fs <= actuator_limit]

  problem = maximize(0, [lifting_condition; feasibility])

  solver = SCSSolver(verbose = 0)
  solve!(problem, solver)

  ret = problem.status == :Optimal

  ret
end

function check_feasible_configuration(chosen_ws, remaining_ws, offset_w, actuator_limit,
  max_num_robots)

  num_remaining = max_num_robots - length(chosen_ws)

  W_chosen = hcat(chosen_ws...)
  W_remaining = hcat(remaining_ws...)

  fs_chosen = Variable(length(chosen_ws))
  fs_remaining = Variable(length(remaining_ws))

  fs_enabled = Variable(length(remaining_ws), :Bin)

  feasibility = Convex.Constraint[]
  if length(remaining_ws) > 0 && length(chosen_ws) > 0
    #println("remaining and chosen")
    feasibility = [
      W_chosen*fs_chosen + W_remaining*fs_remaining + offset_w == 0.0;
      fs_chosen >= 0.0;
      fs_remaining >= 0.0;
      fs_chosen <= actuator_limit;
      fs_remaining <= fs_enabled*actuator_limit;
      sum(fs_enabled) <= num_remaining
    ]
  elseif length(remaining_ws) > 0
    #println("remaining")
    feasibility = [
      W_remaining*fs_remaining + offset_w == 0.0;
      fs_remaining >= 0.0;
      fs_remaining <= fs_enabled*actuator_limit;
      sum(fs_enabled) <= num_remaining
    ]
  else
    #println("chosen")
    feasibility = [
      W_chosen*fs_chosen + offset_w == 0.0;
      fs_chosen >= 0.0;
      fs_chosen <= actuator_limit;
    ]
  end

  problem = minimize(0.0, feasibility)

  # solver was giving incorrect results with presolve set to false
  solver = GLPKSolverMIP(presolve=false, msg_lev=GLPK.MSG_OFF)

  #TT=STDERR
  #out_read, out_write = redirect_stderr()
  #close(out_write)

  try
    solve!(problem, solver)
  catch ex
    #data = readavailable(out_read)
    #close(out_read)

    #redirect_stderr(TT)

    if isa(ex, InterruptException)
      println("Caught Interrupt")
    else
      println("Caught other exception")
      @show ex
    end
    throw(ex)
  end

  ret = problem.status == :Optimal

  ret
end

function milp()
  x = Variable(4, :Bin)

  problem = minimize(sum(x), x>=0.5)

  solver = GLPKSolverMIP(presolve=true)
  solve!(problem, solver)

  ret = problem.optval

  if problem.status == :Infeasible
    ret = 1000
  end

  ret
end

# probability of feasible configuration
function feasibility_probability(belief, chosen_ws, remaining_ws,
  actuator_limit, max_num_robots)
  g = -9.8
  ranges = get_range(belief)

  data = get_data(belief)

  total = 0.0

  for ii = 1:length(ranges[1])
    for jj = 1:length(ranges[2])
      for kk = 1:length(ranges[3])
        p = from_indices(belief, [ii;jj;kk])
        gravity_w = g * p[3] * point_to_wrench(p[1:2])

        feasible = check_feasible_configuration(chosen_ws, remaining_ws,
          gravity_w, actuator_limit, max_num_robots)

        total += feasible * data[ii,jj,kk]
      end
    end
  end

  total
end

# optimality

# expected configuration value

#end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
using Statistics
using Test
using EvoTrees

@testset "EvoTrees" begin

@info "Testing core"
include("core.jl")

@info "Testing MLJ"
include("MLJ.jl")

end

ext: jl


text:
using Documenter, SearchBall

makedocs(
    format = :html,
    # assets = [""],
    sitename = "SearchBall Documentation",
    authors = "Tobias Jagla",
    pages = [
        "Home" => "index.md",
        "GUI" => "gui.md",
        "Strategies" => "strategies.md",
        "Simulation" => "simulation.md",
        "Algorithms" => "algorithms.md"
    ]
)

# deploydocs(
#     repo = "github.com/control13/SearchBall.jl.git",
#     target = "build",
#     julia = "0.6",
#     deps   = nothing,
#     make   = nothing
# )


text:
# TODO: add simulation-modifying convenience wrappers to make bc's incoming. Without that, this module provides essentially no utility above SALT
"""
    module SaturableCPA

Convenience wrappers for solution of the nonlinear SatruableCPA equation.
Exports `HelmholtzSCPA` and `MaxwellSCPA` constructors and extends methods of
`NLsolve`, which must be separately imported.
"""
module SaturableCPA

export HelmholtzSCPA
export MaxwellSCPA

using ..Common
using ..Lasing
using NLsolve
using RecipesBase

import ..Lasing.SALTProblem
import ..Lasing.HelmholtzSALT
import ..Lasing.MaxwellSALT

################################################################################
# SCPA objects and constructors (mostly hooks into SALT constructors)

"""
    SCPAProblem{M,THS}

Saturable CPA object with `M` component fields
"""
struct SCPAProblem{NMODES,M,THS}
    SALT::THS

    SCPAProblem(SALT::SALTProblem{NMODES,M}) where {NMODES,M} = new{NMODES,M,typeof(SALT)}(SALT)
end

function Base.getproperty(mcpa::SCPAProblem, sym::Symbol)
    if sym == :SALT
        getfield(mcpa,:SALT)
    else
        getproperty(getfield(mcpa,:SALT),sym)
    end
end

Base.propertynames(::SCPAProblem,private=false) = propertynames(SALTProblem,private)

HelmholtzSCPA{NMODES} = SCPAProblem{NMODES,1}
"""
    HelmholtzSCPA(sim::Simulation, m::Integer) -> scpa

SPCA object with `m` fields for Helmholtz problem
"""
HelmholtzSCPA(sim::Simulation,m::Integer) = SCPAProblem(HelmholtzSALT{m}(sim))

MaxwellSCPA{NMODES} = SCPAProblem{NMODES,3}
"""
    MaxwellSCPA(sim::Simulation, m::Integer) -> scpa

SPCA object with `m` fields for Maxwell problem
"""
MaxwellSCPA(sim::Simulation,m::Integer) = SCPAProblem(MaxwellSALT{m}(sim))

################################################################################
# Extend NLsolve functions

fnames = (:nlsolve,:fixedpoint)
for fn ∈ fnames
    @eval NLsolve.$(fn)(ms::SCPAProblem; kwargs...) = $(fn)(ms, ms.ωs, ms.ψs; kwargs...)

    @eval begin function NLsolve.$(fn)(ms::SCPAProblem, ωs_init, ψs_init::VectorField; kwargs...)
            ms.m==size(ψs_init,2) || throw("number of modes in ψs_init ($(size(ψs_init,2))) must be the same as given in SCPA ($(ms.m))")
            return $(fn)(ms.SALT,ωs_init,ψs_init; kwargs...)
        end
    end
end

"""
    nlsolve(scpa; kwargs...) -> results

    nlsolve(scpa, init_ωs, init_ψs; kwargs...) -> results

solve nonlinear SCPA equation where `scpa` is a `HelmholtzSCPA` or `MaxwellSCPA`.
Results also stored `scpa`.

`init_ψs` must be a `ScalarField` (for Helmholtz) or `ElectricField` (for Maxwell)
"""
nlsolve

"""
    fixedpoint(scpa; kwargs...) -> results

    fixedpoint(scpa, init_ωs, init_ψs; kwargs...) -> results

Convenience wrapper for `NLsolve`'s `fixedpoint` wrapper.
Solves nonlinear SCPA equation where `scpa` is a `HelmholtzSCPA` or `MaxwellSCPA`.
Results also stored `scpa`.

`init_ψs` must be a `ScalarField` (for Helmholtz) or `ElectricField` (for Maxwell)
"""
fixedpoint

################################################################################
# Pretty Printing

import ..Common.PRINTED_COLOR_LIGHT

function Base.show(io::IO,ms::SCPAProblem)
    print(io,ms.m, " mode")
    ms.m>1 ? print(io,"s") : nothing
    if typeof(ms)<:HelmholtzSCPA
        printstyled(io," HelmholtzSCPA",color = PRINTED_COLOR_LIGHT)
    elseif typeof(ms)<:MAxwellSCPA
        printstyled(io," MaxwellSCPA",color = PRINTED_COLOR_LIGHT)
    end
    print(IOContext(io,:SCPA=>true), ms.SALT)
end

################################################################################
# Plotting
@recipe f(ms::SCPAProblem; by=abs2) = ms, by
@recipe f(by::Function, ms::SCPAProblem) = ms, by
@recipe f(ms::SCPAProblem, by::Function) = ms.SALT, by

@recipe f(sim::Simulation, ms::SCPAProblem; by=abs2) = sim, ms, by
@recipe f(ms::SCPAProblem, sim::Simulation; by=abs2) = sim, ms, by

@recipe f(ms::SCPAProblem, sim::Simulation, by::Function) = sim, ms, by
@recipe f(sim::Simulation, by::Function, ms::SCPAProblem) = sim, ms, by
@recipe f(ms::SCPAProblem, by::Function, sim::Simulation) = sim, ms, by
@recipe f(by::Function, sim::Simulation, ms::SCPAProblem) = sim, ms, by
@recipe f(by::Function, ms::SCPAProblem, sim::Simulation) = sim, ms, by

@recipe f(sim::Simulation, ms::SCPAProblem, by::Function) = sim, ms.SALT, by

end # module


text:
### This file contains convenience functions for constructing LoopSets.

# function strip_unneeded_const_deps!(ls::LoopSet)
#     for op ∈ operations(ls)
#         if isconstant(op) && iszero(length(reducedchildren(op)))
#             op.dependencies = NODEPENDENCY
#         end
#     end
#     ls
# end

function Base.copyto!(ls::LoopSet, q::Expr)
  q.head === :for || throw("Expression must be a for loop.")
  add_loop!(ls, q, 8)
  # strip_unneeded_const_deps!(ls)
end

function add_ci_call!(
  q::Expr,
  @nospecialize(f),
  args,
  syms,
  i,
  valarg = nothing,
  mod = nothing,
)
  call = if f isa Core.SSAValue
    Expr(:call, syms[f.id])
  else
    Expr(:call, f)
  end
  for arg ∈ @view(args[2:end])
    if arg isa Core.SSAValue
      push!(call.args, syms[arg.id])
    else
      push!(call.args, arg)
    end
  end
  if mod !== nothing # indicates it's `vmaterialize(!)`
    push!(call.args, Expr(:call, Expr(:curly, :Val, QuoteNode(mod))), valarg)
  end
  push!(q.args, Expr(:(=), syms[i], call))
end

function substitute_broadcast(
  q::Expr,
  mod::Symbol,
  inline::Bool,
  u₁::Int8,
  u₂::Int8,
  v::Int8,
  threads::Int,
  warncheckarg::Int,
)
  ci = first(Meta.lower(LoopVectorization, q).args).code
  nargs = length(ci) - 1
  ex = Expr(:block)
  syms = [gensym() for _ ∈ 1:nargs]
  configarg = (inline, u₁, u₂, v, true, threads, warncheckarg)
  unroll_param_tup = Expr(:call, lv(:avx_config_val), :(Val{$configarg}()), staticexpr(0))
  for n ∈ 1:nargs
    ciₙ = ci[n]
    ciₙargs = ciₙ.args
    f = first(ciₙargs)
    if ciₙ.head === :(=)
      push!(ex.args, Expr(:(=), f, syms[((ciₙargs[2])::Core.SSAValue).id]))
    elseif isglobalref(f, Base, :materialize!)
      add_ci_call!(ex, lv(:vmaterialize!), ciₙargs, syms, n, unroll_param_tup, mod)
    elseif isglobalref(f, Base, :materialize)
      add_ci_call!(ex, lv(:vmaterialize), ciₙargs, syms, n, unroll_param_tup, mod)
    else
      add_ci_call!(ex, f, ciₙargs, syms, n)
    end
  end
  esc(ex)
end

function LoopSet(q::Expr, mod::Symbol = :Main)
  ls = LoopSet(mod)
  check_inputs!(q, ls.prepreamble)
  contract_pass!(q)
  copyto!(ls, q)
  resize!(ls.loop_order, num_loops(ls))
  ls
end
LoopSet(q::Expr, m::Module) = LoopSet(macroexpand(m, q)::Expr, Symbol(m))

function loopset(q::Expr) # for interactive use only
  ls = LoopSet(q)
  set_hw!(ls)
  ls
end

function check_macro_kwarg(
  arg,
  inline::Bool,
  check_empty::Bool,
  u₁::Int8,
  u₂::Int8,
  v::Int8,
  threads::Int,
  warncheckarg::Int,
)
  ((arg.head === :(=)) && (length(arg.args) == 2)) ||
    throw(ArgumentError("macro kwarg should be of the form `argname = value`."))
  kw = (arg.args[1])::Symbol
  value = (arg.args[2])
  if kw === :inline
    inline = value::Bool
  elseif kw === :unroll
    if value isa Integer
      u₁ = convert(Int8, value)::Int8
    elseif Meta.isexpr(value, :tuple, 2)
      u₁ = convert(Int8, value.args[1])::Int8
      u₂ = convert(Int8, value.args[2])::Int8
    else
      throw(ArgumentError("Don't know how to process argument in `unroll=$value`."))
    end
  elseif kw === :vectorize
    v = convert(Int8, value)
  elseif kw === :check_empty
    check_empty = value::Bool
  elseif kw === :thread
    if value isa Bool
      threads = Core.ifelse(value::Bool, -1, 1)
    elseif value isa Integer
      threads = max(1, convert(Int, value)::Int)
    else
      throw(ArgumentError("Don't know how to process argument in `thread=$value`."))
    end
  elseif kw === :warn_check_args
    warncheckarg = convert(Int, value)::Int
  else
    throw(
      ArgumentError(
        "Received unrecognized keyword argument $kw. Recognized arguments include:\n`inline`, `unroll`, `check_empty`, and `thread`.",
      ),
    )
  end
  inline, check_empty, u₁, u₂, v, threads, warncheckarg
end
function process_args(
  args;
  inline::Bool = false,
  check_empty::Bool = false,
  u₁::Int8 = zero(Int8),
  u₂::Int8 = zero(Int8),
  v::Int8 = zero(Int8),
  threads::Int = 1,
  warncheckarg::Int = 1,
)
  for arg ∈ args
    inline, check_empty, u₁, u₂, v, threads, warncheckarg =
      check_macro_kwarg(arg, inline, check_empty, u₁, u₂, v, threads, warncheckarg)
  end
  inline, check_empty, u₁, u₂, v, threads, warncheckarg
end
# check if the body of loop is a block, if not convert it to a block issue#395
# and check if the range of loop is an enumerate, if it is replace it, issue#393
function check_inputs!(q, prepreamble)
  if Meta.isexpr(q, :for)
    if !Meta.isexpr(q.args[2], :block)
      q.args[2] = Expr(:block, q.args[2])
      replace_enumerate!(q, prepreamble) # must after warp block
    else # maybe inner loops in block
      replace_enumerate!(q, prepreamble)
      for arg in q.args[2].args
        check_inputs!(arg, prepreamble) # check recursively for inner loop
      end
    end
  end
  return q
end
function replace_enumerate!(q, prepreamble)
  looprange = q.args[1]
  if Meta.isexpr(looprange, :block)
    for i = 1:length(looprange.args)
      replace_single_enumerate!(q, prepreamble, i)
    end
  else
    replace_single_enumerate!(q, prepreamble)
  end
  return q
end
function replace_single_enumerate!(q, prepreamble, i = nothing)
  if isnothing(i) # not nest loop
    looprange, body = q.args[1], q.args[2]
  else # nest loop
    looprange, body = q.args[1].args[i], q.args[2]
  end
  @assert Meta.isexpr(looprange, :(=), 2)
  itersyms, r = looprange.args
  if Meta.isexpr(r, :call, 2) && r.args[1] == :enumerate
    _iter = r.args[2]
    if _iter isa Symbol
      iter = _iter
    else # name complex expr
      iter = gensym(:iter)
      push!(prepreamble.args, :($iter = $_iter))
    end
    if Meta.isexpr(itersyms, :tuple, 2)
      indsym, varsym = itersyms.args[1]::Symbol, itersyms.args[2]::Symbol
      _replace_looprange!(q, i, indsym, iter)
      pushfirst!(body.args, :($varsym = $iter[$indsym+firstindex($iter)-1]))
    elseif Meta.isexpr(itersyms, :tuple, 1) # like `for (i,) in enumerate(...)`
      indsym = itersyms.args[1]::Symbol
      _replace_looprange!(q, i, indsym, iter)
    elseif itersyms isa Symbol # if itersyms are not unbox in loop range
      throw(ArgumentError("`for $itersyms in enumerate($r)` is not supported,
        please use `for ($(itersyms)_i, $(itersyms)_v) in enumerate($r)` instead."))
    else
      throw(ArgumentError("Don't know how to handle expression `$itersyms`."))
    end
  end
  return q
end
_replace_looprange!(q, ::Nothing, indsym, iter) =
  q.args[1] = :($indsym = Base.OneTo(length($iter)))
_replace_looprange!(q, i::Int, indsym, iter) =
  q.args[1].args[i] = :($indsym = Base.OneTo(length($iter)))

function turbo_macro(mod, src, q, args...)
  q = macroexpand(mod, q)
  if q.head === :for
    ls = LoopSet(q, mod)
    inline, check_empty, u₁, u₂, v, threads, warncheckarg = process_args(args)
    esc(setup_call(ls, q, src, inline, check_empty, u₁, u₂, v, threads, warncheckarg))
  else
    inline, check_empty, u₁, u₂, v, threads, warncheckarg =
      process_args(args, inline = true)
    substitute_broadcast(q, Symbol(mod), inline, u₁, u₂, v, threads, warncheckarg)
  end
end
"""
    @turbo

Annotate a `for` loop, or a set of nested `for` loops whose bounds are constant across iterations, to optimize the computation. For example:

    function AmulB!(C, A, B)
        @turbo for m ∈ indices((A,C), 1), n ∈ indices((B,C), 2) # indices((A,C),1) == axes(A,1) == axes(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ indices((A,B), (2,1)) # indices((A,B), (2,1)) == axes(A,2) == axes(B,1)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end

The macro models the set of nested loops, and chooses an ordering of the three loops to minimize predicted computation time.

Current limitations:

1. It assumes that loop iterations are independent.
2. It does not perform bounds checks.
3. It assumes that each loop iterates at least once. (Use `@turbo check_empty=true` to lift this assumption.)
4. That there is only one loop at each level of the nest.

It may also apply to broadcasts:

```jldoctest
julia> using LoopVectorization

julia> a = rand(100);

julia> b = @turbo exp.(2 .* a);

julia> c = similar(b);

julia> @turbo @. c = exp(2a);

julia> b ≈ c
true
```

# Extended help

Advanced users can customize the implementation of the `@turbo`-annotated block
using keyword arguments:

```julia
@turbo inline=false unroll=2 thread=4 body
```

where `body` is the code of the block (e.g., `for ... end`).

`thread` is either a Boolean, or an integer.
The integer's value indicates the number of threads to use.
It is clamped to be between `1` and `min(Threads.nthreads(),LoopVectorization.num_cores())`.
`false` is equivalent to `1`, and `true` is equivalent to `min(Threads.nthreads(),LoopVectorization.num_cores())`.

`inline` is a Boolean. When `true`, `body` will be directly inlined
into the function (via a forced-inlining call to `_turbo_!`).
When `false`, it wont force inlining of the call to `_turbo_!` instead, letting Julia's own inlining engine
determine whether the call to `_turbo_!` should be inlined. (Typically, it won't.)
Sometimes not inlining can lead to substantially worse code generation, and >40% regressions, even in very
large problems (2-d convolutions are a case where this has been observed).
One can find some circumstances where `inline=true` is faster, and other circumstances
where `inline=false` is faster, so the best setting may require experimentation. By default, the macro
tries to guess. Currently the algorithm is simple: roughly, if there are more than two dynamically sized loops
or and no convolutions, it will probably not force inlining. Otherwise, it probably will.

`check_empty` (default is `false`) determines whether or not it will check if any of the iterators are empty.
If false, you must ensure yourself that they are not empty, else the behavior of the loop is undefined and
(like with `@inbounds`) segmentation faults are likely.

`unroll` is an integer that specifies the loop unrolling factor, or a
tuple `(u₁, u₂) = (4, 2)` signaling that the generated code should unroll more than
one loop. `u₁` is the unrolling factor for the first unrolled loop and `u₂` for the next (if present),
but it applies to the loop ordering and unrolling that will be chosen by LoopVectorization,
*not* the order in `body`.
`uᵢ=0` (the default) indicates that LoopVectorization should pick its own value,
and `uᵢ=-1` disables unrolling for the correspond loop.

The `@turbo` macro also checks the array arguments using `LoopVectorization.check_args` to try and determine
if they are compatible with the macro. If `check_args` returns false, a fall back loop annotated with `@inbounds`
and `@fastmath` is generated. Note that `VectorizationBase` provides functions such as `vadd` and `vmul` that will
ignore `@fastmath`, preserving IEEE semantics both within `@turbo` and `@fastmath`.
`check_args` currently returns false for some wrapper types like `LinearAlgebra.UpperTriangular`, requiring you to
use their `parent`. Triangular loops aren't yet supported.

Setting the keyword argument `warn_check_args=true` (e.g. `@turbo warn_check_args=true for ...`) in a loop or
broadcast statement will cause it to warn once if `LoopVectorization.check_args` fails and the fallback
loop is executed instead of the LoopVectorization-optimized loop.
Setting it to an integer > 0 will warn that many times, while setting it to a negative integer will warn
an unlimited amount of times. The default is `warn_check_args = 0`.
"""
macro turbo(args...)
  turbo_macro(__module__, __source__, last(args), Base.front(args)...)
end
"""
    @tturbo

Equivalent to `@turbo`, except it adds `thread=true` as the first keyword argument.
Note that later arguments take precendence.

Meant for convenience, as `@tturbo` is shorter than `@turbo thread=true`.
"""
macro tturbo(args...)
  turbo_macro(__module__, __source__, last(args), :(thread = true), Base.front(args)...)
end

function def_outer_reduct_types!(ls::LoopSet)
  for or ∈ ls.outer_reductions
    op = operations(ls)[or]
    pushpreamble!(ls, Expr(:(=), outer_reduct_init_typename(op), typeof_expr(op)))
  end
end
"""
    @_turbo

This macro mostly exists for debugging/testing purposes. It does not support many of the use cases of [`@turbo`](@ref).
It emits loops directly, rather than punting to an `@generated` function, meaning it doesn't have access to type
information when generating code or analyzing the loops, often leading to bad performance.

This macro accepts the `inline` and `unroll` keyword arguments like `@turbo`, but ignores the `check_empty` argument.
"""
macro _turbo(q)
  q = macroexpand(__module__, q)
  ls = LoopSet(q, __module__)
  set_hw!(ls)
  def_outer_reduct_types!(ls)
  esc(Expr(:block, ls.prepreamble, lower_and_split_loops(ls, -1)))
end
macro _turbo(arg, q)
  @assert q.head === :for
  q = macroexpand(__module__, q)
  inline, check_empty, u₁, u₂, v =
    check_macro_kwarg(arg, false, false, zero(Int8), zero(Int8), zero(Int8), 1, 0)
  ls = LoopSet(q, __module__)
  set_hw!(ls)
  def_outer_reduct_types!(ls)
  esc(Expr(:block, ls.prepreamble, lower(ls, u₁ % Int, u₂ % Int, v % Int, -1)))
end

"""
    @turbo_debug

Returns a `LoopSet` object instead of evaluating the loops. Useful for debugging and introspection.
"""
macro turbo_debug(q)
  q = macroexpand(__module__, q)
  ls = LoopSet(q, __module__)
  esc(LoopVectorization.setup_call_debug(ls))
end

# define aliases
const var"@avx" = var"@turbo"
const var"@avxt" = var"@tturbo"
const var"@avx_debug" = var"@turbo_debug"


text: include("shellchron.jl")


text:
#=
BellmanFord:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays

#A is the input matrix, s is the source node, and n is the number of nodes, Output: Shortest path from source node
function bellmanford(A, s, n)

    #vector init +inf
    d = GBVector{Float64}(n)
        for i = 1:n
            d[i] = Inf
        end

    d[s]=0.0
    for i = 1:n
           d = mul(d, A, Semirings.MIN_PLUS, mask=d, desc=Descriptors.S)
            print(d)
    end

    return d

end


text:
function drug_transport_organ!(u, p, t)
    F, Cart, P, V = p
    F*(Cart - u/(P * V))
end

"""
Physiologically based pharmacokinetic model (PBPM)

https://en.wikipedia.org/wiki/Physiologically_based_pharmacokinetic_modelling
"""
function physio_pharma_drug_transport!(du, u, p, t)
    Cart, F, P, V = p
    for i in 1:length(du)
        p_i = (Cart, F[i], P[i], V[i])
        du[i] = drug_transport_organ!(u[i], p_i, t)
    end
    nothing
end


text:
#using NeuralQuantum
using Test

@testset "NeuralQuantum" begin

    println("Testing Operators...")
    @testset "Operators" begin
        include("Operators/operators.jl")
        include("Operators/ising.jl")
    end

    println("Testing machines...")
    @testset "Machines" begin
        include("Machines/test_cached.jl")
        include("Machines/test_grad.jl")
        include("Machines/test_batched.jl")
        include("Machines/test_ndmcomplex.jl")
    end

    println("Testing samplers...")
    @testset "Samplers" begin
        include("Samplers/test_samplers.jl")
    end

    #=
    println("Testing algorithm caches...")
    @testset "Algorithm" begin
        include("Algorithms_caches/Caches_test.jl")
    end

    println("Testing Problems and algorithms...")
    @testset "Problem & algs" begin
        include("Problems/problem_algs.jl")
        include("Problems/problem_operators.jl")
        include("Problems/observables.jl")
        include("Problems/hamiltonian.jl")
    end
    =#
end


text:
using CausalInference
using LightGraphs
using Test
using Random
using Statistics
include("plotdag.jl")

# After examples discussed in chapter 2 of
# Pearl, Judea. Causality. Cambridge University Press, 2009.

Random.seed!(123)
N = 1000
p = 0.01
x = randn(N)
v = x + randn(N)*0.25
w = x + randn(N)*0.25
z = v + w + randn(N)*0.25
s = z + randn(N)*0.25

X = [x v w z s]
C = cor(X)
df = (x=x, v=v, w=w, z=z, s=s)

println("Running Gaussian tests")
@time est_g = pcalg(df, p, gausscitest)

variables = [String(k) for k in keys(df)]
tp = plot_dag(est_g, variables)
save(PDF("estdag"), tp)

# This should be the true DAG

g = DiGraph(5)
d = nv(g)
for (i,j) in [(1,2), (1,3), (2,4), (3,4), (4,5)]
   add_edge!(g, i, j)
end

tp = plot_dag(g, variables)
save(PDF("truedag"), tp)

# Use the oracle to tell about a graph what can be learned from observations
dg = pcalg(d, dseporacle, g)

tp = plot_dag(dg, variables)
save(PDF("equivalencedag"), tp)

@testset "pctest" begin
    @test collect(LightGraphs.edges(est_g)) == collect(LightGraphs.edges(dg))
end


text:
module mmmNAFEMS_R0031_3m
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import Statistics: mean
function test()
  # println("""
  # NAFEMS publication R0031/3 Composite plate test.
  # Simply supported on all four edges.  Uniform transverse  loading.
  # The modeled part is one quarter of the full plate here.
  # """)

  # This is a test recommended by the National Agency for Finite Element Methods
  # and Standards (U.K.): Test R0031/3 from NAFEMS publication R0031, “Composites
  # Benchmarks,” February 1995.
  t0 = time()
  # Skin (face) material parameters
  E1s = 1.0e7*phun("psi")
  E2s = 0.4e7*phun("psi")
  E3s = 0.4e7*phun("psi")
  nu12s = 0.3
  nu13s = 0.3
  nu23s = 0.3
  G12s = 0.1875e7*phun("psi")
  G13s = 0.1875e7*phun("psi")
  G23s = 0.1875e7*phun("psi")
  # Core material parameters
  E1c = 10.0*phun("psi")
  E2c = 10.0*phun("psi")
  E3c = 10e4.*phun("psi")
  nu12c = 0.
  nu13c = 0.
  nu23c = 0.
  G12c = 10.0*phun("psi")
  G13c = 3.0e4*phun("psi")
  G23c = 1.2e4*phun("psi")
  L = 10.0*phun("in") # side of the square plate
  nL = 8 # number of elements along the side of the plate
  tolerance = 0.0001*phun("in")
  xs = collect(linearspace(0.0, L/2, nL+1))
  ys = collect(linearspace(0.0, L/2, nL+1))
  ts = [0.028; 0.75; 0.028]*phun("in")
  nts = [2; 3;  2; ] # number of elements through the thickness
  tmag = 100*phun("psi")

  # Generate mesh
  fens,fes = H8layeredplatex(xs, ys, ts, nts)
  fens,fes = H8toH20(fens,fes)

  MR = DeforModelRed3D
  skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    0.0, 0.0, 0.0)
  corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    0.0, 0.0, 0.0)

  gr = GaussRule(3, 3)

  rl1 = selectelem(fens, fes, label=1)
  skinbot = FDataDict("femm"=>FEMMDeforLinear(MR,
      IntegDomain(subset(fes, rl1), gr), skinmaterial))

  rl3 = selectelem(fens, fes, label=3)
  skintop = FDataDict("femm"=>FEMMDeforLinear(MR,
      IntegDomain(subset(fes, rl3), gr), skinmaterial))

  rl2 = selectelem(fens, fes, label=2)
  core = FDataDict("femm"=>FEMMDeforLinear(MR,
      IntegDomain(subset(fes, rl2), gr), corematerial))

  lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
  lxL2 = selectnode(fens, box=[L/2 L/2 -Inf Inf -Inf Inf], inflate=tolerance)
  ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
  lyL2 = selectnode(fens, box=[-Inf Inf L/2 L/2 -Inf Inf], inflate=tolerance)

  ex0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
  exL2 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lxL2 )
  ey0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>ly0 )
  eyL2 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lyL2 )

  bfes = meshboundary(fes)
  ttopl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 1.0])
  Trac = FDataDict("traction_vector"=>[0.0; 0.0; -tmag],
      "femm"=>FEMMBase(IntegDomain(subset(bfes, ttopl), GaussRule(2, 3))))

  modeldata = FDataDict("fens"=>fens,
   "regions"=>[skinbot, core, skintop], "essential_bcs"=>[ex0, exL2, ey0, eyL2],
   "traction_bcs"=> [Trac]
   )
  modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

  u = modeldata["u"]
  geom = modeldata["geom"]
  lcenter = selectnode(fens, box=[L/2 L/2  L/2 L/2 -Inf Inf], inflate=tolerance)
  cdis = mean(u.values[lcenter, 3])/phun("in")
  # println("Center node displacements $(cdis) [in]; NAFEMS-R0031-3 lists –0.123	[in]")
  # println("")
  @test abs(cdis - (-0.13634800328800462)) < 1.0e-5

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

end
end
using .mmmNAFEMS_R0031_3m
mmmNAFEMS_R0031_3m.test()

module mmtwistedmsh8mmm
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
import LinearAlgebra: norm, cholesky, cross
using Test
import Statistics: mean
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 is
  #   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 = 7;
  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])

  # 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")
  @test abs(theutip[dir]-uex)/uex < 0.0012

  # Write out mesh with displacements
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8")
  modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :Cauchy, "component"=> :xy)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with von Mises stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  vm  = modeldata["postprocessing"]["exported"][1]["field"]
  # println("extremes of vm, nodal: $([minimum(vm.values),   maximum(vm.values)])")
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
  @test norm([minimum(vm.values),   maximum(vm.values)]-[4.78774, 522.126]) < 0.01

  # Write out mesh with von Mises stresses, elementwise
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  vm  = modeldata["postprocessing"]["exported"][1]["field"]
  # println("extremes of vm, elemental: $([minimum(vm.values),   maximum(vm.values)])")
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
  @test norm([minimum(vm.values),   maximum(vm.values)]-[1.85882, 522.126]) < 0.01

  # Write out mesh with von Mises stresses, elementwise
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

end
end
using .mmtwistedmsh8mmm
mmtwistedmsh8mmm.test()

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

  # println("""
  # % Vibration modes of unit cube  of almost incompressible material.
  # % Mean-strain hexahedron.
  # % Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
  # % tetrahedral. International Journal for Numerical Methods in
  # % Engineering 67: 841-867.""")
  t0 = time()

  E = 1*phun("PA");
  nu = 0.499;
  rho= 1*phun("KG/M^3");
  a=1*phun("M"); b=a; h= a;
  n1=8 # How many element edges per side?
  na= n1; nb= n1; nh =n1;
  neigvs=20                   # how many eigenvalues
  omega_shift=(0.1*2*pi)^2;

  fens,fes = H8block(a,b,h, na,nb,nh)

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

  # Make model data
  modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1],
    "omega_shift"=>omega_shift, "neigvs"=>neigvs)

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

  fs = modeldata["omega"]/(2*pi)
  # println("Eigenvalues: $fs [Hz]")
  @test norm(fs-[1.92866e-7, 2.07497e-7, 2.16105e-7, 2.31656e-7, 2.35711e-7, 2.53067e-7, 0.266016, 0.266016, 0.364001, 0.364001, 0.364001, 0.366888, 0.366888, 0.366888, 0.415044, 0.415044, 0.41703, 0.467364, 0.467364, 0.467364]) < 0.0001

  modeldata["postprocessing"] = FDataDict("file"=>"unit_cube_mode",
    "mode"=>10)
  modeldata=FinEtoolsDeforLinear.AlgoDeforLinearModule.exportmode(modeldata)
  # @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  true

end
end
using .mmunitmmccubemmvibrationmmms
mmunitmmccubemmvibrationmmms.test()

module mmtwistedbeamisomm
using FinEtools
using FinEtoolsDeforLinear
using Test
using FinEtoolsDeforLinear.AlgoDeforLinearModule
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()
  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  = H20block(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"=>FEMMDeforLinear(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 abs(theutip[dir]/uex*100-99.85504856450584) < 1.0e-6

  # Write out mesh with displacements
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam")
  modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :Cauchy, "component"=> :xy)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with von Mises stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # 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)] - [6.94796, 451.904]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

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

  # 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)] - [0.493918, 459.106]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # 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)] - [-459.106, -0.493918]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with pressure, 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)] - [-160.396, 160.396]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

end
end
using .mmtwistedbeamisomm
mmtwistedbeamisomm.test()

module mmtwistedbeamoorthomm
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
using FinEtoolsDeforLinear.AlgoDeforLinearModule
import Statistics: mean
function test()
  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  = H20block(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 = MatDeforElastOrtho(MR, E, nu)
  region1 = FDataDict("femm"=>FEMMDeforLinear(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 abs(theutip[dir]/uex*100-99.85504856450584) < 1.0e-6

  # Write out mesh with displacements
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam")
  modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :Cauchy, "component"=> :xy)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with von Mises stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # 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)] - [6.94796, 451.904]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

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

  # 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)] - [0.493918, 459.106]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # 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)] - [-459.106, -0.493918]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with pressure, 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)] - [-160.396, 160.396]) < 1.e-3
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

end
end
using .mmtwistedbeamoorthomm
mmtwistedbeamoorthomm.test()

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

  # println("""
  # Vibration modes of unit cube  of almost incompressible material.
  #
  # This example EXPORTS the model to Abaqus.
  #
  # Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
  # tetrahedral. International Journal for Numerical Methods in
  # Engineering 67: 841-867.
  # """)
  t0 = time()

  E = 1*phun("PA");
  nu = 0.499;
  rho = 1*phun("KG/M^3");
  a = 1*phun("M"); b = a; h =  a;
  n1 = 5;# How many element edges per side?
  na =  n1; nb =  n1; nh  = n1;
  neigvs = 20                   # how many eigenvalues
  OmegaShift = (0.01*2*pi)^2;

  MR = DeforModelRed3D
  fens,fes  = H20block(a,b,h, na,nb,nh)

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

  numberdofs!(u)

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

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

  K =stiffness(femm, geom, u)
  femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
  M =mass(femm, geom, u)
  d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
  d = broadcast(-, d, OmegaShift);
  fs = real(sqrt.(complex(d)))/(2*pi)
  # println("Eigenvalues: $fs [Hz]")
  @test norm(fs - [2.73674e-7, 3.00469e-7, 3.14245e-7, 3.19946e-7, 3.42634e-7, 3.56347e-7, 0.262723, 0.262723, 0.357791, 0.357791, 0.357791, 0.36088, 0.36088, 0.36088, 0.408199, 0.408397, 0.408397, 0.461756, 0.461756, 0.461756]) < 1.0e-3

  mode = 7
  scattersysvec!(u, v[:,mode])
  File =  "unit_cube_modes.vtk"
  vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
  # @async run(`"paraview.exe" $File`)
  try rm(File); catch end

  AE = AbaqusExporter("unit_cube_modes_h20");
  # AE.ios = STDOUT;
  HEADING(AE, "Vibration modes of unit cube  of almost incompressible material.");
  COMMENT(AE, "The  first six frequencies are rigid body modes.");
  COMMENT(AE, "The  first nonzero frequency (7) should be around 0.26 Hz");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  COMMENT(AE, "The hybrid form of the serendipity hexahedron is chosen because");
  COMMENT(AE, "the material is  nearly incompressible.");
  ELEMENT(AE, "c3d20rh", "AllElements", 1, connasarray(fes))
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  DENSITY(AE, rho)
  STEP_FREQUENCY(AE, neigvs)
  END_STEP(AE)
  close(AE)
  nlines = 0
  open("unit_cube_modes_h20.inp") do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 1409
  rm("unit_cube_modes_h20.inp")
end
end
using .muunit_cube_modes_exportmmm
muunit_cube_modes_exportmmm.test()

module mmpipemmPSmorthom
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross, dot
import Statistics: mean

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

function test()

## Thick pipe with internal pressure: plane strain
#

##
# Link to the  <matlab:edit('pub_thick_pipe_ps') m-file>.

## Description
##
# This is a simple modification of the full three-dimensional simulation of
# the tutorial pub_thick_pipe takes advantage of the plane-strain 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.

##
#
# <html>
# <table border=0><tr><td>
# <img src="../docs/pub_thick_pipe_ps.png" width = "30#">
# </td></tr>
# <tr><td>Figure 1. Definition of the geometry of the internally pressurized thick pipe</td></tr>
# </table>
# </html>

##
# 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.

## Solution
#

##
# 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).

nc=3; nt=3;

##
# 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 = DeforModelRed2DStrain

# Create the mesh and initialize the geometry.  First we are going
# to construct the block of elements with the first coordinate
# corresponding to the angle, and the second
# coordinate is the thickness in the radial direction.
anglrrange = 90.0/180*pi;
fens,fes =  Q8block(anglrrange, b-a, nc, nt);

# Extract the boundary  and mark the finite elements on the
# interior surface.
bdryfes = meshboundary(fes);
bcl = selectelem(fens, bdryfes, box=[-Inf,Inf,0.,0.], inflate=tolerance);
internal_fenids = connectednodes(subset(bdryfes,bcl));
# Now  shape the block  into  the actual wedge piece of the pipe.
ayr=fens.xyz;
for i=1:count(fens)
    angl=ayr[i,1]; r=a+ayr[i,2];
    fens.xyz[i,:] = [r*sin(angl),(r*cos(angl))];
end

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

# The symmetry boundary condition  is specified by selecting nodes
# on the plane x=0.
l1 = selectnode(fens; box=[0.0 0.0 -Inf Inf], inflate = tolerance)
setebc!(u, l1, true, 1, 0.0)
# The second symmetry boundary condition is specified by selecting
# nodes on the plane y=0.
l1 = selectnode(fens; box=[-Inf Inf 0.0 0.0], 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)))
function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
  copyto!(forceout, XYZ/norm(XYZ)*press)
  return forceout
end
fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
F2 = distribloads(el1femm, geom, u, fi, 2);

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

femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2)), 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,:];
ur = zeros(FFlt,length(internal_fenids));

for j = 1:length(internal_fenids)
    n = fens.xyz[internal_fenids[j],:];
    n = n'/norm(n);# normal to the cylindrical internal surface
    ur[j] = dot(vec(uv[j,:]),vec(n));
end

# Report the  relative displacement on the internal surface:
# println("(Approximate/true displacement) at the internal surface: $( mean(ur)/urex*100  ) %")
@test abs(mean(ur)/urex*100 - 100) < 0.1

# 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)])
try rm(File); catch end

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

function inspector(idat::MyIData, elnum, conn, xe,  out,  xq)
  function outputRm(c)
    theNormal=c;
    r=norm(theNormal);# distance from the axis of symmetry
    theNormal =theNormal/r;# compute the unit normal vector
    e1p=[theNormal';0.];# local cylind. coordinate  basis vectors
    e3p=[0.,0.,1.]';# this one points along the axis of the cylinder
    e2p=cross(vec(e3p),vec(e1p));# this one is along the hoop direction
    R= [vec(e1p) vec(e2p) vec(e3p)];# transformation matrix for the stress
    return R
  end
  Rm=outputRm(xq)
  tm=zeros(FFlt,3,3)
  stressvtot!(MR, tm, out);# stress in global XYZ
  tpm = Rm'*tm*Rm;#  stress matrix in cylindrical coordinates
  sp=zeros(FFlt,6)
  stressttov!(MR, sp, tpm);# stress vector in cylindr. coord.
  push!(idat.r,norm(xq))
  push!(idat.s,sp[idat.c])
  return idat
end

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

@test norm(idat.s - [-7.44858e5, -3.55143e5, -7.44858e5, -3.55143e5, -2.19961e5,
-1.08612e5, -2.19961e5, -1.08612e5, -58910.9, -12517.6, -58910.9, -12517.6,
-7.44858e5, -3.55143e5, -7.44858e5, -3.55143e5, -2.19961e5, -1.08612e5,
-2.19961e5, -1.08612e5, -58910.9, -12517.6, -58910.9, -12517.6, -7.44858e5,
-3.55143e5, -7.44858e5, -3.55143e5, -2.19961e5, -1.08612e5, -2.19961e5,
 -1.08612e5, -58910.9, -12517.6, -58910.9, -12517.6])/1.0e5 < 1.e-3
# using Plots
# plotly()
#
# # Plot the analytical solution.
# r = linearspace(a,b,100);
# plot(r, radial_stress(r))
# # Plot the computed  integration-point data
# scatter!(idat.r, idat.s, m=:circle, color=:red)
# gui()
end
end
using .mmpipemmPSmorthom
mmpipemmPSmorthom.test()


text:
using Test, Random, LSHFunctions

include(joinpath("..", "utils.jl"))

#==================
Tests
==================#
@testset "LpHash tests" begin
    Random.seed!(RANDOM_SEED)

    @testset "Can construct an ℓ^p distance hash function" begin
        # Construct a hash for L^1 distance
        L1_hash = L1Hash(5; scale = 2)
        @test n_hashes(L1_hash) == 5
        @test L1_hash.scale == 2
        @test L1_hash.power == 1
        @test similarity(L1_hash) == ℓ1
        @test hashtype(L1_hash) == Int32

        # Construct a hash for L^2 distance
        L2_hash = L2Hash(12; scale = 3.4)
        @test n_hashes(L2_hash) == 12
        @test L2_hash.scale == Float32(3.4)
        @test L2_hash.power == 2
        @test similarity(L2_hash) == ℓ2

        # Construct a hash using a specified dtype
        @test isa(L1Hash(1; dtype=Float32), LSHFunctions.LpHash{Float32})
        @test isa(L2Hash(1; dtype=Float64), LSHFunctions.LpHash{Float64})
    end

    @testset "Hashes are correctly computed" begin
        n_hashes = 8
        scale = 2

        hashfn = L2Hash(n_hashes; scale=scale)

        # Test on a single input
        x = randn(8)
        hashes = hashfn(x)
        manual_hashes = floor.(Int32, hashfn.coeff * x ./ scale .+ hashfn.shift)

        @test isa(hashes, Vector{Int32})
        @test hashes == manual_hashes

        # Test on many inputs, simultaneously
        x = randn(8, 128)
        hashes = hashfn(x)
        manual_hashes = floor.(Int32, hashfn.coeff * x ./ scale .+ hashfn.shift)

        @test isa(hashes, Matrix{Int32})
        @test hashes == manual_hashes
    end

    @testset "Hashes have the correct dtype" begin
        hashfn = L1Hash(5)

        # Test 1: Vector{Float64} -> Vector{Int32}
        hashes = hashfn(randn(5))
        @test eltype(hashes) == hashtype(hashfn)
        @test isa(hashes, Vector{Int32})

        # Test 2: Matrix{Float64} -> Matrix{Int32}
        hashes = hashfn(randn(5, 10))
        @test isa(hashes, Matrix{Int32})
    end

    @testset "Nearby points experience more frequent collisions" begin
        hashfn = L2Hash(1024; dtype=Float64, scale=4)

        x1 = randn(128)
        x2 = x1 + 0.05 * randn(length(x1))
        x3 = x1 + 0.50 * randn(length(x1))
        x4 = x1 + 2.00 * randn(length(x1))

        h1, h2, h3, h4 = hashfn(x1), hashfn(x2), hashfn(x3), hashfn(x4)
        @test sum(h1 .== h4) < sum(h1 .== h3) < sum(h1 .== h2)
    end

    @testset "Hash collision frequency matches probability" begin
        hashfn = L2Hash(1024; scale = 4)

        # Dry run
        @test test_collision_probability(hashfn, 0.05)

        # Full test
        @test all(test_collision_probability(hashfn, 0.05) for ii = 1:128)
    end

    @testset "Hash inputs of different sizes" begin
        n_hashes = 10
        for hashfn_type in (:L1Hash, :L2Hash)
            hashfn = eval(hashfn_type)(n_hashes)
            @test size(hashfn.coeff) == (n_hashes, 0)

            hashfn(rand(5))
            @test size(hashfn.coeff) == (n_hashes, 5)

            hashfn(rand(20))
            @test size(hashfn.coeff) == (n_hashes, 20)

            hashfn(rand(10))
            @test size(hashfn.coeff) == (n_hashes, 20)
        end
    end

    @testset "Hash inputs of different sizes with resize_pow2 = true" begin
        n_hashes = 20
        for hashfn_type in (:L1Hash, :L2Hash)
            hashfn = eval(hashfn_type)(n_hashes; resize_pow2=true)
            @test size(hashfn.coeff) == (n_hashes, 0)

            hashfn(rand(4))
            @test size(hashfn.coeff) == (n_hashes, 4)

            hashfn(rand(5))
            @test size(hashfn.coeff) == (n_hashes, 8)

            hashfn(rand(10))
            @test size(hashfn.coeff) == (n_hashes, 16)

            hashfn(rand(8))
            @test size(hashfn.coeff) == (n_hashes, 16)
        end
    end

    @testset "collision_probability works correctly" begin
        hashfn = L1Hash()

        # collision_probability should be 1 for two inputs of distance zero
        x = rand(4)
        @test collision_probability(hashfn, x, x) ≈ 1.0

        # collision_probability with n_hashes=N should be the same as
        # collision_probability with n_hashes=1, raised to the power N
        y = rand(4)
        @test collision_probability(hashfn, x, y; n_hashes=10) ≈
              collision_probability(hashfn, x, y; n_hashes=1)^10
    end
end


text:
@testset "contour_tracing" begin

    for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
        for i = 1:20
            test_image = eval(Symbol("test_image_$(i)"))

            # Call contour tracing algorithm and compare the number of components
            labels = ImageComponentAnalysis.label_components(ContourTracing(), test_image)
            num_components = maximum(labels)

            @test num_components == test_image_results[i]
        end
    end

end


text:
using AlgebraicDecisionDiagrams
using Test

# Boolean Decision Diagrams
include("bdd.jl")
# Algebraic Decision Diagrams with real constants
include("add.jl")
# Parametrized Algebraic Decision Diagrams (ADDs with multilinear expressions as constants)
include("padd.jl")


text:
#__precompile__(true)

module UnsteadyFlowSolvers

import Dierckx: Spline1D, derivative, evaluate

import ForwardDiff

import DelimitedFiles

import Serialization

import NLsolve: nlsolve, not_in_place

import Statistics: mean

import PyPlot: plot, scatter, figure, xlabel, ylabel, xlim, ylim,
    xticks, yticks, subplot, subplot2grid, legend, axis, savefig,
    close, tight_layout

import Plots: @layout

import LaTeXStrings: @L_str

#For use in development and debugging
import Revise

import Printf: @printf

export
    # kinematics types and funtions
    MotionDef,
    KinemPar,
    KinemDef,
    EldUpDef,
    EldUptstartDef,
    ConstDef,
    EldRampReturnDef,
    EldUpIntDef,
    EldUpInttstartDef,
    SinDef,
    CosDef,

    # 2D low-order solver types
    TwoDSurf,
    TwoDOFPar,
    KinemPar2DOF,
    TwoDVort,
    TwoDVVort,
    TwoDFlowField,
    TwoDVFlowField,
    KelvinCondition,
    KelvinKutta,

    # vortex count control utility
    delVortDef,
    delNone,
    delSpalart,

    # utility functions
    simpleTrapz,
    camber_calc,
    find_tstep,
    simpleInterp,
    cleanWrite,

    #2D low-order solver methods
    lautat,
    ldvm,
    ldvmLin,
    ldvm2DOF,

    # 2D plot output functions
    makeForcePlots2D,
    makeVortPlots2D

### source files

# kinematic types
include("kinem.jl")

# utility functions
include("utils.jl")

# vortex count control utility
include("delVort.jl")

# low-order 2D solvers
include("lowOrder2D/typedefs.jl")            # type definitions
include("lowOrder2D/calcs.jl")               # calculation functions
include("lowOrder2D/solvers.jl")             # solver methods
include("lowOrder2D/postprocess.jl")         # postprocessing functions

# 2D plotting functions
include("plots/plots2D.jl")

end


text:
sm = @stochastic_model begin
    @stage 1 begin
        @decision(model, x¹ <= 2)
        @decision(model, w¹)
        @decision(model, y¹)
        @objective(model, Min, 100*x¹ + 3*w¹ + 0.5*y¹)
        @constraint(model, x¹ + w¹ - y¹ == 1)
    end
    @stage 2 begin
        @uncertain d
        @variable(model, x² <= 2)
        @variable(model, w²)
        @decision(model, y²)
        @objective(model, Min, x² + 3*w² + 0.5*y²)
        @constraint(model, y¹ + x² + w² - y² == d)
    end
    @stage 3 begin
        @uncertain d
        @variable(model, x³ <= 2)
        @variable(model, w³)
        @recourse(model, y³)
        @objective(model, Min, x³ + 3*w³)
        @constraint(model, y² + x³ + w³ - y³ == d)
    end
end
s₂₁ = @scenario d = 1.0 probability = 0.5
s₂₂ = @scenario d = 3.0 probability = 0.5

s₃₁ = @scenario d = 1.0 probability = 0.25
s₃₂ = @scenario d = 3.0 probability = 0.25
s₃₃ = @scenario d = 1.0 probability = 0.25
s₃₄ = @scenario d = 3.0 probability = 0.25

sp = instantiate(sm, ([s₂₁,s₂₂], [s₃₁,s₃₂,s₃₃,s₃₄]))


text:
# Example of Parameter algebra operations
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

using InvertibleNetworks, LinearAlgebra

# Input
nx = 16
ny = 16
n_channel = 16
n_hidden = 64
batchsize = 2
X1 = randn(Float32, nx, ny, n_channel, batchsize)
X2 = randn(Float32, nx, ny, n_channel, batchsize)

# Networks
logdet=true
permute="full"
N1 = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet); N1.forward(X1)
N2 = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet); N2.forward(X2)

# Collect parameters
θ1 = get_params(N1)
θ2 = get_params(N2)

# Operations
θ1+θ2
randn(Float32)*θ1
norm.(θ1)
dot.(θ1, θ2)
dot(θ1, θ2)

# Set parameters
set_params!(N1, θ2)
Y1 = N1.forward(X1)[1]
Y1_ = N2.forward(X1)[1]
err_rel = norm(Y1-Y1_)/norm(Y1)


text:
using Rasch
using Optim
using Test
using DelimitedFiles
using LinearAlgebra

@testset "Rasch.jl" begin
    X = readdlm(joinpath(dirname(pathof(Rasch)), "../test/lsat.txt"));
    n, m = size(X)
    r = vec(mapslices(sum, X, dims=2))
    s = Float64.(vec(mapslices(sum, X, dims=1)))
    f = [sum(r.==j) for j in 1:m]
    t = Float64.(vec(mapslices(sum, X .* r, dims=1)))
    q = Float64.(vec(mapslices(sum, X .* r .^2, dims=1)))

    res0, delta_hat0, πs0 =  haberman(m,reshape(s,(m,1)),f,zeros(m-1),Optim.Newton(), Optim.Options(f_tol=1e-6, g_tol=0))
    res1, delta_hat1, πs1 =  haberman(m,hcat(s,t),f,zeros(2m-2),Optim.Newton(), Optim.Options(f_tol=1e-6, g_tol=0))
    res2, delta_hat2, πs2 =  haberman(m,hcat(s,t,q),f,zeros(3m-3),Optim.Newton(), Optim.Options(f_tol=1e-6, g_tol=0)).
end


text:
"""
Bridges `CP.BinPacking{CP.FIXED_CAPACITY_BINPACKING}` to `CP.BinPacking{CP.NO_CAPACITY_BINPACKING}` by adding constraints
on the capacity variables.
"""
struct FixedCapacityBinPacking2BinPackingBridge{T} <: MOIBC.AbstractBridge
    bp::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T}}
    capa::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{FixedCapacityBinPacking2BinPackingBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.BinPacking{CP.FIXED_CAPACITY_BINPACKING, T},
) where {T}
    return MOIBC.bridge_constraint(
        FixedCapacityBinPacking2BinPackingBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{FixedCapacityBinPacking2BinPackingBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.BinPacking{CP.FIXED_CAPACITY_BINPACKING, T},
) where {T}
    # Create the simpler BinPacking set.
    bp_set = CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(s.n_bins, s.n_items, s.weights)
    bp = MOI.add_constraint(model, f, bp_set)

    # Add the capacity constraints.
    f_scalars = MOIU.scalarize(f)
    capa = [
        MOI.add_constraint(model, f_scalars[i], MOI.LessThan(s.capacities[i]))
        for i in 1:s.n_bins
    ]

    return FixedCapacityBinPacking2BinPackingBridge(bp, capa)
end

function MOI.supports_constraint(
    ::Type{FixedCapacityBinPacking2BinPackingBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.BinPacking{CP.FIXED_CAPACITY_BINPACKING, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:FixedCapacityBinPacking2BinPackingBridge})
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{FixedCapacityBinPacking2BinPackingBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOI.get(
    ::FixedCapacityBinPacking2BinPackingBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
    },
) where {T}
    return 1
end

function MOI.get(
    b::FixedCapacityBinPacking2BinPackingBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {T}
    return length(b.capa)
end

function MOI.get(
    b::FixedCapacityBinPacking2BinPackingBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
    },
) where {T}
    return [b.bp]
end

function MOI.get(
    b::FixedCapacityBinPacking2BinPackingBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {T}
    return copy(b.capa)
end


text:
"""
    GroupManifold{𝔽,M<:AbstractManifold{𝔽},O<:AbstractGroupOperation} <: AbstractDecoratorManifold{𝔽}

Decorator for a smooth manifold that equips the manifold with a group operation, thus making
it a Lie group. See [`IsGroupManifold`](@ref) for more details.

Group manifolds by default forward metric-related operations to the wrapped manifold.

# Constructor

    GroupManifold(manifold, op)
"""
struct GroupManifold{𝔽,M<:AbstractManifold{𝔽},O<:AbstractGroupOperation} <:
       AbstractDecoratorManifold{𝔽}
    manifold::M
    op::O
end

@inline function active_traits(f, M::GroupManifold, args...)
    return merge_traits(
        IsGroupManifold(M.op),
        active_traits(f, M.manifold, args...),
        IsExplicitDecorator(),
    )
end

decorated_manifold(G::GroupManifold) = G.manifold

(op::AbstractGroupOperation)(M::AbstractManifold) = GroupManifold(M, op)
function (::Type{T})(M::AbstractManifold) where {T<:AbstractGroupOperation}
    return GroupManifold(M, T())
end

function inverse_retract(
    ::TraitList{<:IsGroupManifold},
    G::GroupManifold,
    p,
    q,
    method::GroupLogarithmicInverseRetraction,
)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = log_lie(G, pinvq)
    return translate_diff(G, p, Identity(G), Xₑ, conv)
end

function inverse_retract!(
    ::TraitList{<:IsGroupManifold},
    G::GroupManifold,
    X,
    p,
    q,
    method::GroupLogarithmicInverseRetraction,
)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = log_lie(G, pinvq)
    return translate_diff!(G, X, p, Identity(G), Xₑ, conv)
end

function is_point(
    ::TraitList{<:IsGroupManifold},
    G::GroupManifold,
    e::Identity,
    te=false;
    kwargs...,
)
    ie = is_identity(G, e; kwargs...)
    (te && !ie) && throw(DomainError(e, "The provided identity is not a point on $G."))
    return ie
end

function is_vector(
    t::TraitList{<:IsGroupManifold},
    G::GroupManifold,
    e::Identity,
    X,
    te=false,
    cbp=true;
    kwargs...,
)
    if cbp
        ie = is_identity(G, e; kwargs...)
        (te && !ie) && throw(DomainError(e, "The provided identity is not a point on $G."))
        (!te) && return ie
    end
    return is_vector(G.manifold, identity_element(G), X, te, false; kwargs...)
end

Base.show(io::IO, G::GroupManifold) = print(io, "GroupManifold($(G.manifold), $(G.op))")


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# TODO: Upstream, LoopVectorization
#       At the time of writing, LoopVectorization.jl cannot handle this kind of
#       loop optimally when passing our custom functions `cons2entropy` and
#       `entropy2cons`. Thus, we need to insert the physics directly here to
#       get a significant runtime performance improvement.
function entropy_projection!(cache, u::StructArray{<:SVector{4, <:Union{Float32, Float64}}},
                             mesh::VertexMappedMesh, equations::CompressibleEulerEquations2D, dg::DGMultiFluxDiff{<:Polynomial})
  rd = dg.basis
  @unpack Vq = rd
  @unpack VhP, u_values, entropy_var_values = cache
  @unpack projected_entropy_var_values, entropy_projected_u_values = cache

  apply_to_each_field(mul_by!(Vq), u_values, u)

  # The following is semantically equivalent to
  # @threaded for i in Base.OneTo(length(u_values))
  #   entropy_var_values[i] = cons2entropy(u_values[i], equations)
  # end
  # but much more efficient due to explicit optimization via `@turbo` from
  # LoopVectorization.jl.
  @unpack gamma, inv_gamma_minus_one = equations
  let
    rho_values, rho_v1_values, rho_v2_values, rho_e_values = StructArrays.components(u_values)
    w1_values, w2_values, w3_values, w4_values = StructArrays.components(entropy_var_values)

    @turbo thread=true for i in indices(
        (rho_values, rho_v1_values, rho_v2_values, rho_e_values,
         w1_values, w2_values, w3_values, w4_values))
      rho    = rho_values[i]
      rho_v1 = rho_v1_values[i]
      rho_v2 = rho_v2_values[i]
      rho_e  = rho_e_values[i]

      # The following is basically the same code as in `cons2entropy`
      v1 = rho_v1 / rho
      v2 = rho_v2 / rho
      v_square = v1^2 + v2^2
      p = (gamma - 1) * (rho_e - 0.5 * rho * v_square)
      s = log(p) - gamma * log(rho)
      rho_p = rho / p

      w1_values[i] = (gamma - s) * inv_gamma_minus_one - 0.5 * rho_p * v_square
      w2_values[i] = rho_p * v1
      w3_values[i] = rho_p * v2
      w4_values[i] = -rho_p
    end
  end

  # "VhP" fuses the projection "P" with interpolation to volume and face quadrature "Vh"
  apply_to_each_field(mul_by!(VhP), projected_entropy_var_values, entropy_var_values)

  # The following is semantically equivalent to
  # @threaded for i in Base.OneTo(length(projected_entropy_var_values)) #eachindex(projected_entropy_var_values)
  #   entropy_projected_u_values[i] = entropy2cons(projected_entropy_var_values[i], equations)
  # end
  # but much more efficient due to explicit optimization via `@turbo` from
  # LoopVectorization.jl.
  let
    rho_values, rho_v1_values, rho_v2_values, rho_e_values = StructArrays.components(entropy_projected_u_values)
    w1_values, w2_values, w3_values, w4_values = StructArrays.components(projected_entropy_var_values)

    @turbo thread=true for i in indices(
        (rho_values, rho_v1_values, rho_v2_values, rho_e_values,
         w1_values, w2_values, w3_values, w4_values))

      # The following is basically the same code as in `entropy2cons`
      # Convert to entropy `-rho * s` used by
      # - See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
      #   [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
      # instead of `-rho * s / (gamma - 1)`
      gamma_minus_one = gamma - 1
      w1 = gamma_minus_one * w1_values[i]
      w2 = gamma_minus_one * w2_values[i]
      w3 = gamma_minus_one * w3_values[i]
      w4 = gamma_minus_one * w4_values[i]

      # s = specific entropy, eq. (53)
      s = gamma - w1 + (w2^2 + w3^2) / (2 * w4)

      # eq. (52)
      rho_iota = (gamma_minus_one / (-w4)^gamma)^(inv_gamma_minus_one) * exp(-s * inv_gamma_minus_one)

      # eq. (51)
      rho_values[i]    = -rho_iota * w4
      rho_v1_values[i] =  rho_iota * w2
      rho_v2_values[i] =  rho_iota * w3
      rho_e_values[i]  =  rho_iota * (1 - (w2^2 + w3^2) / (2 * w4))
    end
  end
end

end # @muladd


text:
"Return the GxB type suffix of a Julia datatype as a String."
function suffix(T)
    if T == Bool
        return "BOOL"
    elseif T == Int8
        return "INT8"
    elseif T == UInt8
        return "UINT8"
    elseif T == Int16
        return "INT16"
    elseif T == UInt16
        return "UINT16"
    elseif T == Int32
        return "INT32"
    elseif T == UInt32
        return "UINT32"
    elseif T == Int64
        return "INT64"
    elseif T == UInt64
        return "UINT64"
    elseif  T == Float32
        return "FP32"
    elseif T == Float64
        return "FP64"
    elseif T == ComplexF32
        return "FC32"
    elseif T == ComplexF64
        return "FC64"
    else
        throw(ArgumentError("Not a valid GrB data type"))
    end
end

function towrappertype(T)
    if T == Bool
        return :Bool
    elseif T == Int8
        return :Int8
    elseif T == UInt8
        return :UInt8
    elseif T == Int16
        return :Int16
    elseif T == UInt16
        return :UInt16
    elseif T == Int32
        return :Int32
    elseif T == UInt32
        return :UInt32
    elseif T == Int64
        return :Int64
    elseif T == UInt64
        return :UInt64
    elseif  T == Float32
        return :Cfloat
    elseif T == Float64
        return :Cdouble
    elseif T == ComplexF32
        return :GxB_FC32_t
    elseif T == ComplexF64
        return :GxB_FC64_t
    else
        throw(ArgumentError("Not a valid GrB data type"))
    end
end

"Load a global constant from SSGrB, optionally specify the resulting pointer type."
function load_global(str, type::Type{Ptr{T}} = Ptr{Nothing}) where {T}
    x =
    try
        dlsym(libgraphblas_handle[], str)
    catch e
        @warn "Symbol not available " * str
        return type(C_NULL)
    end
    return unsafe_load(cglobal(x, type))
end

load_global(str, type) = load_global(str, Ptr{type})

isGxB(name) = occursin("GxB", name)
isGrB(name) = occursin("GrB", name)

function splitconstant(str)
    return String.(split(str, "_"))
end

function tolist(vec)
    vec = sort(string.(vec))
    s = ""
    for i in 1:(length(vec) - 1)
        s = s * "`$(vec[i])`, "
    end
    s = s * "`$(vec[end])`"
    return s
end


text:
using Test
using ProbabilisticCircuits

include("../helper/pc_equals.jl")

@testset "Load and save a small JPC" begin

    function test_my_circuit(pc)

        @test pc isa ProbCircuit

        # Testing number of nodes and parameters
        @test  9 == num_parameters(pc)
        @test 20 == num_nodes(pc)

        # Testing Read Parameters
        EPS = 1e-7
        or1 = children(children(pc)[1])[2]
        @test abs(or1.log_probs[1] - (-1.6094379124341003)) < EPS
        @test abs(or1.log_probs[2] - (-1.2039728043259361)) < EPS
        @test abs(or1.log_probs[3] - (-0.916290731874155))  < EPS
        @test abs(or1.log_probs[4] - (-2.3025850929940455)) < EPS

        or2 = children(children(pc)[1])[1]
        @test abs(or2.log_probs[1] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[2] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[3] - (-2.3025850929940455))  < EPS
        @test abs(or2.log_probs[4] - (-0.35667494393873245)) < EPS

        @test abs(pc.log_probs[1] - (0.0)) < EPS
        @test check_parameter_integrity(pc)
    end

    # first load circuit from PSDD file
    paths = (zoo_psdd_file("little_4var.psdd"), zoo_vtree_file("little_4var.vtree"))
    formats = (PsddFormat(), VtreeFormat())
    pc1 = read(paths, StructProbCircuit, formats)

    mktempdir() do tmp

        # write as a unstructured logic circuit
        jpc_path = "$tmp/example.jpc"
        write(jpc_path, pc1)

        # read as a unstructured logic circuit
        pc2 = read(jpc_path, ProbCircuit)

        test_my_circuit(pc2)
        test_pc_equals(pc1, pc2)

        # write with vtree
        vtree_path = "$tmp/example.vtree"
        paths = (jpc_path, vtree_path)
        write(paths, pc1)

        # read as a structured probabilistic circuit
        pc3 = read(paths, StructProbCircuit)

        @test pc3 isa StructProbCircuit
        test_my_circuit(pc3)
        test_pc_equals(pc1, pc3)
        @test vtree(pc1) == vtree(pc3)

        # read/write compressed
        write("$jpc_path.gz", pc1)
        pc2 = read("$jpc_path.gz", ProbCircuit)

        test_my_circuit(pc2)
        test_pc_equals(pc1, pc2)

        # read/write compressed structured
        paths = ("$jpc_path.gz", vtree_path)
        write(paths, pc1)
        pc3 = read(paths, StructProbCircuit)

        @test pc3 isa StructProbCircuit
        test_my_circuit(pc3)
        test_pc_equals(pc1, pc3)
        @test vtree(pc1) == vtree(pc3)

    end

end

@testset "Can save JPCs with nonbinary multiplications" begin

    pc1 = fully_factorized_circuit(ProbCircuit, 10)

    mktempdir() do tmp

        jpc_path = "$tmp/example.jpc"
        write(jpc_path, pc1)

        pc2 = read(jpc_path, ProbCircuit)

        test_pc_equals(pc1, pc2)
    end

end