text:
module TestExamples3D

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "3d")

# Run basic tests
@testset "Examples 3D" begin
  @testset "parameters.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters.toml"),
            l2   = [0.00015975754755823664],
            linf = [0.001503873297666436])
  end
  @testset "parameters.toml with restart and t_end=2" begin
    Trixi.run(joinpath(EXAMPLES_DIR, "parameters.toml"))
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters.toml"),
            l2   = [0.00017800012466353434],
            linf = [0.001452075263740804],
            t_end = 2, restart = true, restart_filename = "out/restart_000040.h5")
  end
  @testset "parameters.toml with initial_conditions_sin_periodic" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters.toml"),
            l2   = [0.002727292086517533],
            linf = [0.024833049753677727],
            initial_conditions=Trixi.initial_conditions_sin_periodic)
  end
  @testset "parameters.toml with initial_conditions_constant" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters.toml"),
            l2   = [9.770171014620371e-16],
            linf = [2.4424906541753444e-15],
            initial_conditions=Trixi.initial_conditions_constant)
  end
  @testset "parameters.toml with initial_conditions_linear_z and periodicity=false" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters.toml"),
            l2   = [6.607840408143593e-16],
            linf = [5.773159728050814e-15],
            initial_conditions=Trixi.initial_conditions_linear_z,
            boundary_conditions=Trixi.boundary_conditions_linear_z, periodicity=false)
  end
  @testset "parameters_source_terms.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_source_terms.toml"),
            l2   = [0.010323099666828388, 0.00972876713766357, 0.00972876713766343, 0.009728767137663324, 0.015080409341036285],
            linf = [0.034894880154510144, 0.03383545920056008, 0.033835459200560525, 0.03383545920054587, 0.06785780622711979])
  end
  @testset "parameters_source_terms.toml with split_form" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_source_terms.toml"),
            l2   = [0.010323099666828388, 0.00972876713766357, 0.00972876713766343, 0.009728767137663324, 0.015080409341036285],
            linf = [0.034894880154510144, 0.03383545920056008, 0.033835459200560525, 0.03383545920054587, 0.06785780622711979],
            volume_integral_type = "split_form")
  end
  @testset "parameters_eoc_test_euler.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_euler.toml"),
            l2   = [0.000363723832448333, 0.00039555684672049366, 0.0003955568467203738, 0.00039555684672064724, 0.0007811604790242773],
            linf = [0.002400072140187337, 0.0029635489437536133, 0.0029635489437540574, 0.0029635489437565, 0.007191455734479657])
  end
  @testset "parameters_eoc_test_coupled_euler_gravity.toml with resid_tol = 1.0e-4, t_end = 0.2" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_eoc_test_coupled_euler_gravity.toml"),
            l2   = [0.00042767972112699913, 0.00047204316046796835, 0.00047204316046784795, 0.0004720431604680035, 0.0010987015429634586, 0.00012296598036447797, 0.0005745341792812197, 0.0005745341792812442, 0.0005745341792812238],
            linf = [0.0034966337661186397, 0.0037643976198782347, 0.003764397619878901, 0.0037643976198780127, 0.008370354378078648, 0.0010129211321238465, 0.0024406779290754455, 0.002440677929075438, 0.0024406779290755756],
            resid_tol = 1.0e-4, t_end = 0.2)
  end
  @testset "parameters_mortar.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_mortar.toml"),
            l2   = [0.0018461483161353273],
            linf = [0.017728496545256434])
  end
  @testset "parameters_mortar_euler.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_mortar_euler.toml"),
            l2   = [0.0019011097431965655, 0.0018289464087588392, 0.0018289464087585998, 0.0018289464087588862, 0.003354766311541738],
            linf = [0.011918594206950184, 0.011808582644224241, 0.011808582644249999, 0.011808582644239785, 0.02464803617735356])
  end
  @testset "parameters_amr.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_amr.toml"),
            l2   = [9.773858425669403e-6],
            linf = [0.0005853874124926092])
  end
  @testset "parameters_amr_euler.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_amr_euler.toml"),
            l2   = [0.00382819196730758, 0.0038281919673075725, 0.0038281919673075746, 0.0038281919673075738, 0.0057422879509614905],
            linf = [0.07390560349428554, 0.07390560349428577, 0.07390560349428621, 0.07390560349428643, 0.11085840524143098],
            t_end=0.1)
  end
  @testset "parameters_amr_blob.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_amr_blob.toml"),
            l2   = [0.04641288111176061, 0.25587674204071276, 0.03358077695134774, 0.033580776951347646, 1.0552376897521447],
            linf = [9.773023537783775, 48.58018285242446, 9.56936427666669, 9.56936427666669, 170.97251488955624],
            t_end=0.2)
  end
  @testset "parameters_ec.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec.toml"),
            l2   = [0.025101741317688664, 0.01655620530022176, 0.016556205300221737, 0.016549388264402515, 0.09075092792976944],
            linf = [0.43498932208478724, 0.2821813924028202, 0.28218139240282025, 0.2838043627560838, 1.5002293438086647])
  end
  @testset "parameters_ec.toml with initial_conditions=Trixi.initial_conditions_constant" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec.toml"),
            l2   = [5.717218008425079e-16, 6.088971423170968e-16, 6.23130776282275e-16, 7.29884557381127e-16, 5.167198077601542e-15],
            linf = [3.885780586188048e-15, 4.454769886308441e-15, 3.219646771412954e-15, 4.884981308350689e-15, 4.440892098500626e-14],
            initial_conditions=Trixi.initial_conditions_constant)
  end
  @testset "parameters_ec.toml with flux_chandrashekar" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec.toml"),
            l2   = [0.025105743648126774, 0.016571417754430256, 0.01657141775443023, 0.016565202090289916, 0.09077232065771225],
            linf = [0.4349225166034201, 0.27945714200874, 0.2794571420087401, 0.28021366413271664, 1.5240679700745954],
            surface_flux=flux_chandrashekar, volume_flux=flux_chandrashekar)
  end
  @testset "parameters_ec.toml with flux_kennedry_gruber" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec.toml"),
            l2   = [0.025120431810845507, 0.016599310737401483, 0.01659931073740148, 0.016592567464138185, 0.090856457771812],
            linf = [0.43120500632996794, 0.28419288751363336, 0.2841928875136334, 0.28583515705222146, 1.515485025725378],
            surface_flux=flux_kennedy_gruber, volume_flux=flux_kennedy_gruber)
  end
  @testset "parameters_ec.toml with flux_shima_etal" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec.toml"),
            l2   = [0.025099944530993942, 0.016561611274319134, 0.016561611274319127, 0.01655478190136039, 0.09076538812894279],
            linf = [0.43472962954165273, 0.2824065323711477, 0.2824065323711474, 0.28409419760015847, 1.4995295774522692],
            surface_flux=flux_shima_etal, volume_flux=flux_shima_etal)
  end
  @testset "parameters_density_pulse.toml with source_terms_harmonic" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_density_pulse.toml"),
            l2   = [0.05719652660597408, 0.0571965266059741, 0.05719652660597407, 0.05719652660597409, 0.08579478990896279],
            linf = [0.27375961853433606, 0.27375961853433517, 0.27375961853433384, 0.2737596185343343, 0.4106394278015033],
            source_terms=Trixi.source_terms_harmonic,
            extra_analysis_quantities=["l2_error_primitive", "linf_error_primitive"])
  end
  @testset "parameters_taylor_green_vortex.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_taylor_green_vortex.toml"),
            l2   = [0.0003494971047256544, 0.03133386380969968, 0.031333863809699644, 0.04378595081016185, 0.015796569210801217],
            linf = [0.0013934701399120897, 0.07284947983025436, 0.07284947983025408, 0.12803234075782724, 0.07624639122292365],
            t_end = 0.5)
  end
  @testset "parameters_shock_capturing.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_shock_capturing.toml"),
            l2   = [0.025558219399128387, 0.01612806446620796, 0.016128064466207948, 0.016120400619198158, 0.09208276987000782],
            linf = [0.3950327737713353, 0.26324766244272796, 0.2632476624427279, 0.2634129727753079, 1.371321006006725])
  end
  @testset "parameters_shock_capturing.toml with initial_conditions_sedov_blast_wave" begin
    # OBS! This setup does not run longer but crashes (also the parameters do not make sense) -> only for testing the IC!
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_shock_capturing.toml"),
            l2   = [0.03627060784392582, 0.05178777376859809, 0.05178777376859804, 0.05178777376859802, 0.23043996953698023],
            linf = [0.9307998233177583, 1.4326649193439467, 1.4326649193439467, 1.4326649193439467, 12.80585041235138],
            initial_conditions=Trixi.initial_conditions_sedov_blast_wave, cfl=0.25, shock_alpha_max=1.0, t_end=0.1)
  end
  @testset "parameters_shock_capturing.toml with initial_conditions_sedov_self_gravity" begin
    # OBS! This setup does not run longer but crashes (also the parameters do not make sense) -> only for testing the IC!
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_shock_capturing.toml"),
            l2   = [0.04846527000320781, 0.051787773760055514, 0.051787773760055486, 0.05178777376005548, 0.23043996953467236],
            linf = [0.9307979866990295, 1.4326649193456429, 1.4326649193456429, 1.4326649193456429, 12.805850412386896],
            initial_conditions=Trixi.initial_conditions_sedov_self_gravity, cfl=0.25, shock_alpha_max=1.0, t_end=0.1)
  end
  @testset "parameters_shock_capturing.toml with amr_indicator=`sedov_self_gravity`" begin
    # OBS! This setup does not make much practical sense. It is only added to exercise the
    # `sedov_self_gravity` AMR indicator, which in its original configuration is too expensive for
    # CI testing
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_shock_capturing.toml"),
            l2   = [0.022890627324485553, 0.013353127563710173, 0.013353127563710156, 0.013352146925243637, 0.08348119429398775],
            linf = [0.3769085161148348, 0.380358363190641, 0.3803583631906434, 0.380366775575835, 1.3162027903728162],
            n_steps_max=10, amr_interval=1, amr_indicator="sedov_self_gravity", max_refinement_level=4)
  end
  @testset "parameters_hyp_diff_llf.toml with initial_refinement_level=2" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_hyp_diff_llf.toml"),
            l2   = [0.0015303292770225546, 0.011314166522881952, 0.011314166522881981, 0.011314166522881947],
            linf = [0.022634590339093097, 0.10150613595329361, 0.10150613595329361, 0.10150613595329361],
            initial_refinement_level=2)
  end
  @testset "parameters_hyp_diff_llf.toml with initial_refinement_level=2, surface_flux=Trixi.flux_upwind)" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_hyp_diff_llf.toml"),
            l2   = [0.0015377708559180534, 0.011376842329542572, 0.011376842329542624, 0.0113768423295426],
            linf = [0.02271542063004106, 0.10191067906109286, 0.10191067906109552, 0.10191067906109286],
            initial_refinement_level=2, surface_flux=Trixi.flux_upwind)
  end
  @testset "parameters_hyp_diff_nonperiodic.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_hyp_diff_nonperiodic.toml"),
            l2   = [0.00022868324220593294, 0.0007974310370259415, 0.0015035143239197598, 0.0015035143239198418],
            linf = [0.0016329580288680923, 0.0029870270738030775, 0.009177053066089513, 0.009177053066084184])
  end
  @testset "parameters_ec_mhd.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec_mhd.toml"),
            l2   = [0.01921453037426997, 0.01924853398980921, 0.01924853398980923, 0.019247118340533328, 0.08310482412935676, 0.010362656540935251, 0.010362656540935237, 0.010364587080559528, 0.00020760700572485828],
            linf = [0.2645851360519166, 0.33611482816103344, 0.33611482816103466, 0.36952265576762666, 1.230825809630423, 0.09818527443798974, 0.09818527443798908, 0.10507242371450054, 0.008456471524217968])
  end
  @testset "parameters_ec_mhd.toml with initial_conditions=Trixi.initial_conditions_constant" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_ec_mhd.toml"),
            l2   = [4.850506049646793e-16, 2.4804155700127237e-15, 3.579471462379534e-15, 2.7395862184339726e-15, 2.4916602560342516e-14, 1.669368799061149e-15, 1.4052897861706032e-15, 1.0685989093080367e-15, 1.1611070325375158e-15],
            linf = [3.552713678800501e-15, 1.4710455076283324e-14, 2.3814283878209608e-14, 2.6423307986078726e-14, 1.6342482922482304e-13, 1.1546319456101628e-14, 1.0880185641326534e-14, 1.4099832412739488e-14, 1.1483287543575534e-14],
            atol = 1000*eps(),
            initial_conditions=Trixi.initial_conditions_constant)
  end
  @testset "parameters_alfven_wave.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_alfven_wave.toml"),
            l2   = [0.0038729054515012624, 0.00903693761037057, 0.0041729297273898815, 0.01160504558506348, 0.006241548790045999, 0.009227641613254402, 0.0034580608435846143, 0.011684993365513006, 0.0022068452165023645],
            linf = [0.012628629484152443, 0.03265276295369954, 0.012907838374176334, 0.044746702024108326, 0.02796611265824822, 0.03453054781110626, 0.010261557301859958, 0.044762592434299864, 0.010012319622784436])
  end
  @testset "parameters_mortar_alfven_wave.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_mortar_alfven_wave.toml"),
            l2   = [0.0021484102061835623, 0.006826504155492453, 0.0030653111370061784, 0.008735898256361025, 0.0051601878379492335, 0.007157480202233399, 0.0028291977973972948, 0.008815052614117018, 0.0022321821323698257],
            linf = [0.012956870409227328, 0.05529249146399706, 0.020854504834048836, 0.05898012498637771, 0.03162799656904003, 0.05512773554440975, 0.017941374395225362, 0.060061114374191496, 0.013036070296136178],
            t_end = 0.25)
  end
  # 3D Orszag-Tang included to exercise all terms in the HLL flux
  @testset "parameters_orszag_tang.toml" begin
    test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_orszag_tang.toml"),
            l2   = [0.0043911605751115424, 0.04144735653371165, 0.04150129965650717, 0.04150353600000829, 0.036931197750736805, 0.021125598820694595, 0.032956068087418154, 0.03296235602392588, 6.318083915607208e-6],
            linf = [0.01789383976134809, 0.08496187610572214, 0.08909116075943745, 0.08505952838326755, 0.10443373959204932, 0.05387852204182135, 0.08812990990777562, 0.07804874749131957, 8.138512446081734e-5],
            t_end = 0.06)
  end  # too expensive for CI
  # @testset "parameters_sedov_shock_capturing.toml with n_steps_max = 2" begin
  #   test_trixi_run(joinpath(EXAMPLES_DIR, "parameters_sedov_shock_capturing.toml"),
  #           l2   = [0.00015213881280510253, 0.001481110249423103, 0.0014811102494231387, 0.001481110249423187, 0.002940437008367858],
  #           linf = [0.03254534843490764, 0.38932044051654113, 0.38932044051654097, 0.38932044051654097, 1.050399588579145],
  #           n_steps_max = 2)
  # end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
module twisted_beam_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule
using FinEtools.MeshExportModule

function twisted_beam_algo()
    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  = 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(IntegData(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, IntegData(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,:],1)
    println("displacement  = $(theutip[dir]) as compared to converged $uex")
    println("normalized displacement  = $(theutip[dir]/uex*100) %")

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

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

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

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

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

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

    println("Done")
    true
end # twisted_beam_algo

function twisted_beam_algo_stress()
    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 ut-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 = 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(IntegData(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, IntegData(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,:],1)
    println("displacement  = $(theutip[dir]) as compared to converged $uex")
    println("normalized displacement  = $(theutip[dir]/uex*100) %")

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

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

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

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

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

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

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

    # 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)])")

    # Write out mesh with principal stresses, elementwise
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-press-ew",
    "quantity"=> :pressure, "component"=> 1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    ps  = modeldata["postprocessing"]["exported"][1]["field"]
    println("extremes of pressure: $([minimum(ps.values),   maximum(ps.values)])")

    println("Done")
    true
end # twisted_beam_algo_stress

function twisted_beam_export()
    println("""
    Refer to twisted_beam.jl.

    This example EXPORTS the model to Abaqus. Import the .inp file
    into Abaqus and run the job.
    """)
    E = 0.29e8;
    nu = 0.22;
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl = 2; nt = 1; nw = 1; ref = 5;
    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(IntegData(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, IntegData(fes, GaussRule(3,2)),
    material))

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

    AE = AbaqusExporter("twisted_beam");
    HEADING(AE, "Twisted beam example");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "c3d8rh", "AllElements", 1, region1["femm"].integdata.fes.conn)
    ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdata.fes), flux1["femm"].integdata.fes.conn)
    NSET_NSET(AE, "l1", l1)
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    SURFACE_SECTION(AE, "TractionElements")
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3)
    DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
    END_STEP(AE)
    close(AE)

    true
end # twisted_beam_export

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

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

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

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

    # Traction on the opposite edge
    boundaryfes  =   meshboundary(fes);
    Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
    el1femm  = FEMMBase(IntegData(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, IntegData(fes, GaussRule(3,2)),
    material))

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

    AE = AbaqusExporter("twisted_beam");
    # AE.ios = STDOUT;
    HEADING(AE, "Twisted beam example");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "c3d8rh", "AllElements", 1, region1["femm"].integdata.fes.conn)
    ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdata.fes), flux1["femm"].integdata.fes.conn)
    NSET_NSET(AE, "l1", l1)
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    SURFACE_SECTION(AE, "TractionElements")
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2)
    BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3)
    DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
    END_STEP(AE)
    close(AE)

end # twisted_beam_export_nb

function twisted_beam_msh8()
    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(IntegData(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, IntegData(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,:],1)
    println("displacement  = $(theutip[dir]) as compared to converged $uex")

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

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

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

    # 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)])")

    # 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)])")

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

    println("Done")
    true
end # twisted_beam_msh8

function twisted_beam_msh8_algo_stress()
    println("""
    The initially twisted cantilever beam is one of the standard test    problems for verifying the finite-element accuracy [1]. The beam is        clamped at one end and loaded either with unit in-plane or        unit out-of-plane force at the other. The centroidal axis of the beam s
    straight at the undeformed  configuration, while its cross-sections are
    twisted about the centroidal axis from 0 at the clamped end to pi/2 at
    the free end.

    Reference:
    Zupan D, Saje M (2004) On "A proposed standard set of problems to test
    finite element accuracy": the twisted beam. Finite Elements in Analysis
    and Design 40: 1445-1451.
    """)
    E = 0.29e8;
    nu = 0.22;
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl = 2; nt = 1; nw = 1; ref = 4;
    p =   1/W/t;
    #  Loading in the Z direction
    loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
    #   Loading in the Y direction
    #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
    tolerance  = t/1000;

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

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

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

    # Traction on the opposite edge
    boundaryfes  =   meshboundary(fes);
    Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
    el1femm  = FEMMBase(IntegData(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, IntegData(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,:],1)
    println("displacement  = $(theutip[dir]) as compared to converged $uex")
    println("normalized displacement  = $(theutip[dir]/uex*100) %")

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

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

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

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

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

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

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

    # 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)])")

    # Write out mesh with principal stresses, elementwise
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam-press-ew",
    "quantity"=> :pressure, "component"=> 1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    ps  = modeldata["postprocessing"]["exported"][1]["field"]
    println("extremes of pressure: $([minimum(ps.values),   maximum(ps.values)])")

    println("Done")
    true
end # twisted_beam_msh8_algo_stress

function allrun()
    println("#####################################################")
    println("# twisted_beam_algo ")
    twisted_beam_algo()
    println("#####################################################")
    println("# twisted_beam_algo_stress ")
    twisted_beam_algo_stress()
    println("#####################################################")
    println("# twisted_beam_export ")
    twisted_beam_export()
    println("#####################################################")
    println("# twisted_beam_export_nb ")
    twisted_beam_export_nb()
    println("#####################################################")
    println("# twisted_beam_msh8 ")
    twisted_beam_msh8()
    println("#####################################################")
    println("# twisted_beam_msh8_algo_stress ")
    twisted_beam_msh8_algo_stress()
    return true
end # function allrun

end # module twisted_beam_examples

ext: jl


text:
"""
TestTupleStructs is necessary as a module because we define new types and want to be able to include this file multiple times without generating redefinition of constant errors.

These tests make a test file test.jld in order to ensure that the resulting types are able to be stored in JLD format.

"""
module TestTupleStructs
using Base.Test
using TupleStructs

ex = quote
    type S1
        a::Int64;
        b::Int
    end
end
eval(ex)
p = S1(1,3)
@test p.a == 1
@test p.b == 3

ts = tostruct(:S, Tuple{Int, Int})
eval(ts)
q = S(1,3)
@test q.f1 == 1
@test q.f2 == 3

@show ts = tostruct(:T, Tuple{Int64, Int, ASCIIString, Tuple{Int, Float64}})
eval(ts)
@show t= T(1, 2, "hello", (3, 0.0))
@test t.f1 == 1
@test t.f2 == 2
@test t.f3 == "hello"
@test t.f4 == (3,0.0)

function isequ(t::TestTupleStructs.T, s::TestTupleStructs.T)
    return (t.f1 == s.f1) && (t.f2 == s.f2) && (t.f3 == s.f3) && (t.f4 == s.f4)
end

import Base.==

function ==(t::TestTupleStructs.S1, s::TestTupleStructs.S1)
    return (t.a == s.a) && (t.b == s.b)
end
function ==(t::TestTupleStructs.S, s::TestTupleStructs.S)
    return (t.f1 == s.f1) && (t.f2 == s.f2)
end

using JLD

t2 = t
q1 = q
p1 = p
@save "test.jld" t
@save "test.jld" p
@save "test.jld" q

@load "test.jld"
@test isequ(t2,t)
@test p == p1
@test q == q1

end


text:
@testset "Mixture density" begin

    @testset "Limit concentrations leads to pure substance?" begin
        water = Water()
        brine860 = BrineNaCl(8.6)
        brine920 = BrineNaCl(9.2)
        barite = Barite()

        Ps = LinRange(101325.0, 100.0E6, 5)
        Ts = LinRange(283.15, 413.15, 5)
        @testset "Water + Barite mixture" begin
            FList = [water]
            SLits = [barite]
            fs = [1.0, 0.0]

            WatBar0 = DrillFluid(FList, SLits, fs = fs)

            fs = [0.0, 1.0]
            WatBar100 = DrillFluid(FList, SLits, fs = fs)

            for P in Ps
                for T in Ts
                    rho_water = rho(P, T, water)
                    rho_barite = rho(P, T, barite)
                    rho_mix0 = rho(P, T, WatBar0)
                    rho_mix100 = rho(P, T, WatBar100)

                    @test rho_water == rho_mix0

                    @test rho_barite == rho_mix100

                end
            end

        end

        @testset "Brine 8.6ppg + Barite mixture" begin
            FList = [brine860]
            SLits = [barite]
            fs = [1.0, 0.0]

            BrBar0 = DrillFluid(FList, SLits, fs = fs)

            fs = [0.0, 1.0]
            BrBar100 = DrillFluid(FList, SLits, fs = fs)

            for P in Ps
                for T in Ts
                    rho_Br = rho(P, T, brine860)
                    rho_barite = rho(P, T, barite)
                    rho_mix0 = rho(P, T, BrBar0)
                    rho_mix100 = rho(P, T, BrBar100)

                    @test rho_Br == rho_mix0

                    @test rho_barite == rho_mix100

                end
            end

        end

        @testset "Brine 9.2ppg + Barite mixture" begin
            FList = [brine920]
            SLits = [barite]
            fs = [1.0, 0.0]

            BrBar0 = DrillFluid(FList, SLits, fs = fs)

            fs = [0.0, 1.0]
            BrBar100 = DrillFluid(FList, SLits, fs = fs)

            for P in Ps
                for T in Ts
                    rho_Br = rho(P, T, brine920)
                    rho_barite = rho(P, T, barite)
                    rho_mix0 = rho(P, T, BrBar0)
                    rho_mix100 = rho(P, T, BrBar100)

                    @test rho_Br == rho_mix0

                    @test rho_barite == rho_mix100

                end
            end

        end

    end

end


text:
# Helpers for benchmark

type BenchmarkTable
    name::ASCIIString

    colnames::Vector{ASCIIString}
    rownames::Vector{ASCIIString}

    colmap::Dict{ASCIIString, Int}
    rowmap::Dict{ASCIIString, Int}
    rows::Vector{Vector{Float64}}

    function BenchmarkTable(name::ASCIIString, colnames::Vector{ASCIIString})
        colmap = Dict{ASCIIString, Int}([colnames[i] => i for i in 1 : length(colnames)])
        new(name, colnames, ASCIIString[],
            colmap, Dict{ASCIIString, Int}(), Array(Vector{Float64}, 0))
    end
end

nrows(tab::BenchmarkTable) = length(tab.rownames)
ncolumns(tab::BenchmarkTable) = length(tab.colnames)

getindex(tab::BenchmarkTable, i::Int, j::Int) = tab.rows[i][j]
row(tab::BenchmarkTable, i::Int) = tab.rows[i]

function add_row!(tab::BenchmarkTable, rowname::ASCIIString, row::Vector{Float64})
    if haskey(tab.rowmap, rowname)
        error("The rowname $rowname has existed.")
    end
    nc = ncolumns(tab)
    if length(row) != nc
        error("Invalid row length.")
    end

    nr = nrows(tab)
    push!(tab.rownames, rowname)
    tab.rowmap[rowname] = nr + 1

    push!(tab.rows, row)
    tab
end

function _leftalign_string(s::ASCIIString, len::Int)
    slen = length(s)
    slen < len ? string(s, repeat(" ", len - slen)) :
    slen > len ? s[1:len] : s
end

function _rightalign_string(s::ASCIIString, len::Int)
    slen = length(s)
    slen < len ? string(repeat(" ", len - slen), s) :
    slen > len ? s[1:len] : s
end

function Base.show(io::IO, tab::BenchmarkTable)
    println(io, "BenchmarkTable: $(tab.name)")

    # compute lengths
    lr = maximum([length(x) for x in tab.rownames])
    lc = maximum([length(x) for x in tab.colnames])
    lc = max(lc, 8)

    m = nrows(tab)
    n = ncolumns(tab)

    # print headers
    hlen = lr + (lc + 3) * n + 3
    println(io, repeat("=", hlen))
    print(io, repeat(" ", lr + 2))
    print(io, "| ")

    for cname in tab.colnames
        print(io, _leftalign_string(cname, lc))
        print(io, " | ")
    end
    println(io)
    println(io, repeat("-", hlen))

    # print rows
    for i in 1 : m
        rname = tab.rownames[i]
        row = tab.rows[i]

        print(io, _leftalign_string(rname, lr))
        print(io, "  | ")
        for j in 1 : n
            vs = @sprintf("%8.4f", row[j])
            print(io, _rightalign_string(vs, lc))
            print(io, " | ")
        end
        println(io)
    end
    println(io, repeat("-", hlen))
end


text:
using GLVisualize, GLAbstraction, GeometryTypes, Reactive, Colors, GLFW
function circlescale(v0, pressed)
    scale = v0
    if isempty(pressed)
        for i in eachindex(scale)
            scale[i] = Vec2f0(20)
        end
    else
        ls = linspace(1.0, 2.0, length(scale))
        for i in eachindex(scale)
            scale[i] = Vec2f0(20*ls[i])
        end
    end
    scale
end
mpos(x) = translationmatrix(Vec3f0(x, 0))
function to_color(button)
    if button == GLFW.MOUSE_BUTTON_LEFT
        return GLVisualize.default(RGBA)
    elseif button == GLFW.MOUSE_BUTTON_MIDDLE
        return RGBA{Float32}(0.9,0.1, 0.3, 1.0)
    elseif button == GLFW.MOUSE_BUTTON_MIDDLE
        return RGBA{Float32}(0.9,0.1, 0.9, 1.0)
    end
    RGBA{Float32}(0.6,0.6,0.6,1.0)
end

function mouse_color(pressed)
    if isempty(pressed)
        return RGBA{Float32}(0,0,0,0.8)
    else
        return to_color(first(pressed))
    end
end

function add_mouse(window)
    N = 2
    points = fill(Point2f0(0), N)

    @materialize mouseposition, mouse_buttons_pressed = window.inputs

    model_matrix = map(mpos, mouseposition)

    scales0 = fill(Vec2f0(20), N)

    scale = foldp(circlescale, scales0, mouse_buttons_pressed)
    color = map(mouse_color, mouse_buttons_pressed)

    _view(visualize(
        (Circle(Point2f0(0), 20f0), points),
        color=RGBA{Float32}(0,0,0,0), stroke_width=2f0,
        stroke_color=color, scale=scale,
        model=model_matrix
    ), window, camera=:fixed_pixel)

end


text:
using UnitTestDesign
using Test
using SafeTestsets

@testset "UnitTestDesign.jl" begin
    include("test_combinations.jl")
    include("test_coverage_matrix.jl")
    include("test_coverage_set.jl")
    include("test_excursions.jl")
    include("test_factorial_interface.jl")
    include("test_full_factorial.jl")
    include("test_greedy_tuples.jl")
    include("test_parameter_order.jl")
end


text:
# Sensitivty Analysis
random_set = collect(.200:.005:.300)
final_MAWP_out = zeros(10000)

for i in 1:size(final_MAWP_out,1)
    print("this is iteration", i,"\n")
# For all assessments - determine the inspection data grid
M1 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M2 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.100, rand(random_set), rand(random_set), rand(random_set), rand(random_set), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M3 = [0.300, 0.300, 0.300, 0.300, 0.300, rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M4 = [0.300, 0.300, 0.300, 0.300, 0.300, rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), rand(random_set), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M5 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M6 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
CTPGrid = hcat(M6,M5,M4,M3,M2,M1) # build in descending order
global CTPGrid = rotl90(CTPGrid) # rotate to correct orientation
CTPGrid[4,6]
# Level 1 fit for service
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]
equipment_group = "piping" # "vessel", "tank"
flaw_location = "external" # "External","Internal"
pipe_code = "B31.3" # "B31.3", "B31.8"
FCA_string = "external"
metal_loss_categorization = "LTA" # "LTA" or "Groove-Like Flaw"
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"
tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.
trd = .3 # uniform thickness away from the local metal loss location established by thickness measurements at the time of the assessment.
FCAml = 0.07 # Future Corrosion Allowance applied to the region of metal loss.
FCA = 0.0 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).
LOSS = 0.0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.
Do = 3.5 # Outside Diameter
D = Do - 2*(tnom)
P = 1480.0 # internal design pressure.
S = 20000.0 # allowable stress.
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.
F = 0.6 # B31.8 design factor
T = 1.0 # temperature factor B31.8
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .
t = trd - LOSS - FCAml # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).
spacings = 0.5 # spacings determine by visual inspection to adequately ccategorizse the corrosion -----------+ may add to CTP_Grid function for plotting purposes
# Flaw dimensions
s = 6.0 # longitudinal extent or length of the region of local metal loss based on future corroded thickness,
c = 2.0 # circumferential extent or length of the region of local metal loss (see Figure 5.2 and Figure 5.10), based on future corroded thickness, tc .
Ec = 1.0 # circumferential weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
El = 1.0 # longitudinal weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9

# For all assessments determine far enough from structural discontinuity
# Flaw-To-Major Structural Discontinuity Spacing
L1msd = [12.0] # distance to the nearest major structural discontinuity.
L2msd = [12.0] # distance to the nearest major structural discontinuity.
L3msd = [12.0] # distance to the nearest major structural discontinuity.
L4msd = [12.0] # distance to the nearest major structural discontinuity.
L5msd = [12.0] # distance to the nearest major structural discontinuity.
Lmsd = minimum([L1msd,L2msd,L3msd,L4msd,L5msd])
if (Lmsd[1] >= (1.8*(sqrt(D*(trd - LOSS - FCA)))))
    print("Satisfied - Flaw is located far enough from structural discontinuity\n")
    lmsd_satisfied = 1
else
    print("Not satisfied - Flaw is too close to a structural discontinuity - Conduct a level 3 assessment\n")
    lmsd_satisfied = 0
end

# Groove Like Flaw dimensions
gl = .05 # length of the Groove-Like Flaw based on future corroded condition.
gw = .4 # width of the Groove-Like Flaw based on future corroded condition.
gr = 0.1 # radius at the base of a Groove-Like Flaw based on future corroded condition.
β = 40.0 # see (Figure 5.4) :: orientation of the groove-like flaw with respect to the longitudinal axis or a parameter to compute an effective fracture toughness for a groove being evaluated as a crack-like flaw, as applicable.

# Perform level 1 assessment
if (part5_applicability[1] == 1 && lmsd_satisfied == 1) # begin level 1 assessment
    #let part_5_lta_output = Array{Any,2},
    part_5_lta_output = Part5LTALevel2(CTPGrid; annex2c_tmin_category=annex2c_tmin_category, equipment_group=equipment_group, flaw_location=flaw_location, FCA_string=FCA_string, metal_loss_categorization=metal_loss_categorization, units=units, Lmsd=Lmsd, tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, tsl=tsl, t=t,spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)
    #end # let end
    part_5_lta_output
    final_MAWP_out[i] = part_5_lta_output[10,2]
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 1 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 1)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
end
end # end random for loop

#output
using DataFrames
using CSV
using Statistics
out = DataFrame(hcat(final_MAWP_out))

labels = ["Average MAWPr", "Maximum MAWPr", "Minimum MAWPr", "Actual MAWPr (known grid)"]
average_MAWPr = mean(final_MAWP_out)
max_MAWPr = maximum(final_MAWP_out)
min_MAWPr = minimum(final_MAWP_out)
actual_MAWPr = part_5_lta_output[10,2]

mawp_stats = [average_MAWPr, max_MAWPr, min_MAWPr, actual_MAWPr]
stats_out = hcat(labels,mawp_stats)

CSV.write("C:/Users/Andrew.Bannerman/OneDrive - Shell/Documents/API CODES/level_2_LTA_grid_assumption_MAWPr_out_case_3.csv", out;delim=',')

# s and c random
s = 6.0 # longitudinal extent or length of the region of local metal loss based on future corroded thickness,
c = 2.0 # circumferential extent or length of the region of local metal loss (see Figure 5.2 and Figure 5.10), based on future corroded thickness, tc .
Ec = 1.0 # circumferential weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
El = 1.0 # longitudinal weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency

s_random = collect(6.0:.5:20.0)
c_random = collect(2.0:.5:12.0)

step_9_satisfied = zeros(5000)
criteria = zeros(5000)
c_param = zeros(5000)

for i in 1:5000
    c = rand(c_random)
    s = rand(s_random)
if (c <= (2*s)*(El/Ec)) # eq (5.13)  # original = 38
    print("eq 5.13 satisfied - no further evaluation is required\n")
    step_9_satisfied[i] = 1.0
    criteria[i] = ((2*s)*(El/Ec))
    c_param[i] = c
else
    print("eq 5.13 not satisfied - Proceed to STEP 9.2\n")
    step_9_satisfied[i] = 0.0
    criteria[i] = ((2*s)*(El/Ec))
    c_param[i] = c
end # end step 9.1
end # end loop


text:
function hello()
    println("Hello world, here I am!")
end

x = 3


text:
"""

```
contrib(s::Selection, atom_contributions::Array{Float64}, selection)
```

Extract the contribution of a given atom type selection from the
solute or solvent atomic contributions to the MDDF.

`s` here is the solute or solvent selection (type `ComplexMixtures.Selection`)
`atom_contributions` is the `R.solute_atom` or `R.solvent_atom` arrays of the `Result` structure,
and the last argument is the selection of atoms from the solute to be considered, given as a list
of indexes, list of atom names, vector of `PDBTools.Atom`s, or a `PDBTools.Residue`.

"""
function contrib(s::Selection, atom_contributions::Array{Float64}, indexes::Vector{Int})
    nbins = size(atom_contributions, 1)
    c = zeros(nbins)
    # If the selection is a single molecule, the indexes are anything
    if s.nmols == 1
        for it in indexes
            ind = findfirst(isequal(it),s.index)
            if isnothing(ind)
                error("Index $it of input list not found in selection indexes list.")
            end
            c += atom_contributions[:, ind]
        end
    # If more than one molecule, the index must correspond to an atom within one molecule
    else
        for it in indexes
            if it > s.natomspermol
                error(
                    "The index list contains atoms with indexes greater than the number of atoms of the molecule.",
                )
            end
            c += atom_contributions[:, it]
        end
    end
    return c
end

#
# If a list of atom names is provided
#
function contrib(s::Selection, atom_contributions::Array{Float64}, names::Vector{String})
    indexes = Vector{Int}(undef, 0)
    for name in names
        index = findall(isequal(name), s.names)
        if length(index) == 0
            error(" Atom in input list is not part of solvent (or solute): $name")
        end
        append!(indexes, index)
    end
    return contrib(s, atom_contributions, indexes)
end

#
# If a list of atoms of PDBTools.Atom is provided
#
function contrib(
    s::Selection,
    atom_contributions::Array{Float64},
    atoms::Vector{PDBTools.Atom};
    warning = true,
)
    (warning && s.nmols > 1) && warning_nmols_types()
    indexes = [atom.index for atom in atoms]
    # Check which types of atoms belong to this selection
    selected_types = which_types(s, indexes, warning = warning)
    return contrib(s, atom_contributions, selected_types)
end

#
# If a residue of type PDBTools.Residue is provided
#
function contrib(
    s::Selection,
    atom_contributions::Array{Float64},
    residue::Residue;
    warning = true,
)
    (warning && s.nmols > 1) && warning_nmols_types()
    indexes = collect(residue.range)
    # Check which types of atoms belong to this selection
    selected_types = which_types(s, indexes, warning = warning)
    return contrib(s, atom_contributions, selected_types)
end

function warning_nmols_types()
    println(
        "WARNING: there is more than one molecule in this selection. " *
        "Contributions are summed over all atoms of the same type.",
    )
end


text:
# Tests SerializedMsgFormat in async mode operation.
# Both server and client run in the same process.
# In this mode the transport can be in-memory.
include("srvr.jl")
include("clnt.jl")

function test_serialized_msgformat()
    run_srvr(JuliaWebAPI.SerializedMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), true)
    run_clnt(JuliaWebAPI.SerializedMsgFormat(), JuliaWebAPI.ZMQTransport("127.0.0.1", 9999, ZMQ.REQ, false))
end

function test_inproc_transport()
    run_srvr(JuliaWebAPI.DictMsgFormat(), JuliaWebAPI.InProcTransport(:juliawebapi), true)
    run_clnt(JuliaWebAPI.DictMsgFormat(), JuliaWebAPI.InProcTransport(:juliawebapi))
end

# run tests if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runsermsgformat") && test_serialized_msgformat()
!isempty(ARGS) && (ARGS[1] == "--runinproctransport") && test_inproc_transport()


text:
using CitationNetworks
using LightGraphs#master
using LightGraphs.SimpleGraphs
using DataFrames
using SimpleWeightedGraphs
using MetaGraphs, GraphDataFrameBridge
using Test

A = [
 0  0  1  0  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0  0  0
 0  0  0  0  1  0  1  0  0  0  0
 0  0  0  0  0  1  0  0  0  1  1
 0  0  0  0  0  0  0  1  1  0  0
 0  0  0  0  0  0  0  0  1  1  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
]

g = SimpleDiGraph(A)

ew_spc, vw_spc, tf_spc = weights_spc(g, normalize = false)

# MetaDiGraph

eldf = DataFrame(
  from = [src(e) for e in edges(g)],
  to = [dst(e) for e in edges(g)],
  weight = ew_spc
)

mg = MetaDiGraph(eldf, :from, :to, weight = :weight)

@testset "SPC weights" begin
    spc_true = [3, 3, 4, 4, 2, 2, 1, 1, 2, 2, 1, 1]

    @test all(ew_spc .== spc_true)
    #@test all((Matrix(weights(mg)) .> 0) .== (A .> 0)) # check default weight in MetaGraphs
    @test tf_spc == 10.0
end

@testset "local forward main path" begin
    flmp1 = mainpath(mg, ForwardLocal())
    flmp2 = mainpath(mg, 1, ForwardLocal())
    flmp3 = mainpath(mg, [1,2], weights(g), ForwardLocal())

    set_mp_prop!(mg, flmp1, :flmp)

    flmp_true = [2, 4, 6, 9, 10]
    flmp_weights_true = [4, 2, 1, 1]
    edges_on_flmp_true = [false, false, true, false, false, true, false, false, false, false, true, true]

    @test [get_prop(mg, e, :flmp) for e in edges(mg)] == edges_on_flmp_true
    @test flmp1.weights == flmp_weights_true
    @test collect(filter_vertices(mg, :flmp, true)) == flmp_true
    @test mainpath(mg, 1, ForwardLocal()).weights == [3, 4, 2, 2]
end

@testset "local backward main path" begin

    blmp1 = mainpath(mg, BackwardLocal())
    blmp2 = mainpath(mg, 10, BackwardLocal())
    blmp3 = mainpath(mg, [10, 7], weights(mg), BackwardLocal())

    set_mp_prop!(mg, blmp1, :blmp)

    blmp_true = [1, 2, 3, 5, 7, 8, 9]
    blmp_weights_true = [2, 2, 2, 3, 3, 4]
    edges_on_blmp_true = [true, true, false, true, true, false, false, false, true, true, false, false]

    @test [get_prop(mg, e, :blmp) for e in edges(mg)] == edges_on_blmp_true
    @test blmp1.weights == blmp_weights_true
    @test collect(filter_vertices(mg, :blmp, true)) == blmp_true
    @test mainpath(mg, 7, BackwardLocal()).weights == [2, 3, 3]
end

@testset "standard global main path" begin

    sgmp1 = mainpath(mg, StandardGlobal())

    set_mp_prop!(mg, sgmp1, :sgmp)

    sgmp_true = [1, 2, 3, 5, 8, 9]
    sgmp_weights_true = [2, 2, 4, 3, 3]
    edges_on_sgmp_true = [true, true, false, true, false, false, false, false, true, true, false, false]

    @test [get_prop(mg, e, :sgmp) for e in edges(mg)] == edges_on_sgmp_true
    @test sgmp1.weights == sgmp_weights_true
    @test collect(filter_vertices(mg, :sgmp, true)) == sgmp_true
end


text:
using AbstractTrees

abstract type FSNode end

struct File <: FSNode
    path::String
end

AbstractTrees.children(::File) = ()

struct Directory <: FSNode
    path::String
end

function AbstractTrees.children(d::Directory)
    contents = readdir(d.path)
    children = Vector{Union{Directory,File}}(undef,length(contents))
    for (i,c) in enumerate(contents)
        path = joinpath(d.path,c)
        children[i] = isdir(path) ? Directory(path) : File(path)
    end
    return children
end

# in real life this probably wouldn't be useful, but it's convenient for testing
AbstractTrees.nodevalue(n::FSNode) = n.path

AbstractTrees.printnode(io::IO, d::Directory) = print(io, basename(d.path))
AbstractTrees.printnode(io::IO, f::File) = print(io, basename(f.path))

function mk_tree_test_dir(f, parentdir=tempdir(); prefix="jl_")
    mktempdir(parentdir; prefix) do path
        cd(path) do
            open(io -> write(io, "test1"), "f1"; write=true)
            mkdir("A")
            open(io -> write(io, "test2"), joinpath("A","f2"); write=true)
            mkdir("B")
            f(path)
        end
    end
end


text:
function zad31a(n,k)
    return faktorijel(n)/(faktorijel(k)*faktorijel(n-k))
end

function faktorijel(n)
    if n == 2
        return n
    end
    return n * fakt(n-1)
end

function zad31b(n,k)
    if n == k || k == 0
        return 1
    end
    return zad31b(n-1,k) + zad31b(n-1,k-1)
end

function zad31c(n,k)
    if k == 0
        return 1
    end
    return zad31c(n,n-1) * (n-k+1)/k

end


text:
const mu_dmu_dict = Dict(
    "uniform"=>[
        0.01068 0.01017 0.00877 0.00981 0.01028 0.0 0.00984 0.0099 0.00955 0.01017 0.01054;
        0.00987 0.00953 0.00898 0.00913 0.01002 0.0 0.01034 0.01027 0.00905 0.00942 0.01055;
        0.00963 0.00852 0.00801 0.0089 0.00909 0.00972 0.00865 0.00859 0.00869 0.00914 0.00934;
        0.01034 0.00943 0.00879 0.00893 0.00982 0.0 0.01003 0.00886 0.00828 0.00955 0.01016;
        0.01049 0.00923 0.00969 0.00962 0.01039 0.0 0.0099 0.00877 0.00956 0.00937 0.00998;
        0.0 0.01101 0.0 0.0 0.0 0.0 0.01028 0.00981 0.00907 0.00994 0.01056;
        0.0106 0.00949 0.00964 0.0101 0.01029 0.0 0.0 0.0 0.01042 0.0 0.0;
        0.00923 0.00924 0.00817 0.00923 0.0097 0.0 0.01026 0.0098 0.00868 0.01029 0.01074;
        0.0093 0.00913 0.00834 0.00849 0.00957 0.0 0.00946 0.00925 0.00858 0.00889 0.00994;
        0.01005 0.00954 0.00918 0.00987 0.00904 0.01095 0.00953 0.00873 0.00938 0.00917 0.01026;
        0.01039 0.01007 0.0099 0.01002 0.01025 0.0 0.01013 0.00991 0.0093 0.00985 0.01064],
    "random_state_variant"=>[
        0.0731 0.06118 0.04944 0.04841 0.06033 0.0 0.02796 0.02958 0.03639 0.03986 0.03665;
        0.02722 0.02042 0.02538 0.02442 0.01502 0.0 0.02416 0.02196 0.01886 0.02062 0.02693;
        0.00801 0.01088 0.01536 0.01625 0.01471 0.01353 0.0154 0.01263 0.00435 0.00436 0.00808;
        0.00538 0.00642 0.0103 0.01385 0.01405 0.0 0.00623 0.00492 0.00236 0.00263 0.00168;
        0.00292 0.00392 0.00934 0.01188 0.01428 0.0 0.00087 0.00121 0.00183 0.0022 0.00178;
        0.0 0.00074 0.0 0.0 0.0 0.0 0.00055 0.00107 0.00148 0.00202 0.00211;
        0.00151 0.00091 0.00146 0.00235 0.0032 0.0 0.0 0.0 0.00116 0.0 0.0;
        0.00111 0.00125 0.00146 0.00161 0.00111 0.0 0.00213 0.00212 0.00179 0.00262 0.00299;
        0.00057 0.0007 0.00094 0.00082 0.00099 0.0 0.00184 0.00194 0.00229 0.00215 0.00249;
        0.00069 0.00056 0.00054 0.00059 0.00073 0.00141 0.00153 0.00161 0.0012 0.00157 0.00154;
        0.00031 0.00055 0.00059 0.00079 0.00092 0.0 0.00127 0.00077 0.00067 0.00068 0.0005],
    "favor_down"=>[
        0.09853  0.09109  0.08863  0.09136  0.09845  0.0      0.095    0.0893   0.08658  0.08963  0.09608;
        0.00731  0.00712  0.00639  0.0067   0.00764  0.0      0.00719  0.00679  0.00667  0.00699  0.00676;
        0.00054  0.00035  0.00046  0.00041  0.00042  0.00043  0.00055  0.00049  0.00045  0.00048  0.00053;
        4.0e-5   6.0e-5   7.0e-5   5.0e-5   5.0e-5   0.0      2.0e-5   5.0e-5   6.0e-5   3.0e-5   5.0e-5;
        2.0e-5   1.0e-5   0.0      1.0e-5   0.0      0.0      0.0      0.0      1.0e-5   1.0e-5   1.0e-5;
        0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      1.0e-5   0.0;
        0.0      2.0e-5   4.0e-5   2.0e-5   2.0e-5   0.0      0.0      0.0      1.0e-5   0.0      0.0;
        0.0      0.0      0.0      1.0e-5   0.0      0.0      0.0      0.0      0.0      0.0      0.0;
        0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0;
        0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0;
        0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0      0.0
    ],
    "corners"=>[
        0.0105   0.01014  0.00938  0.01074  0.01053  0.0      0.01033  0.00939  0.00931  0.00997  0.01063;
        0.00964  0.00976  0.00898  0.00882  0.01014  0.0      0.01051  0.00866  0.00893  0.00933  0.00965;
        0.00943  0.0092   0.00829  0.00898  0.00919  0.01025  0.00878  0.00921  0.00842  0.0091   0.00937;
        0.00958  0.00883  0.00877  0.00959  0.01033  0.0      0.0095   0.00897  0.00847  0.00881  0.00973;
        0.01015  0.00913  0.00941  0.01018  0.01122  0.0      0.01021  0.00961  0.00887  0.00919  0.00985;
        0.0      0.01096  0.0      0.0      0.0      0.0      0.01077  0.00971  0.00938  0.00951  0.01048;
        0.01027  0.00919  0.00942  0.00946  0.01037  0.0      0.0      0.0      0.01045  0.0      0.0;
        0.00935  0.00907  0.00863  0.00934  0.01012  0.0      0.01044  0.00994  0.00846  0.00992  0.01077;
        0.0089   0.00901  0.00792  0.00849  0.00922  0.0      0.01012  0.00911  0.00897  0.00959  0.00979;
        0.00952  0.00992  0.00936  0.00925  0.00942  0.01066  0.00913  0.00917  0.00908  0.01018  0.01046;
        0.0109   0.00944  0.00974  0.01007  0.01044  0.0      0.01028  0.01003  0.01017  0.00946  0.01023;
    ],
    "random_state_weight_variant"=>[
        0.00041  0.00577  0.00095  0.00069  0.00026  0.0      0.02119  0.02458  0.0358   0.037    0.02396;
        0.00019  0.00021  0.0008   0.00125  0.0014   0.0      0.01625  0.0125   0.00618  0.00641  0.008;
        0.00019  0.00108  0.00142  0.00146  0.00109  0.00223  0.01206  0.01124  0.0023   0.00296  0.00255;
        0.0034   0.01087  0.00341  0.00189  0.00177  0.0      0.10153  0.0996   0.00163  0.00399  0.00424;
        0.01356  0.00385  0.00166  0.00197  0.00215  0.0      0.00694  0.00552  0.00366  0.00812  0.00888;
        0.0      0.00045  0.0      0.0      0.0      0.0      0.03168  0.00552  0.00392  0.01613  0.01213;
        0.00068  0.00156  0.00171  0.00236  0.11574  0.0      0.0      0.0      0.00146  0.0      0.0;
        0.00183  0.00241  0.00199  0.00109  0.00074  0.0      0.00229  0.00219  0.00219  0.00476  0.00675;
        0.01277  0.00538  0.00238  0.00169  0.00196  0.0      0.00173  0.00169  0.00193  0.00479  0.00835;
        0.00277  0.00289  0.00356  0.00533  0.00318  0.00307  0.00201  0.00092  0.0018   0.00549  0.08244;
        0.00113  0.00269  0.00384  0.00472  0.0049   0.0      0.0033   0.05392  0.00326  0.00579  0.00972;
    ]

)


text:
@kwargs function mul!(
    C::GBMatrix,
    A::GBMatOrTrans,
    B::GBMatOrTrans;
    op::SemiringUnion = Semirings.PLUS_TIMES_SEMIRING
)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    op = getoperator(op, eltype(C))
    A, desc, B = _handletranspose
    libgb.GrB_mxm(C, mask, accum, op, A, B, desc)
end

@kwargs function mul!(
    w::GBVector,
    u::GBVector,
    A::GBMatOrTrans;
    op::SemiringUnion = Semirings.PLUS_TIMES_SEMIRING
)
    size(u) == size(A, 1) || throw(DimensionMismatch())
    size(w) == size(A, 2) || throw(DimensionMismatch())
    op = getoperator(op, eltype(w))
    u, desc, A = _handletranspose
    libgb.GrB_vxm(w, mask, accum, op, u, A, desc)
end

@kwargs function mul!(
    w::GBVector,
    A::GBMatOrTrans,
    u::GBVector;
    op::SemiringUnion = Semirings.PLUS_TIMES_SEMIRING
)
    size(u) == size(A, 2) || throw(DimensionMismatch())
    size(w) == size(A, 1) || throw(DimensionMismatch())
    op = getoperator(op, eltype(w))
    A, desc, u = _handletranspose
    libgb.GrB_mxv(w, mask, accum, op, A, u, desc)
end

@kwargs function mul(
    A::GBVecMatOrTrans,
    B::GBVecMatOrTrans;
    op::SemiringUnion = Semirings.PLUS_TIMES_SEMIRING
)
    op = getoperator(op, optype(eltype(A), eltype(B)))
    t = juliatype(ztype(op))
    if A isa GBVector && B isa GBMatOrTrans
        C = GBVector{t}(size(B, 2))
    elseif A isa GBMatOrTrans && B isa GBVector
        C = GBVector{t}(size(A, 1))
    elseif A isa GBMatOrTrans && B isa GBMatOrTrans
        C = GBMatrix{t}(size(A, 1), size(B, 2))
    end
    A, desc, B = _handletranspose(A, desc, B)
    mul!(C, A, B; op, mask, accum, desc)
    return C
end


text:
cd(@__DIR__); runapprox = true; allNGWPs = false; include("setups/sunflower.jl")
gr(dpi = 200)

##
VM_NGWP = vm_ngwp(𝚽, GP_dual)
f = matread("../datasets/sunflower_barbara_voronoi.mat")["f_trouser_voronoi"]
G_Sig.f = reshape(f, (N, 1))
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
println("================ftrouser-VM-NGWP-top-basis-vectors=================")
for i in 2:17
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual, dr, dc)
    print("(j, k, l) = ($(j), $(k), $(l))  ")
    if j == jmax
        print("φ_{$(GP_dual.ind[dr]-1)}")
    end
    println()
    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
                cbar = false, clims = (-0.15, 0.15))
    savefig(plt,
        "../figs/SunFlower_ftrouser_DAG_VM_NGWP_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
using HypothesisTests, Base.Test
using HypothesisTests: default_tail

sim_data_h0=[
    0.297287984535462;0.382395967790608;-0.597634476728231;-0.0104452446373756;
    -0.839026854388764;0.311111338498334;2.29508782383731;-2.26708634880053;
    0.529965576166746;0.431421526422912;0.583708287568779;0.963271605038191;
    0.458790955053717;-0.522336757421508;0.408395838324752;-0.0504512299336653;
    -0.693653643803886;-1.77383325340806;0.120765964937439;-0.757633037720528;
    -1.72975618139069;0.795948622004616;0.670061981256062;0.550852377500761;
    -0.0633745924295606;1.33694278225092;-0.0731486333377393;-0.745464316640756;
    -1.22005512853465;-0.0531773375276925;-0.165136357802834;-2.11537024842731;
    -0.0667679865065892;1.22246184846806;0.567695459738444;0.763498928323387;
    0.378598879853350;-0.645597056636597;-0.664646244391031;-1.80302632421089;
    -1.28904760384974;-0.335270079518118;0.0704676013611580;0.341793767379823;
    1.73516566605438;1.29991596431019;0.206364081402617;-1.00885995553760;
    -0.850055670315398;1.12940943588433;-1.34853345649183;1.24029090042906;
    -0.0556470932073243;0.759390213945460;-0.0302668866467625;-0.361082667756094;
    -2.00730320253256;-0.568859127346638;-1.14962752439700;1.88062217463119;
    -1.47062454138738;-0.535210970165732;-0.963543959825501;-1.38510806767591;
    0.134300237561270;-0.616116666008819;-1.71998935631974;0.320769279987411;
    -1.44737289788151;-0.492271188069566;-0.413087383983183;1.00192851245595;
    1.17468132597511;-1.32697034506205;-1.59321486804407;-0.438709975758797;
    1.22579749799787;1.48909027806961;2.15943318710085;-0.909479018037503;
    1.14782030061250;-0.505072196522373;-0.267067035218130;1.49908239393073;
    0.797304160175390;-0.171640443307097;-0.469079609930096;0.217624021208600;
    0.359146407602350;0.320024715708058;0.259216328055686;0.459695985041628;
    0.209986360042143;1.52317539117726;1.25979460361220;-0.165283101114261;
    0.156194494887636;-0.814315294376699
    ]

t = HypothesisTests.BoxPierceTest(sim_data_h0,2,1)

@test t.n == 98
@test t.lag == 2
@test t.dof == 1
@test t.Q ≈ 1.233942980734545
@test pvalue(t) ≈ 0.2666415904008932
@test default_tail(t) == :right
show(IOBuffer(), t)

t = HypothesisTests.LjungBoxTest(sim_data_h0,5,2)

@test t.n == 98
@test t.lag == 5
@test t.dof == 2
@test t.Q ≈ 3.2090126519163626
@test pvalue(t) ≈ 0.36050846449240337
@test default_tail(t) == :right
show(IOBuffer(), t)

sim_data_h1 = [
    0.297287984535462;0.739141549233162;0.200148986990924;0.00799107498178517;
    -0.889482260507899;-0.758664696605681;1.65153486606286;-0.0576451005433905;
    -0.0346690043041814;0.407112251420912;1.08264369056513;2.14031035829007;
    2.70237027783226;2.07841446849019;2.09178211716330;1.83616297011523;
    0.882207285185406;-1.26603340222014;-1.66313630328235;-2.37358658099330;
    -4.07911918759795;-3.38691842881493;-2.17050437704247;-1.03767734630572;
    -0.657436094883687;0.859322672282215;1.15526940186602;0.383062163913809;
    -1.10696135239789;-1.49644960957930;-1.62878748357863;-3.62098034584788;
    -3.92330815645045;-2.39921383551812;-1.13436869594816;0.122020643841029;
    0.865334261247033;0.356197863707534;-0.496809086316101;-2.50605658690248;
    -4.14727278223788;-4.56018044213284;-4.15756709452688;-3.27923261341258;
    -0.952643341682653;1.14051373831478;1.86077356988515;0.881914206830145;
    -0.349990693084770;0.444846342133559;-0.709720638006125;0.255172232181642;
    0.463475776812484;1.23900947646595;1.31750175206863;0.848216591786478;
    -1.38469381800938;-2.48495668649384;-3.71616740278679;-1.83329170276481;
    -2.55572436386911;-3.05209269597923;-3.85933788583984;-5.10068572188995;
    -4.82872126295472;-4.88037646498750;-6.12782473541832;-5.56850746301833;
    -6.29123443287801;-6.37120026861768;-6.17115737646099;-4.49210025871193;
    -2.36449177154091;-2.81673039329756;-4.26394380853887;-4.71042342801617;
    -3.14752747305988;-0.874815661197385;2.05391263558195;1.81766084302006;
    2.71283952156198;2.20503697644599;1.56512548004846;2.71572187705509;
    3.58663276862696;3.31760231592873;2.43605333859629;2.14560733274553;
    2.20305920531810;2.32001356226611;2.38231484117959;2.62246972577731;
    2.64225557862103;3.90714116778931;5.15568733137306;4.84939934619662;
    4.42876751091166;3.04538591485831
    ]

t = HypothesisTests.BoxPierceTest(sim_data_h1,3)

@test t.n == 98
@test t.lag == 3
@test t.dof == 0
@test t.Q ≈ 176.16899390632153
@test pvalue(t) ≈ 5.925682677971866e-38
show(IOBuffer(), t)

t = HypothesisTests.LjungBoxTest(sim_data_h1,12)

@test t.n == 98
@test t.lag == 12
@test t.dof == 0
@test t.Q ≈ 271.8894095341075
@test pvalue(t) ≈ 3.6622231247462687e-51
show(IOBuffer(), t)


text:
"""
        waspas(decisionMat, weights, fns, lambda)

Apply WASPAS (Weighted Aggregated Sum Product ASsessment ) for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alterntives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.
 - `lambda::Float64`: joint criterion. 0<=lambda<=1, default=0.5.

# Description
waspas() applies the WASPAS method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::WASPASResult`: WASPASResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [3        12.5        2        120        14        3;
       5        15        3        110        38        4;
       3        13        2        120        19        3;
       4        14        2        100        31        4;
       3        15        1.5        125        40        4]
5×6 Array{Float64,2}:
 3.0  12.5  2.0  120.0  14.0  3.0
 5.0  15.0  3.0  110.0  38.0  4.0
 3.0  13.0  2.0  120.0  19.0  3.0
 4.0  14.0  2.0  100.0  31.0  4.0
 3.0  15.0  1.5  125.0  40.0  4.0

julia> df = makeDecisionMatrix(decmat)
5×6 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6
     │ Float64  Float64  Float64  Float64  Float64  Float64
─────┼──────────────────────────────────────────────────────
   1 │     3.0     12.5      2.0    120.0     14.0      3.0
   2 │     5.0     15.0      3.0    110.0     38.0      4.0
   3 │     3.0     13.0      2.0    120.0     19.0      3.0
   4 │     4.0     14.0      2.0    100.0     31.0      4.0
   5 │     3.0     15.0      1.5    125.0     40.0      4.0

julia> weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201];

julia> fns = [maximum, minimum, minimum, maximum, minimum, maximum];

julia> lambda = 0.5;

julia> result = wpm(df, weights, fns, lambda);

julia> result.scores
5-element Array{Float64,1}:
 0.8050212165665626
 0.7750597051081832
 0.770180748518019
 0.7964243424353943
 0.7882389370890854

 julia> result.bestIndex
 1
```

# References
Zavadskas, E. K., Turskis, Z., Antucheviciene, J., & Zakarevicius, A. (2012). Optimization of Weighted Aggregated Sum Product Assessment. Elektronika Ir Elektrotechnika, 122(6), 3-6. https://doi.org/10.5755/j01.eee.122.6.1810
Aytaç Adalı, E. & Tuş Işık, A.. (2017). Bir Tedarikçi Seçim Problemi İçin SWARA ve WASPAS Yöntemlerine Dayanan Karar Verme Yaklaşımı. International Review of Economics and Management, 5 (4) , 56-77. DOI: 10.18825/iremjournal.335408
"""
function waspas(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1}, lambda::Float64=0.5):WASPASResult

    row, col = size(decisionMat)
    normalizedDecisionMat = similar(decisionMat)
    w = unitize(weights)
    colminmax = zeros(Float64, col)
    @inbounds for i in 1:col
        colminmax[i] = decisionMat[:, i] |> fns[i]
        if fns[i] == maximum
            normalizedDecisionMat[:, i] = decisionMat[:, i] ./ colminmax[i]
        elseif fns[i] == minimum
            normalizedDecisionMat[:, i] = colminmax[i] ./ decisionMat[:, i]
    end
    end
    scoreMat = similar(normalizedDecisionMat)
    for i in 1:col
        scoreMat[:, i] = normalizedDecisionMat[:, i].^w[i]
    end

    scoresWPM = zeros(Float64, row)
    for i in 1:row
        scoresWPM[i] = prod(scoreMat[i, :])
    end

    scoresWSM = w * normalizedDecisionMat |> rowsums

    scoreTable = [scoresWSM  scoresWPM]

    l = unitize([lambda, 1 - lambda])

    scores = l' .*  scoreTable |> rowsums

    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = WASPASResult(
        decisionMat,
        normalizedDecisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
module TestLogging
using Test
using Pia.Logging

@testset "Logging" begin
    @testset "basic" begin
        let buffer = IOBuffer()
            @test Logging.log(1, 2, 3; duplex=buffer) === nothing
            @test length(Logging.Uncategorized.messages) == 1
            seek(buffer, 0)
            @test Logging.Uncategorized.messages[1].message == "1 2 3"
            @test read(buffer, String) == "1 2 3\n"

            take!(buffer)
            seek(buffer, 0)
            @test Logging.log(1, 2, 3; duplex=buffer, joiner=", ") === nothing
            @test length(Logging.Uncategorized.messages) == 2
            seek(buffer, 0)
            @test Logging.Uncategorized.messages[2].message == "1, 2, 3"
            @test read(buffer, String) == "1, 2, 3\n"
        end
    end

    @testset "custom category" begin
        let cat = LogCategory("TestCategory", nothing)
            @test cat(1, 2, 3) === nothing
            @test Logging.log(cat, 1, 2, 3; joiner=", ") === nothing
            @test length(cat.messages) == 2 && cat.messages[1].message == "1 2 3" && cat.messages[2].message == "1, 2, 3"
        end

        let buffer = IOBuffer(), cat = LogCategory("TestCategory", buffer)
            @test cat(1, 2, 3) === nothing
            @test buffer.size > 0
            seek(buffer, 0)
            @test read(buffer, String) == "1 2 3\n"
        end
    end

    @testset "purge" begin
        let cat = LogCategory("TestCategory", nothing)
            cat(1, 2, 3)
            @assert length(cat.messages) == 1
            @test empty!(cat) === cat
            @test length(cat.messages) == 0
        end

        let cat = LogCategory("TestCategory", nothing)
            cat(1, 2, 3)
            sleep(0.1)
            t0 = time()
            cat(2, 3, 4)
            @assert length(cat.messages) == 2
            @test empty!(cat, older_than=t0) === cat
            @test length(cat.messages) == 1 && cat.messages[1].message == "2 3 4"
        end
    end
end

end # module TestLogging


text:
using Documenter
using FlightMechanics

makedocs(authors="Alejandro Saez Mollejo",
         sitename = "FlightMechanics.jl",
         pages = [
         "Home" => "index.md",
         "API Reference" => ["api-reference.md"],
         "Basics" => [
           "isa.md",
           "coordinate-systems.md",
           "dynamic_models.md",
           "anemometry.md",
           "kinematics.md"
         ],
         "Simulator" => [
           "models.md",
           "aircrafts.md",
           "trimmer.md"
         ]
         ])

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/AlexS12/FlightMechanics.jl.git",
    julia = "0.7"
)


text:
function interchain_stats(chains)

    d = Dict()

    for param in chains.name_map.parameters
        μ = std(mean(chains[param], dims=1))
        σ = std(std(chains[param], dims=1))
        push!(d, param => Dict(:μ => μ, :σ => σ))
    end

    return d

end


text:
module Utils

"""
    apply_chmod_recursively(root_dir::String; dir_mode::Integer, file_mode::Integer)
"""
function apply_chmod_recursively(root_dir::String;
                                 dir_mode::Integer,
                                 file_mode::Integer)
    chmod(root_dir, dir_mode)
    for (root, dirs, files) in walkdir(root_dir)
        for dir in dirs
            full_dir_path = joinpath(root, dir)
            isdir(full_dir_path) || throw(ArgumentError("unexpected error"))
            chmod(full_dir_path, dir_mode)
        end
        for file in files
            full_file_path = joinpath(root, file)
            isfile(full_file_path) || throw(ArgumentError("unexpected error"))
            chmod(full_file_path, file_mode)
        end
    end
    return nothing
end

"""
    make_world_readable_recursively(root_dir::String)
"""
function make_world_readable_recursively(root_dir::String)
    return apply_chmod_recursively(root_dir; dir_mode=0o555, file_mode=0o444)
end

"""
    make_world_writeable_recursively(root_dir::String)
"""
function make_world_writeable_recursively(root_dir::String)
    return apply_chmod_recursively(root_dir; dir_mode=0o777, file_mode=0o666)
end

end # module


text:
# This is the first Stan model in Statistical Rethinking Edition 1 (page 249).
# In Edition 2 (page 242)

# \toc

# ## Data

import CSV
import Random

using DataFrames
using Turing
using TuringModels

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "rugged.csv")
df = CSV.read(data_path, DataFrame)

df.log_gdp = log.(df.rgdppc_2000)
dropmissing!(df)

df = select(df, :log_gdp, :rugged, :cont_africa);

# DataFrame `df` is shown Section [df](#df).

# ## Model

@model function model_fn(y, x₁, x₂)
  σ ~ truncated(Cauchy(0, 2), 0, Inf)
  βAR ~ Normal(0, 10)
  βR ~ Normal(0, 10)
  βA ~ Normal(0, 10)
  α ~ Normal(0, 100)

  μ = α .+ βR * x₁ .+ βA * x₂ .+ βAR * x₁ .* x₂
  y ~ MvNormal(μ, σ)
end

model = model_fn(df.log_gdp, df.rugged, df.cont_africa);

# ## Output

chns = sample(model, NUTS(), 1000)

# \defaultoutput{}

# ## Original output

"""
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean         SD       Naive SE       MCSE      ESS
    a  9.22360053 0.139119116 0.0021996664 0.0034632816 1000
   bR -0.20196346 0.076106388 0.0012033477 0.0018370185 1000
   bA -1.94430980 0.227080488 0.0035904578 0.0057840746 1000
  bAR  0.39071684 0.131889143 0.0020853505 0.0032749642 1000
sigma  0.95036370 0.052161768 0.0008247500 0.0009204073 1000

Quantiles:
          2.5%       25.0%       50.0%      75.0%        97.5%
    a  8.95307475  9.12719750  9.2237750  9.31974000  9.490234250
   bR -0.35217930 -0.25334425 -0.2012855 -0.15124725 -0.054216855
   bA -2.39010825 -2.09894500 -1.9432550 -1.78643000 -1.513974250
  bAR  0.13496995  0.30095575  0.3916590  0.47887625  0.650244475
sigma  0.85376115  0.91363250  0.9484920  0.98405750  1.058573750
""";

# ## df

df


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -6.563131467197746 <= q <= 9.310824610097802)

                     add_NL_constraint(m, :(log(1 + exp(0.5323461867533106 + 0.6233251355511138*log(1 + exp(0.9117705078314118 + -0.4325489455874667*$(x[1]) + 0.44461795780922087*$(x[2]) + -0.310674968881143*$(x[3]) + 0.3874148238299093*$(x[4]))) + 0.13135733952437967*log(1 + exp(0.4793544889478012 + 0.21125743039643652*$(x[1]) + -0.625324840825729*$(x[2]) + -0.9496988785997882*$(x[3]) + 0.8662818947925528*$(x[4]))) + -0.0038770043009219712*log(1 + exp(0.7417071748063861 + 0.06375566706279434*$(x[1]) + -0.03899151888275032*$(x[2]) + 0.7226457302550271*$(x[3]) + 0.006005128865761655*$(x[4]))))) + log(1 + exp(0.48577545639109676 + -0.6684827900810446*log(1 + exp(0.9117705078314118 + -0.4325489455874667*$(x[1]) + 0.44461795780922087*$(x[2]) + -0.310674968881143*$(x[3]) + 0.3874148238299093*$(x[4]))) + 0.9262921195438665*log(1 + exp(0.4793544889478012 + 0.21125743039643652*$(x[1]) + -0.625324840825729*$(x[2]) + -0.9496988785997882*$(x[3]) + 0.8662818947925528*$(x[4]))) + 0.29789654203123295*log(1 + exp(0.7417071748063861 + 0.06375566706279434*$(x[1]) + -0.03899151888275032*$(x[2]) + 0.7226457302550271*$(x[3]) + 0.006005128865761655*$(x[4]))))) + log(1 + exp(0.2863540742777313 + -0.41164742912319063*log(1 + exp(0.9117705078314118 + -0.4325489455874667*$(x[1]) + 0.44461795780922087*$(x[2]) + -0.310674968881143*$(x[3]) + 0.3874148238299093*$(x[4]))) + 0.6047572744817185*log(1 + exp(0.4793544889478012 + 0.21125743039643652*$(x[1]) + -0.625324840825729*$(x[2]) + -0.9496988785997882*$(x[3]) + 0.8662818947925528*$(x[4]))) + -0.7133357903136108*log(1 + exp(0.7417071748063861 + 0.06375566706279434*$(x[1]) + -0.03899151888275032*$(x[2]) + 0.7226457302550271*$(x[3]) + 0.006005128865761655*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions