text:
using MLJLinearModels, StableRNGs
using BenchmarkTools, Random, LinearAlgebra
DO_COMPARISONS = false; include("../testutils.jl")

n, p = 50_000, 500
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)

# =============== #
# RIDGE FUNCTIONS #
#
# Hv!           ✅
# =============== #

# No fit_intercept
s = R.scratch(X; i=false)
λ = 0.5
ridge = RidgeRegression(λ; fit_intercept=false)
Hv! = R.Hv!(ridge, X, y, s)
v   = randn(p)
Hv  = similar(v)

# Sept 13, 2019 :: 15.300 ms (0 allocations: 0 bytes)
# May 12, 2020 :: 15.785 ms (22 allocations: 432 bytes)
#              :: 15.332 ms (22 allocations: 432 bytes)
#              :: 15.198 ms (22 allocations: 432 bytes)
@btime Hv!($Hv, $v)

# With fit_intercept
s = R.scratch(X; i=true)
ridge = RidgeRegression(λ)
Hv! = R.Hv!(ridge, X, y, s)
v  = randn(p+1)
Hv = similar(v)

# Sept 13, 2019 :: 26.246 ms (5 allocations: 4.22 KiB)
# (decrease in perf due to views, and need a sum on an array)
# May 12, 2020 :: 120.892 ms (34 allocations: 4.77 KiB) NOTE != rng
#              :: 118.999 ms (35 allocations: 4.81 KiB) NOTE != rng
#              :: 114.046 ms (34 allocations: 768 bytes)
@btime Hv!($Hv, $v)

# =============== #
# LASSO FUNCTIONS #
#
# smooth_fg!
# =============== #

lasso = LassoRegression(0.5, fit_intercept=false)
s = R.scratch(X; i=false)
smooth_fg! = R.smooth_fg!(lasso, X, y, s)
v = randn(p)
g = similar(v)

# Sept 13, 2019 :: 15.508 ms (3 allocations: 390.72 KiB)
# May 12, 2020 :: 16.659 ms (32 allocations: 391.34 KiB)
#              :: 15.868 ms (32 allocations: 391.34 KiB)
@btime smooth_fg!($g, $v);

((X, _, _), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.1)

λ  = 50
lr = LassoRegression(λ)
J  = objective(lr, X, y1)

fista = FISTA()
ista = ISTA()

# <old> initial: ~650ms ; 1725 allocations;  104.24 MiB
# using cache: ; ~650ms ; 1666 allocations;   95.02 MiB
# May 12, 2020 (note RNG is different)
# -- 780.603 ms (4482 allocations: 97.43 MiB)
# -- 772.273 ms (4488 allocations: 98.20 MiB)
@btime θ_fista = fit($lr, $X, $y1, solver=$fista)

# using cache ;  ~590ms ; 1584 allocations;   91.36 MiB
# May 12, 2020 (note RNG is different)
# -- 655.937 ms (4119 allocations: 91.63 MiB)
# -- 644.725 ms (4125 allocations: 92.40 MiB)
@btime θ_ista  = fit($lr, $X, $y1, solver=$ista)

rr = RidgeRegression(λ)
θ_ridge = fit(rr, X, y1)
cg = CG()

# <old> initial: ~ 220ms; 100 allocations; 3.5 MiB
# using cache:   ~ 216ms;  76 allocations; 458.5KiB
# May 12, 2020 :: 982.551 ms (496 allocations: 467.69 KiB)
#              :: 1.005 s (509 allocations: 1.22 MiB)
@btime θ_ridge_cg = fit($rr, $X, $y1, solver=$cg)


text:
# The various columns in the Base_GD file, with rough English translations

# Year  Month  Energy generated?  Energy Desposited?  Withdrawn energy?     Amount of energy used     Company   Tariff number
# A       B          C                     D               E                         F                 G               H
# ANNO, MES,  ENERGIA_GENERADA,  ENERGIA_DEPOSITADA, ENERGIA_RETIRADA  IMPORTE_POR_ENERGIA_RETIRADA,  EMPRESA, CODIGO_TARIFA,

# Tariff description   Total consumption
#         I                 J                K                       L               M                 N               O
# CODIGO_TARIFA1,     TOTAL_CONSUMO_KWH, TOTAL_IMPORTE_ENE, Suma de DISTRITO, Recuento de DISTRITO,  Provincia1,  CONSUMO_NATURAL,  SECTOR,  Data_Check

using DataFrames
using PyPlot
using GLM

# Tariff types:
# 1: Residential
# 2: Residencial horaria
# 4: Commercial y services mon-mica
# 5: Comercios y servicios con potencia
# 6: Industrial mon-mica
# 7: Industrial con potencia
# 8: Preferencial mon-mica
# 9: Preferencial con potencia
# 10: Promocional mon—mica
# 11: Promocional con potencia
# 12: Media tensi—n a

tariff_to_name = Dict([
        (1, "Residential"),
        (2, "Residencial horaria"),
        (4, "Commercial y services mon-mica"),
        (5, "Comercios y servicios con potencia"),
        (6, "Industrial mon-mica"),
        (7, "Industrial con potencia"),
        (8, "Preferencial mon-mica"),
        (9, "Preferencial con potencia"),
        (10, "Promocional mon—mica"),
        (11, "Promocional con potencia"),
        (12, "Media tensi—n a")
        ])

tariff_categories = ["Residential", "Commercial Industrial", "Medium Voltage"]

# Map from individual tariffs to the tariff category they're within
tariff_category_mappings = Dict([
        (1, tariff_categories[1]),
        (4, tariff_categories[2]),
        (5, tariff_categories[2]),
        (6, tariff_categories[2]),
        (7, tariff_categories[2]),
        (12, tariff_categories[3])
        ])

# Obtain the PV output timeseries to use for all installation and modeling calculations
if !@isdefined pv_output
    if @isdefined(cnfl)
        pv_output = monte_carlo_solar_output(cnfl=cnfl)
    else
        pv_output = monte_carlo_solar_output()
    end
end

# Since the installed PV Systems are very unlikely to be fully efficient, we will add a loss factor of 30%
adjusted_pv_output = pv_output * 0.7

# Obtain total PV system output for each month, since utility bill data is monthly
output_by_month = Array{Float64}(undef,12)
days_per_month = [0,31,28,31,30,31,30,31,31,30,31,30,31] # First entry for looping purposes only
for i=1:12
    output_by_month[i] = sum(adjusted_pv_output[sum(days_per_month[1:i])*24+1:sum(days_per_month[1:i+1])*24])
end

# This function cleans up the data for a particular utility+tariff
function create_consumption_and_installation_arrays(utility_bill_df)
    # Only keep entries with valid ENERGIA_GENERADA and CONSUMO_NATURAL
    cleaned_df = filter(row -> (!ismissing(row.ENERGIA_GENERADA) && !ismissing(row.CONSUMO_NATURAL)), utility_bill_df)

    # Only keep entries with PV system output > 100 kWh per month (this would mean it's a roughly 1kW system)
    string_to_float(str) = tryparse(Float64, str)
    cleaned_df[:CONSUMO_NATURAL] = map(string_to_float, cleaned_df[:CONSUMO_NATURAL])
    cleaned_df = filter(row -> (row.CONSUMO_NATURAL > 100), cleaned_df)

    # Construct consumption and installation arrays
    # Only keep values where installation is >= 1
    consumption_array = Array{Float64}(undef,0)
    installation_array = Array{Float64}(undef,0)
    for r in eachrow(cleaned_df)
        month = r.MES
        generation = r.ENERGIA_GENERADA
        installation = generation / output_by_month[month]
        if installation >= 1
            push!(installation_array, installation)
            push!(consumption_array, r.CONSUMO_NATURAL)
        end
    end

    return [consumption_array, installation_array]
end

# This function takes in consumption and installation arrays, and plots them
function plot_consumption_and_installation(consumption_and_installation, tariff_name, company_name)
    plt.figure()
    # Plot out the actual PV system sizes that people have based on their energy bills
    the_plot = plot(consumption_and_installation[1], consumption_and_installation[2], ".")
    ylabel("PV System Capacity [kW]")
    xlabel("Consumer Monthly Energy use [kWh]")
    grid("on");
    title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_name))
    return the_plot
end

# Function plotting the estimated “true” data versus the “economically rational” model
function plot_base_GD_vs_economically_rational(consumption_and_installation, tariff_name, company_name, consumption, model_predictions; model_descriptions=[], x_max=nothing, y_max=nothing, linreg=false)
    # Plot out the actual PV system sizes that people have based on their energy bills
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.scatter(consumption_and_installation[1], consumption_and_installation[2], marker=".", c = "#ffdf22", label = string("Actual PV System Installation for ", tariff_name, " Consumer"))

    # If desired, plot a linear regression of that data
    if (linreg)
        cons_and_install = DataFrame()
        cons_and_install[:CONS] = consumption_and_installation[1]
        cons_and_install[:INSTALL] = consumption_and_installation[2]
        data_b, data_m = coef(lm(@formula(INSTALL ~ CONS), cons_and_install))
        xlims = ax1.get_xlim()
        x_vals = xlims[1]:(xlims[2]-xlims[1])/100:xlims[2]
        y_vals = data_m * x_vals .+ data_b
        ax1.plot(x_vals, y_vals, "--", c="k", label = string("Least-squares regression for actual consumer installation"))
    end

    # Plot out the installation predicted by the economically rational model
    if (length(model_descriptions) > 0)
        for i in 1:length(model_descriptions)
            description = model_descriptions[i]
            model_prediction = model_predictions[i]
            ax1.plot(consumption, model_prediction, label = string("Optimal PV System for ", tariff_name, " Consumer: ", description))
        end
    else
        ax1.plot(consumption, model_predictions, c = "r", label = string("Optimal PV System for ", tariff_name, " Consumer"))
    end
    if x_max !== nothing
        plt.xlim(0,x_max)
    end
    if y_max !== nothing
        plt.ylim(0,y_max)
    end

    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    ylabel("PV System Capacity [kW]")
    xlabel("Consumer Monthly Energy use [kWh]")
    grid("on");
    title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_name))
end

# Function plotting all tariffs of a single tariff category, with a model prediction
function plot_single_tariff_category_per_company_with_model_prediction(company_data, tariff_category, company_name, consumption, model_predictions; model_descriptions=[], x_max=nothing, y_max=nothing, linreg=false)
    company_data_split_by_tariff_category = Array{DataFrame}(undef,3)

    # Loop through all tariffs that we care about
    for (tariff_num, tariff_category) in tariff_category_mappings
        next_value = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        category_index = findfirst(t -> t == tariff_category, tariff_categories)
        if !isassigned(company_data_split_by_tariff_category, category_index)
            company_data_split_by_tariff_category[category_index] = next_value
        else
            append!(company_data_split_by_tariff_category[category_index], next_value)
        end
    end

    category_index = findfirst(t -> t == tariff_category, tariff_categories)

    plot_base_GD_vs_economically_rational(create_consumption_and_installation_arrays(
                company_data_split_by_tariff_category[category_index]), tariff_categories[category_index], company_name, consumption, model_predictions; model_descriptions=model_descriptions, x_max=x_max, y_max=y_max, linreg=linreg)
end

# Function plotting each tariff of a single category individually
function plot_all_tariffs_per_company(company_data, company_name)
    company_data_split_by_tariff_category = Array{DataFrame}(undef,3)

    # Loop through all tariffs that we care about
    for (tariff_num, tariff_category) in tariff_category_mappings
        next_value = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        category_index = findfirst(t -> t == tariff_category, tariff_categories)
        if !isassigned(company_data_split_by_tariff_category, category_index)
            company_data_split_by_tariff_category[category_index] = next_value
        else
            append!(company_data_split_by_tariff_category[category_index], next_value)
        end
    end

    # Loop through and plot each tariff category
    for i in 1:3
        plot_consumption_and_installation(create_consumption_and_installation_arrays(
                company_data_split_by_tariff_category[i]), tariff_categories[i], company_name)
    end
end

# Function plotting a full tariff category with two regressions, divided at a pivot point
function plot_tariff_category_with_two_regressions(company_data, tariff_category, company_name, consumption, model_predictions, inflection_point=1500)

    fig = plt.figure()
    ax1 = fig.add_subplot(111)

    # Collect the tariff-relevant data
    tariff_data = Array{DataFrame}(undef,1)
    for (tariff_num, tariff_cat) in tariff_category_mappings
        if tariff_cat != tariff_category
            continue
        end
        tariff_df = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        if !isassigned(tariff_data, 1)
            tariff_data[1] = tariff_df
        else
            append!(tariff_data[1], tariff_df)
        end
    end

    # Get least-squares regression of the data
    t_consumption, t_installation = create_consumption_and_installation_arrays(tariff_data[1])
    combined = DataFrame()
    combined[:CONSUMPTION] = t_consumption
    combined[:INSTALLATION] = t_installation

    function map_to_float(str)
        try
            convert(Float64, str)
        catch
            return(NA)
        end
    end

    combined[:CONSUMPTION] = map(map_to_float, combined[:CONSUMPTION])
    combined[:INSTALLATION] = map(map_to_float, combined[:INSTALLATION])

    # Specific filter for data cleaning
    cons_upper_bound = company_name=="CNFL" ? 80000 : 400000
    combined = filter(row -> (row.CONSUMPTION < cons_upper_bound), combined)

    combined_1 = filter(row -> (row.CONSUMPTION <= inflection_point), combined)
    data_b_1, data_m_1 = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), combined_1))
    t_consumption_1 = combined_1[:CONSUMPTION]
    t_installation_1 = combined_1[:INSTALLATION]
    ax1.scatter(t_consumption_1, t_installation_1, c="g", marker=".", alpha=0.2, label = string("Actual PV System Installation for a ", tariff_category, " Consumer using less than ", inflection_point, " kWh per month"))

    combined_2 = filter(row -> (row.CONSUMPTION > inflection_point), combined)

    data_b_2, data_m_2 = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), combined_2))
    t_consumption_2 = combined_2[:CONSUMPTION]
    t_installation_2 = combined_2[:INSTALLATION]
    ax1.scatter(t_consumption_2, t_installation_2, c="m", marker=".", alpha=0.2, label = string("Actual PV System Installation for a ", tariff_category, " Consumer using more than ", inflection_point, " kWh per month"))

    x_vals = ax1.get_xlim()
    x_vals_1 = x_vals[1]:(inflection_point-x_vals[1])/100:inflection_point
    x_vals_2 = inflection_point:(x_vals[2]-inflection_point)/100:x_vals[2]
    y_vals_1 = data_m_1 * x_vals_1 .+ data_b_1
    y_vals_2 = data_m_2 * x_vals_2 .+ data_b_2
    ax1.plot(x_vals_1, y_vals_1, "--", c="g", label = string("Least-squares regression for actual installation of consumers using less than ", inflection_point, " kWh per month"))
    ax1.plot(x_vals_2, y_vals_2, "--", c="m", label = string("Least-squares regression for actual installation of consumers using more than ", inflection_point, " kWh per month"))

    # Run a regression on the model data
    ax1.plot(consumption, model_predictions, c = "r", label = string("Optimal PV System for ", tariff_category, " Consumer"))
    title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_category))
    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
#     colors = ["g", "c", "m", "y", "k"]
end

# Function plotting a full tariff category with multiple regressions
function plot_tariff_category_with_many_regressions(company_data, tariff_category, company_name, consumption, model_predictions, inflection_points=[])
    colors = ["g", "c", "m", "y", "k"]

    # Collect the tariff-relevant data
    tariff_data = Array{DataFrame}(undef,1)
    for (tariff_num, tariff_cat) in tariff_category_mappings
        if tariff_cat != tariff_category
            continue
        end
        tariff_df = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        if !isassigned(tariff_data, 1)
            tariff_data[1] = tariff_df
        else
            append!(tariff_data[1], tariff_df)
        end
    end

    # Get least-squares regression of the data
    t_consumption, t_installation = create_consumption_and_installation_arrays(tariff_data[1])
    combined = DataFrame()
    combined[:CONSUMPTION] = t_consumption
    combined[:INSTALLATION] = t_installation

    function map_to_float(str)
        try
            convert(Float64, str)
        catch
            return(NA)
        end
    end

    combined[:CONSUMPTION] = map(map_to_float, combined[:CONSUMPTION])
    combined[:INSTALLATION] = map(map_to_float, combined[:INSTALLATION])

    # Remove values where consumption is less than 30kW a month
    # This is a minimum requirement in CR, numbers lower than 30 indicate they're paying for but not using electricity
    combined = filter(row -> (row.CONSUMPTION >= 30), combined)
    early_consumption = company_name == "CNFL" ? 300 : 400
    max_early_install = company_name == "CNFL" ? 5 : 10
    combined = filter(row -> !(row.CONSUMPTION < early_consumption && row.INSTALLATION > max_early_install), combined)
    late_consumption = company_name == "CNFL" ? 8000 : 20000
    combined = filter(row -> (row.CONSUMPTION < late_consumption), combined)

    if length(inflection_points) == 0 || inflection_points[1] != 0
        pushfirst!(inflection_points, 0)
    end

    for (index, point) in enumerate(inflection_points)
        fig = plt.figure()
        ax1 = fig.add_subplot(111)
        color = popfirst!(colors)
        next_point = index == length(inflection_points) ? Inf : inflection_points[index + 1]
        combined_new = filter(row -> (point <= row.CONSUMPTION < next_point), combined)

        if size(combined_new)[1] == 0
            continue
        end
        data_b, data_m = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), combined_new))
        t_consumption_new = combined_new[:CONSUMPTION]
        t_installation_new = combined_new[:INSTALLATION]
        ax1.scatter(t_consumption_new, t_installation_new, c=color, marker=".", alpha=0.1, label = string("Actual PV System Installation for a ", tariff_category, " Consumer using between ", point, " and ", next_point, " kWh per month"))

        next_point = next_point == Inf ? ax1.get_xlim()[2] : next_point

        x_vals = point:(next_point-point)/100:next_point
        y_vals = data_m * x_vals .+ data_b
        ax1.plot(x_vals, y_vals, "--", c=color, label = string("Least-squares regression for actual installation of consumers using between ", point, " and ", next_point, " kWh per month"))
        # Run a regression on the model data
        ax1.plot(consumption, model_predictions, c = "r", label = string("Optimal PV System for ", tariff_category, " Consumer"))

        title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_category))
        legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    end
end

# Function to plot with a separate linear regression for each tariff within a tariff category.
function plot_with_subtariff_wise_regression(company_data, tariff_category, company_name, consumption, model_predictions)
    colors = ["g", "c", "m", "y", "k"]
    # Loop through all tariffs that we care about
    for (tariff_num, tariff_cat) in tariff_category_mappings
        if tariff_cat != tariff_category
            continue
        end
        fig = plt.figure()
        subtariff_values = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        # Get least-squares regression of the data
        t_consumption, t_installation = create_consumption_and_installation_arrays(subtariff_values)
        combined = DataFrame()
        combined[:CONSUMPTION] = t_consumption
        combined[:INSTALLATION] = t_installation

        function map_to_float(str)
            try
                convert(Float64, str)
            catch
                return(NA)
            end
        end

        combined[:CONSUMPTION] = map(map_to_float, combined[:CONSUMPTION])
        combined[:INSTALLATION] = map(map_to_float, combined[:INSTALLATION])

        # Specific filters for data cleaning
        if (tariff_num == 4 || tariff_num == 7)
            combined = filter(row -> (row.CONSUMPTION < 50000), combined)
        end

        color = popfirst!(colors)
        new_axis = fig.add_subplot(111)
        new_axis.scatter(combined[:CONSUMPTION], combined[:INSTALLATION], c=color, marker=".", alpha=0.2, label = string("Actual PV System Installation for a subtariff ", tariff_num, " Consumer"))
        xlims = (min(combined[:CONSUMPTION]...), max(combined[:CONSUMPTION]...))

        linear_model = lm(@formula(INSTALLATION ~ CONSUMPTION), combined)
        data_b, data_m = coef(linear_model)
        x_vals = xlims[1]:(xlims[2]-xlims[1])/100:xlims[2]
        y_vals = data_m * x_vals .+ data_b
        new_axis.plot(x_vals, y_vals, "--", c=color, label = string("Least-squares regression for actual installation of subtariff ", tariff_num, " consumers"))

        # Also includes standard error (standard deviation) bars on the line of best fit through the data
        stderr_b, stderr_m = stderror(linear_model)
        y_vals_2 = ((data_m + stderr_m) * x_vals .+ (data_b + stderr_b))
        new_axis.plot(x_vals, y_vals_2, "--", c="gray", label = string("Error on Least-squares regression"))
        y_vals_3 = ((data_m - stderr_m) * x_vals .+ (data_b - stderr_b))
        new_axis.plot(x_vals, y_vals_3, "--", c="gray", label = string("Error on Least-squares regression"))

        new_axis.plot(consumption, model_predictions, c = "r", label = string("Optimal PV System for ", tariff_category, " Consumer"))
        title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_to_name[tariff_num]))
        legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    end

end

# Function to plot a whole tariff category, each tariff being a different color
function plot_segmented_tariff_category_with_regression(company_data, tariff_category, company_name, consumption, model_predictions, regression_limit=0, individual_regressions=false)
    colors = ["g", "c", "m", "y", "k"]

    individual_tariffs = Dict{Integer,DataFrame}()
    # Loop through all tariffs that we care about
    for (tariff_num, tariff_cat) in tariff_category_mappings
        # Only look at tariffs within our specified tariff category
        if tariff_cat != tariff_category
            continue
        end
        tariff_df = filter(row -> (!ismissing(row.CODIGO_TARIFA) && row.CODIGO_TARIFA == string(tariff_num)), company_data)
        individual_tariffs[tariff_num] = tariff_df
    end

    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    title(string("PV System Capacity for ", company_name, " Consumer with Tariff ", tariff_category))

    all_tariff_data = Array{DataFrame}(undef,1)

    individual_consumption_and_installations = []

    for (tariff_num, tariff_df) in individual_tariffs
        if !isassigned(all_tariff_data, 1)
            all_tariff_data[1] = tariff_df
        else
            append!(all_tariff_data[1], tariff_df)
        end
        t_consumption, t_installation = create_consumption_and_installation_arrays(tariff_df)
        plot_color = pop!(colors)
        # Plot out the actual PV system sizes that people have based on their energy bills
        ax1.scatter(t_consumption, t_installation, c=plot_color, marker=".", alpha=0.5, label = string("Actual PV System Installation for a ", tariff_category, " Number ", tariff_num, " Consumer"))
        if individual_regressions == true
            push!(individual_consumption_and_installations, (t_consumption, t_installation, tariff_num, plot_color))
        end
    end

    # Get least-squares regression of the data
    all_con, all_ins = create_consumption_and_installation_arrays(all_tariff_data[1])
    combined = DataFrame()
    combined[:CONSUMPTION] = all_con
    combined[:INSTALLATION] = all_ins

    function map_to_float(str)
        try
            convert(Float64, str)
        catch
            return(NA)
        end
    end

    combined[:CONSUMPTION] = map(map_to_float, combined[:CONSUMPTION])
    combined[:INSTALLATION] = map(map_to_float, combined[:INSTALLATION])

    # Run regression on the entire dataset, and plot.
    data_b, data_m = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), combined))
    x_vals = ax1.get_xlim()
    x_vals = x_vals[1]:x_vals[2]-x_vals[1]/100:x_vals[2]
    y_vals = data_m * x_vals .+ data_b
    ax1.plot(x_vals, y_vals, "--", c="b", label = string("Least-squares regression for actual installation"))

    # If indicated, also run a regression on a subset of the data
    if regression_limit != 0
        limited_combined = filter(row -> (row.CONSUMPTION <= regression_limit), combined)
        limited_data_b, limited_data_m = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), limited_combined))
        y_vals = limited_data_m * x_vals .+ limited_data_b
        ax1.plot(x_vals, y_vals, "--", c="g", label = string("Least-squares regression for actual installation of consumers using less than ", regression_limit, " kWh per month"))
    end

    # Plot out the installation predicted by the economically rational model
    ax1.plot(consumption, model_predictions, c = "r", label = string("Optimal PV System for ", tariff_category, " Consumer"))

    # Run a regression on the model data
    model_df = DataFrame()
    model_df[:CONSUMPTION] = consumption
    model_df[:PREDICTIONS] = model_predictions

    model_df[:CONSUMPTION] = map(map_to_float, model_df[:CONSUMPTION])
    model_df[:PREDICTIONS] = map(map_to_float, model_df[:PREDICTIONS])

    model_data_b, model_data_m = coef(lm(@formula(PREDICTIONS ~ CONSUMPTION), model_df))
    y_vals = model_data_m * x_vals .+ model_data_b
    ax1.plot(x_vals, y_vals, "--", c="c", label = string("Least-squares regression for model prediction"))
    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
    ylabel("PV System Capacity [kW]")
    xlabel("Consumer Monthly Energy use [kWh]")
    grid("on");

    # Optionally include per-tariff regressions
    if individual_regressions == true
        num_figures = length(individual_consumption_and_installations)
        fig = plt.figure(figsize=(num_figures, 5*num_figures))
        colors = ["g", "c", "m", "y", "k"]
        for i in 1:length(individual_consumption_and_installations)
            new_axis = fig.add_subplot(num_figures, 1, i)
            for j in 1:length(individual_consumption_and_installations)
                t_consumption, t_installation, tariff_num, plot_color = individual_consumption_and_installations[j]
                new_axis.scatter(t_consumption, t_installation, c=plot_color, marker=".", alpha=0.5)
            end
            t_consumption, t_installation, tariff_num, plot_color = individual_consumption_and_installations[i]
            combined = DataFrame()
            combined[:CONSUMPTION] = t_consumption
            combined[:INSTALLATION] = t_installation
            # Run regression on the specific sub-tariff, and plot.
            data_b, data_m = coef(lm(@formula(INSTALLATION ~ CONSUMPTION), combined))
            x_vals = new_axis.get_xlim()
            x_vals = x_vals[1]:x_vals[2]-x_vals[1]/100:x_vals[2]
            y_vals = data_m * x_vals .+ data_b
            new_axis.plot(x_vals, y_vals, "--", c=plot_color, label = string("Least-squares regression for actual installation for tariff number ", tariff_num))
            legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);
            ylabel("PV System Capacity [kW]")
            xlabel("Consumer Monthly Energy use [kWh]")
            grid("on");
        end
    end

end

ext: jl

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
module ExternalCmds

    include("utils.jl")
    include("tee_file.jl")
    include("run.jl")

    export run_cmd, tee_run, run_bash

end


text:
module Extremes

using Distributions, DataFrames, Dates
using Optim
using LinearAlgebra, Statistics
using Mamba, ProgressMeter
using Gadfly

import CSV
import ForwardDiff, SpecialFunctions

import Distributions.quantile
import Statistics.var
import Base.length, Base.maximum, Base.sum

include("utils.jl")
include("structures.jl")
include("parameterestimation.jl")
include("data.jl")
include("validationplots.jl")

export

    # Variable type
    Variable, VariableStd, DataItem,

    # Cluster type
    Cluster,

    # Generic types
    EVA,
    fittedEVA,

    # Extreme value analysis type
    BlockMaxima,
    ThresholdExceedance,

    # Fitted extreme value analysis model
    pwmEVA,
    MaximumLikelihoodEVA,
    BayesianEVA,

    # Data related functions
    getcluster,

    # Fitting functions
    gevfit,
    gevfitbayes,
    gevfitpwm,
    gpfit,
    gpfitbayes,
    gpfitpwm,
    cint,

    # Other functions
    parametervar,

    # Return level
    ReturnLevel,
    returnlevel,

    # Diagnostic plots
    probplot_std_data,
    probplot_std,
    qqplot_std_data,
    qqplot_std,
    diagnosticplots_std,

    probplot_data,
    probplot,
    qqplot_data,
    qqplot,
    returnlevelplot_data,
    returnlevelplot,
    histplot_data,
    histplot,
    diagnosticplots,
    mrlplot,
    mrlplot_data

end # module


text:
import YAML
config = YAML.load(open("config.yaml"))

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using MultivariateStats
using HDF5

dv = DataVis(config["data_file"], true)[5]

# Conjugation is necessary for the SMA and ALMA
visibilities.conj!(dv)

qq = get_qq(dv)

# Find the indices that sort these in increasing order
ind = sortperm(qq)
qq = qq[ind]
uu = dv.uu[ind]
vv = dv.vv[ind]

# Points at which to create images
parameters = readdlm("parameters.dat")
npoints = size(parameters, 1)

fid = h5open("grid_VV.hdf5", "r")
# Load the sampled visibilities from HDF5
reals = read(fid, "real")
imags = read(fid, "imag")
VV = reals + imags .* im # Complex visibility

close(fid)

nvis = size(reals, 1)

println("Dimensions of input matrix ", size(reals))

# suppose Xtr and Xte are training and testing data matrix,
# with each observation in a column

# Need to do the reals and imaginaries separately, since MultivariateStats does not work w/ Complex128

mean_estimate = mean(reals, 2)
std_estimate = std(reals, 2)

norm_reals = (reals .- mean_estimate) ./ std_estimate
# norm_reals = reals .- mean_estimate

# train a PCA model
M = fit(PCA, norm_reals; mean=0, pratio=0.999)

dump(M)

eigen = projection(M)

println("projection", size(eigen))

ncomp = outdim(M)

# a (npoints, ncomp) matrix
weights = Array(Float64, (npoints, ncomp))

# Calculate the weights that lead to reconstruction of the visibilities
for i=1:npoints
  for j=1:ncomp
    weights[i,j] = sum(reals[:,i] .* eigen[:,j])
  end
end

# Now try reconstructing the spectra
reconstructed_reals = Array(Float64, (nvis, npoints))
for i=1:npoints
  reconstructed_reals[:,i] = std_estimate .* ((eigen * squeeze(weights[i,:], 1)) .+ mean_estimate)
  # reconstructed_reals[:,i] = (eigen * squeeze(weights[i,:], 1)) .+ mean_estimate
end

ms=0.1

# Compare the reconstructed spectra vs the reals.
import PyPlot.plt
using LaTeXStrings

fig, ax = plt[:subplots](nrows=3, figsize=(10, 6), sharex=true)

for i=1:3
  ax[1][:plot](qq, reals[:,i], ".", ms=ms)
  ax[2][:plot](qq, reconstructed_reals[:,i], ".", ms=ms)
  ax[3][:plot](qq, reals[:,i] - reconstructed_reals[:,i], ".", ms=ms)
end

ax[1][:set_xlabel](L"$qq\quad k\lambda$")
ax[1][:grid]()
ax[2][:grid]()
ax[3][:grid]()

plt[:savefig]("visibilities_reconstructed.png", dpi=600)

fig, ax = plt[:subplots](nrows=3, figsize=(10, 6), sharex=true)

ax[1][:plot](qq, mean_estimate, ".", ms=ms)
ax[2][:plot](qq, std_estimate, ".", ms=ms)

for i=1:size(eigen, 2)
  ax[3][:plot](qq, eigen[:,i], ".", ms=ms)
end

ax[3][:set_xlabel](L"$qq\quad k\lambda$")

ax[1][:grid]()
ax[2][:grid]()
ax[3][:grid]()

ax[1][:set_ylabel]("Mean")
ax[2][:set_ylabel]("Std")
ax[3][:set_ylabel]("Real")

plt[:savefig]("visibilities_PCA_real.png", dpi=600)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
using MCMCBenchmarkSuite, Test, Random

@testset "Regression Tests " begin
    Random.seed!(22385)
    path = pathof(MCMCBenchmarkSuite)
    include(joinpath(path,
      "../Linear_Regression/Linear_Regression_Models.jl"))
    Nd = 50
    β0 = 1
    β = fill(.5,1)
    σ = 1
    Nreps = 1
    ProjDir = @__DIR__
    cd(ProjDir)
    samplers = (
      CmdStanNUTS(CmdStanConfig,ProjDir),
      AHMCNUTS(AHMCregression,AHMCconfig),
      DHMCNUTS(sampleDHMC))
    options = (Nsamples=2000, Nadapt=1000, delta=.8, Nd=Nd)
    results = benchmark(samplers, simulateRegression, Nreps;options...)
    @test isa(results, DataFrame)
    # @test results[!,Symbol("B[1]_mean")][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ β[1] atol = .1
    # @test results[!,Symbol("B[1]_mean")][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ β[1] atol = .1
    # @test results[!,Symbol("B[1]_mean")][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ β[1] atol = .1
    # @test results[!,:B0_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ β0[1] atol = .1
    # @test results[!,:B0_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ β0[1] atol = .1
    # @test results[!,:B0_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ β0[1] atol = .1
    # @test results[!,:sigma_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ σ atol = .1
    # @test results[!,:sigma_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ σ atol = .1
    # @test results[!,:sigma_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ σ atol = .1
    # @test results[!,Symbol("B[1]_sampler_rhat")][1] ≈ 1 atol = .03
    # @test results[!,:B0_sampler_rhat][1] ≈ 1 atol = .03
    # @test results[!,:sigma_sampler_rhat][1] ≈ 1 atol = .03
    isdir("tmp") && rm("tmp", recursive=true)
end


text:
using VannModels
using IcanProj
using JLD2
using PyPlot
using DataFrames

function snow_sim(time, prec, tair)

    ntimes = size(prec, 1)
    nregion = size(prec, 2)

    frac_lus = DataFrame(dummy = fill(1/nregion, nregion))

    mdata = TinSnow(24.0, time[1], frac_lus)

    swe = fill(0.0, size(prec))

    for i in 1:ntimes

        mdata.p_in .= prec[i,:]
        mdata.tair .= tair[i,:]

        run_timestep(mdata)

        swe[i,:] .= mdata.swe[1,:]

    end

    return swe

end

opt = get_options()

watershed = opt["stat_sel"][3]

# Load results

file = joinpath(opt["eval_folder"], watershed, "res_1km.jld2")

@load file res_vic

# Get forcing variables

time, prec = get_variable(res_vic, "prec", average = false)

time, tair = get_variable(res_vic, "air_temp", average = false)

# Run snow model

swe = snow_sim(time, prec, tair)


text:
@doc raw"""
    Models.Eggbox()

Eggbox/Egg carton likelihood function

```math
z(x, y) = \left[a + \cos\frac{x}{b} \cdot \cos\frac{x}{b} \right]^5
```

# Examples
```jldoctest
julia> model, lnZ = Models.Eggbox();

julia> lnZ
235.88
```
"""
function Eggbox()
    tmax = 5π

    # uniform prior from 0, 1
    prior(X) = X
    function loglike(X)
        a = cos(tmax * (2 * first(X) - 1) / 2)
        b = cos(tmax * (2 * last(X) - 1) / 2)
        return (2 + a * b)^5
    end

    lnZ = 235.88 # where do we get this from??
    return NestedModel(loglike, prior), lnZ
end


text:
module test_julia_macro_invocation

using Test

macro m()
    (__module__, __source__)
end

(m, node) = @m
@test m === test_julia_macro_invocation
@test node === LineNumberNode(9, Symbol(@__FILE__))

end # module test_julia_macro_invocation

module test_julia_macro_emoji

using Test

macro 🤔(f)
    # println(f)
    f
end

function 🍞(x)
    x, x
end

function 🚗(x)
    x
end

f = @🤔 🍞 ∘ 🚗
@test f('🌀') == ('🌀', '🌀')

end # module test_julia_macro_emoji

module test_julia_macro_var_str

using Test

struct A
    var"local"
end

a = A(1)
@test a.local == a.var"local" == getproperty(a, :local)

end # module test_julia_macro_var_str


text:
#Maximum Compatibility Pairing
using JuMP,Gurobi

#define the model
m =Model(Gurobi.Optimizer)

#Let n represents the number of cities to be paired
n=4

#declaring variables
@variable(m, x[1:n,1:n]>=0)

#define the compatibility coefficients
c =   [	80 65 83 77;
	54 87 61 66 ;
	92 45 53 59 ;
	70 61 81 76 ]

#define the objective function
@objective(m, Max, sum(c[i,j]* x[i,j] for i in 1:n for j in 1: n)) #objective is to maximize total compatibility pairing

#define constraints
for j = 1: n
  @constraint(m ,sum(x[i,j] for i in 1: n) ==1)
end

for i = 1: n
  @constraint(m ,sum(x[i,j] for j in 1: n) ==1 )
end

print(m)

#Printing the optimal solutions obtained
println("Optimal Solutions:")
status = JuMP.optimize!(m)
println(" Objective value: ",JuMP.objective_value(m))

println("Decision variable values:")
println("x = ", JuMP.value.(x))

max_stars_repo_name: Oyekanmi/Optimization-Julia

max_issues_repo_name: Oyekanmi/Optimization-Julia

max_forks_repo_name: Oyekanmi/Optimization-Julia


text:
#####################
# Downward Flow circuits
#####################

"Representation of downward flows with `in_progress` Bool for effective book keeping"

abstract type DownFlowΔNode{O,F} <: DecoratorΔNode{UpFlowΔNode{O,F}} end
abstract type DownFlowInnerNode{O,F} <: DownFlowΔNode{O,F} end

abstract type DownFlow⋀{O,F} <: DownFlowInnerNode{O,F} end
abstract type DownFlow⋁{O,F} <: DownFlowInnerNode{O,F} end

# use this version of DownFlow⋀ when ⋀ gates are unique to their parent ⋁ to save space (compact=true)
mutable struct DownFlow{F}
    downflow::F
    in_progress::Bool
end

struct DownFlowLeaf{O,F} <: DownFlowΔNode{O,F}
    origin::UpFlowLeafNode{O,F}
    downflow::DownFlow{F}
end

struct DownFlow⋀Compact{O,F} <: DownFlow⋀{O,F}
    origin::UpFlow⋀{O,F}
    children::Vector{<:DownFlowΔNode{<:O,F}}
    cached_downflow_sinks::Vector{DownFlow{F}}
end

# use this version of DownFlow⋀ when ⋀ gates can have multiple parents to save computation (compact=false)
struct DownFlow⋀Cached{O,F} <: DownFlow⋀{O,F}
    origin::UpFlow⋀{<:O,<:F}
    children::Vector{<:DownFlowΔNode{<:O,F}}
    downflow::DownFlow{F}
end

# this version of DownFlow⋁ is efficient for arity-1 gates
struct DownFlow⋁Compact{O,F} <: DownFlow⋁{O,F}
    origin::UpFlow⋁{<:O,<:F}
    child::DownFlowΔNode{<:O,F}
    cached_downflow_sinks::Vector{DownFlow{F}}
end

struct DownFlow⋁Cached{O,F} <: DownFlow⋁{O,F}
    origin::UpFlow⋁{<:O,<:F}
    children::Vector{<:DownFlowΔNode{<:O,F}}
    downflow::DownFlow{F}
end

const DownFlowΔ{O,F} = AbstractVector{<:DownFlowΔNode{O,F}}

const FlowΔNode{O,F} = DownFlowΔNode{O,F}
const FlowΔ{O,F} = AbstractVector{<:FlowΔNode{O,F}}

#####################
# traits
#####################

@inline GateType(::Type{<:DownFlowLeaf{O}}) where O = GateType(O)
@inline GateType(::Type{<:DownFlow⋀}) = ⋀Gate()
@inline GateType(::Type{<:DownFlow⋁}) = ⋁Gate()

const HasDownFlow = Union{DownFlow⋁Cached, DownFlow⋀Cached, DownFlowLeaf}

#####################
# constructors and conversions
#####################

"""Construct a downward flow circuit from a given upward flow circuit"""
function DownFlowΔ(circuit::UpFlowΔ{O,F}, opts = flow_opts★)::DownFlowΔ{O,F} where {O,F}

    m = flow_length(circuit)
    fmem  = () -> some_vector(eltype(F), m)

    cache = Dict{ΔNode, DownFlowΔNode}()
    sizehint!(cache, length(circuit)*4÷3)

    flow_node(::LeafGate, n::ΔNode) = DownFlowLeaf{O,F}(n, DownFlow(fmem(),false))

    flow_node(::⋀Gate, n::ΔNode) = begin
        n_children = map(c -> cache[c], children(n))
        @assert length(n_children)>=2
        sinks = downflow_sinks(n_children)
        if opts.compact⋀ && length(sinks) <= opts.max_factors
            DownFlow⋀Compact{O,F}(n, n_children, sinks)
        else
            DownFlow⋀Cached{O,F}(n, n_children, DownFlow(fmem(),false))
        end
    end

    flow_node(::⋁Gate, n::ΔNode) = begin
        n_children = map(c -> cache[c], children(n))
        @assert length(n_children)>=1
        sinks = downflow_sinks(n_children)
        if opts.compact⋁ && length(n_children)==1 && length(sinks) <= opts.max_factors
            DownFlow⋁Compact{O,F}(n, n_children[1], sinks)
        else
            DownFlow⋁Cached{O,F}(n, n_children, DownFlow(fmem(),false))
        end
    end

    map(circuit) do node
        fnode = flow_node(GateType(node), node)
        cache[node] = fnode
        fnode
    end
end

"""Construct a up and down flow circuit from a given other circuit"""
function FlowΔ(circuit::Δ, m::Int, ::Type{El}, opts = flow_opts★) where El
    up_flow_circuit = UpFlowΔ(circuit, m, El, opts)
    DownFlowΔ(up_flow_circuit, opts)
end

#####################
# methods
#####################

@inline literal(n::DownFlowLeaf)::Lit  = literal(n.origin)
@inline constant(n::DownFlowLeaf)::Bool = constant(n.origin)
@inline children(n::DownFlow⋁Compact) = [n.child]
@inline children(n::Union{DownFlow⋀,DownFlow⋁Cached}) = n.children

@inline downflow(n::HasDownFlow) = n.downflow.downflow

@inline downflow_sinks(n::DownFlowInnerNode) = n.cached_downflow_sinks
@inline downflow_sinks(n::HasDownFlow) = [n.downflow]
# @inline downflow_sinks(n::DownFlowLeaf{O,F}) where {O,F} = Vector{DownFlow{F}}()
@inline function downflow_sinks(ns::Vector{<:DownFlowΔNode{O,F}})::Vector{DownFlow{F}} where {O,F}
    flatmap(c -> downflow_sinks(c)::Vector{DownFlow{F}}, ns, DownFlow{F}[])
end

pr(n::DownFlowΔNode) = pr(origin(n))
pr_factors(n::DownFlowΔNode) = pr_factors(origin(n))

flow_length(circuit::DownFlowΔ) = length(downflow_sinks(circuit[end])[1].downflow)
origin_flow_length(circuit::DownFlowΔ) = flow_length(origin(circuit[end]))

#####################

# change the size of the flow vectors throughout this circuit (e.g., change minibatch size)
function resize_flows(circuit::DownFlowΔ{F}, size::Int) where {F}
    # check if resize is needed
    if size != flow_length(circuit)
        for n in circuit
            resize_path_flow(n,size)
        end
    end
end

@inline resize_path_flow(n::DownFlowΔNode, size) = () #do nothing
@inline resize_path_flow(n::HasDownFlow, size) = resize!(n.downflow.downflow, size)

#####################

function pass_down(circuit::DownFlowΔ{O,F}) where {F,O}
    resize_flows(circuit, flow_length(origin(circuit)))
    for n in circuit
        reset_downflow_in_progress(n)
    end
    for downflow in downflow_sinks(circuit[end])
        # initialize root flows to 1
        downflow.downflow .= one(eltype(F))
    end
    for n in Iterators.reverse(circuit)
        pass_down_node(n)
    end
end

reset_downflow_in_progress(n::DownFlowΔNode) = () # do nothing
reset_downflow_in_progress(n::HasDownFlow) = (n.downflow.in_progress = false)

pass_down_node(n::DownFlowΔNode) = () # do nothing
pass_down_node(n::DownFlowLeaf) = () # too costly? @assert pr(n) ≈ downflow(n) "Postitive leaf flows should be equal to their probabilities: $(pr(n)) ≆ $(downflow(n))"

function pass_down_node(n::DownFlow⋀Cached)
    for c in n.children
        for sink in downflow_sinks(c)
            if !sink.in_progress
                assign(sink.downflow, downflow(n))
                sink.in_progress = true
            else
                accumulate_val(sink.downflow, downflow(n))
            end
        end
    end
end

function pass_down_node(n::DownFlow⋁Cached)
    if num_children(n) == 1
        for sink in downflow_sinks(n.children[1])
            if !sink.in_progress
                assign(sink.downflow, downflow(n))
                sink.in_progress = true
            else
                accumulate_val(sink.downflow, downflow(n))
            end
        end
    else
        for c in n.children
            pr_fs = pr_factors(c)
            for sink in downflow_sinks(c)
                if !sink.in_progress
                    assign_prod_normalized(sink.downflow, pr(n), downflow(n), pr_fs)
                    sink.in_progress = true
                else
                    accumulate_prod_normalized(sink.downflow, pr(n), downflow(n), pr_fs)
                end
            end
        end
    end
end

#####################

function pass_up_down(circuit::DownFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {F,O}
    pass_up(origin(circuit), data)
    pass_down(circuit)
end


text:
using GUIFrameworkPlayground
const G = GUIFrameworkPlayground

b = G.SDL2GUIBackend()
#b = G.OpenGLGUIBackend()

win = G.window(b,800,600)

for i=1:300

    x,y = G.mouse_position(win)

    G.clear!(win)
    G.rectangle(win,x,y,200,100)
    G.text(win,x,y,"rectangle")

    ev = G.events(win)
    for i=1:length(ev)

        if typeof(ev[i]) == KeyEvent
            G.text(win,10,20*i+10, string(ev[i]))
        end
    end

    G.update!(win)
    sleep(1/30)
end

G.destroy(win)


text:
module mcompshell0
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.1
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 2000*phun("mm")
    nx = ny = 8
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    npairs = 1
    thickness = 10*phun("mm");
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_0_$p", mater, thickness/npairs/2, 0))
        push!(plies, CM.Ply("ply_90_$p", mater, thickness/npairs/2, 90))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = Float64[-Inf Inf 0 0 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0.1*phun("MPa")*thickness, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    l1 = selectelem(fens, bfes, box = Float64[-Inf Inf ay ay 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, -0.1*phun("MPa")*thickness, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show  maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])
    @test maximum(dchi.values[:, 3]) ./phun("mm") ≈ 0.22004349767718365
    true
end
end
using .mcompshell0
mcompshell0.test()

module mcompshell1
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.1
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 2000*phun("mm")
    nx = ny = 8
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    npairs = 10
    thickness = 10*phun("mm");
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_0_$p", mater, thickness/npairs/2, 0))
        push!(plies, CM.Ply("ply_90_$p", mater, thickness/npairs/2, 90))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = Float64[-Inf Inf 0 0 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0.1*phun("MPa")*thickness, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    l1 = selectelem(fens, bfes, box = Float64[-Inf Inf ay ay 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, -0.1*phun("MPa")*thickness, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show  maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])
    @test maximum(dchi.values[:, 3]) ./phun("mm") ≈ 0.010276258430770539
    true
end
end
using .mcompshell1
mcompshell1.test()

module mcompshell2
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.3
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 2000*phun("mm")
    nx = ny = 8
    # ASFD/9310
    E1 = 37880*phun("MPa")
    E2 = 9407*phun("MPa")
    G12 = 3405*phun("MPa")
    G13 = G12
    nu12 = 0.299;
    nu23 = 0.422
    G23 = 3308*phun("MPa")
    ply_thickness = 1*phun("mm");
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[
    CM.Ply("p", mater, ply_thickness, +45),
    CM.Ply("p", mater, ply_thickness, -45),
    CM.Ply("p", mater, ply_thickness, 0),
    CM.Ply("p", mater, ply_thickness, 0),
    CM.Ply("p", mater, ply_thickness, -45),
    CM.Ply("p", mater, ply_thickness, +45),
    ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.3", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), length(plies)*ply_thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [1, 2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = Float64[0 0 -Inf Inf 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-200*phun("N/mm"), 0, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    l1 = selectelem(fens, bfes, box = Float64[ax ax -Inf Inf 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[200*phun("N/mm"), 0, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-3.3-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

     # The book: The target value is 1.667 mm at either of the loaded edges; so
     # the total displacement is approximately 2*1.667=3.334. The Abaqus model
     # on the other hand gives 1.658*2 = 3.316.
    @test maximum(dchi.values[:, 1]) ./phun("mm") ≈ 3.3161098653904544
    true
end
end
using .mcompshell2
mcompshell2.test()

module mcompshell3
# Example  from Section 6.17.1
# Lazarus Teneketzis Tenek, John Argyris, Finite Element Analysis for Composite Structures, 1998
# Clamped isotropic plate under uniform transverse pressure
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test_composite()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule
    E = 200*phun("GPa");
    nu = 0.3;
    L = 1.0*phun("m");
    t_radius_ratio = 0.001
    thickness = L * t_radius_ratio;
    q = 1.0*phun("kilo*Pa");
    D = E*thickness^3/12/(1-nu^2)
    wcref = 0.00126*q*L^4/D
    nx = ny = 4
    ply_thickness = thickness;
    tolerance = L/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E, nu)
    plies = CM.Ply[
    CM.Ply("p", mater, ply_thickness, 0),
    ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("Tenek, Argyris 6.17.1", plies, mcsys)

    fens, fes = T3block(L,L,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), CM.thickness(layup)), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped Condition
    l1 = connectednodes(meshboundary(fes))
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i])
        # @show minimum(dchi.values[:, i])
    end
    # vtkwrite("plate-6.17.1-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test isapprox(-minimum(dchi.values[:, 3]), wcref, rtol=0.02)
    true
end,
function test_homogeneous()
    formul = FEMMShellT3FFModule
    E = 200*phun("GPa");
    nu = 0.3;
    L = 1.0*phun("m");
    t_radius_ratio = 0.001
    thickness = L * t_radius_ratio;
    q = 1.0*phun("kilo*Pa");
    D = E*thickness^3/12/(1-nu^2)
    wcref = 0.00126*q*L^4/D
    nx = ny = 4
    ply_thickness = thickness;
    tolerance = L/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E, nu)

    fens, fes = T3block(L,L,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped Condition
    l1 = connectednodes(meshboundary(fes))
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i])
        # @show minimum(dchi.values[:, i])
    end
    # vtkwrite("plate-6.17.1-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test isapprox(-minimum(dchi.values[:, 3]), wcref, rtol=0.02)
    true
end
end
using .mcompshell3
mcompshell3.test_composite()
mcompshell3.test_homogeneous()

module mcompshell4
# Example  from Section 6.17.1
# Lazarus Teneketzis Tenek, John Argyris, Finite Element Analysis for Composite Structures, 1998
# Clamped isotropic plate under uniform transverse pressure
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test_composite()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule
    E = 200*phun("GPa");
    nu = 0.3;
    L = 1.0*phun("m");
    t_radius_ratio = 0.001
    thickness = L * t_radius_ratio;
    q = 1.0*phun("kilo*Pa");
    D = E*thickness^3/12/(1-nu^2)
    wcref = 0.00126*q*L^4/D
    nx = ny = 4
    ply_thickness = thickness/3;
    tolerance = L/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E, nu)
    plies = CM.Ply[
    CM.Ply("p", mater, ply_thickness, 0),
    CM.Ply("p", mater, ply_thickness, 45),
    CM.Ply("p", mater, ply_thickness, 90),
    ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("Tenek, Argyris 6.17.1", plies, mcsys)

    fens, fes = T3block(L,L,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), CM.thickness(layup)), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped Condition
    l1 = connectednodes(meshboundary(fes))
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i])
        # @show minimum(dchi.values[:, i])
    end
    # vtkwrite("plate-6.17.1-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test isapprox(-minimum(dchi.values[:, 3]), wcref, rtol=0.02)
    true
end

end
using .mcompshell4
mcompshell4.test_composite()

module mcompshell5
# Example  from Section 6.17.2
# Lazarus Teneketzis Tenek, John Argyris, Finite Element Analysis for Composite Structures, 1998
# Simply supported sandwich plate.
# The deflection provided in the textbook of 31.45 mm is probably wrong.
# The deflection of 2.026 mm was obtained with Abaqus.
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule
    Ef = 6.8*phun("GPa");
    nuf = 0.3;
    Ec = 480*phun("MPa");
    nuc = 0.3
    L = 0.1384*phun("m");
    tf = 0.001*phun("m");
    tc = 0.015*phun("m");
    q = 1.0*phun("MPa");
    nx = ny = 140
    tolerance = L/nx/100
    CM = CompositeLayupModule

    plies = CM.Ply[
    CM.Ply("face_1", CM.lamina_material(Ef, nuf), tf, 0),
    CM.Ply("core", CM.lamina_material(Ec, nuc), tc, 0),
    CM.Ply("face_2", CM.lamina_material(Ef, nuf), tf, 0),
    ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("Tenek, Argyris 6.17.2", plies, mcsys)

    fens, fes = T3block(L,L,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), CM.thickness(layup)), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[L L 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
        # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
        # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Face load
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    # for i in 1:3
    #     @show maximum(dchi.values[:, i]) ./phun("mm")
    #     @show minimum(dchi.values[:, i]) ./phun("mm")
    # end
    # vtkwrite("plate-6.17.2-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test isapprox(-minimum(dchi.values[:, 3]), 2.026e-3, rtol=0.02)
    true
end

end
using .mcompshell5
mcompshell5.test()

module mcompshell6
# Similar definition of the plate as in:
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.1
# but the loading is transverse, pressure of magnitude 1.0
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 1000*phun("mm") # note the changed dimension
    nx = ny = 50
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    q = 1.0*phun("MPa")
    npairs = 1
    thickness = 10*phun("mm");
    # With these inputs, the Abaqus-verified solution is -1502.0324040447633

    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_0_$p", mater, thickness/npairs/2, 0))
        push!(plies, CM.Ply("ply_90_$p", mater, thickness/npairs/2, 90))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show  maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -1502.0324040447633
    true
end
end
using .mcompshell6
mcompshell6.test()

module mcompshell7
# Similar definition of the plate as in mcompshell6 but the layup is [-45/45]
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 1000*phun("mm") # note the changed dimension
    nx = ny = 50
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    q = 1.0*phun("MPa")
    npairs = 1
    thickness = 10*phun("mm");
    # With these inputs, the Abaqus-verified solution is -1176.8346659735575

    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    for p in 1:npairs
        push!(plies, CM.Ply("ply_-45_$p", mater, thickness/npairs/2, -45))
        push!(plies, CM.Ply("ply_+45_$p", mater, thickness/npairs/2, +45))
    end
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1, 2, 3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show  maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-$npairs-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -1176.8346659735575
    true
end
end
using .mcompshell7
mcompshell7.test()

module mcompshell8
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.1
# Variation which tests the extension-bending coupling: unsymmetric cross-ply
# layup should bend one way under compression, and another way under tension.
# Two plies, 90 degrees at the bottom. Applied tensile load in the zero
# direction (the reference direction of the composite layup) should bend the
# plate downwards (in the negative z direction). And vice versa.

using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(load_multiplier)
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 2000*phun("mm")
    nx = ny = 8
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    thickness = 10*phun("mm");
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    push!(plies, CM.Ply("ply_90", mater, thickness/2, 90))
    push!(plies, CM.Ply("ply_0", mater, thickness/2, 0))
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[ax ax 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = Float64[0 0 -Inf Inf 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0.1*phun("MPa")*load_multiplier*thickness, 0, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    l1 = selectelem(fens, bfes, box = Float64[ax ax -Inf Inf 0 0 ], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-0.1*phun("MPa")*load_multiplier*thickness, 0, 0, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    if load_multiplier < 0
        @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -0.22004349767718365
    else
        @test maximum(dchi.values[:, 3]) ./phun("mm") ≈ 0.22004349767718365
    end
    true
end
end
using .mcompshell8
mcompshell8.test(+1)
mcompshell8.test(-1)

module mcompshell9
# As mcompshell8, but the plate is oriented in the yz plane, with the loading in
# the z direction.

using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test(load_multiplier)
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    ax = ay = 2000*phun("mm")
    nx = ny = 8
    # ASFD/9310
    E1 = 133860*phun("MPa")
    E2 = 7706*phun("MPa")
    G12 = 4306*phun("MPa")
    G13 = G12
    nu12 = 0.301;
    nu23 = 0.396
    G23 = 2760*phun("MPa")
    thickness = 10*phun("mm");
    tolerance = ax/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[]
    push!(plies, CM.Ply("ply_90", mater, thickness/2, 90))
    push!(plies, CM.Ply("ply_0", mater, thickness/2, 0))
    mcsys = CSys(Float64[
        0 0 1;
        1 0 0;
        0 1 0])
    layup = CM.CompositeLayup("example_3.1", plies, mcsys)

    fens, fes = T3block(ax,ay,nx,ny);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        x, y, z = fens.xyz[i, :]
        fens.xyz[i, :] .= (z, x, y)
    end
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Pin one of the corners
    l1 = selectnode(fens; box = Float64[0 0 0 0 0 0], inflate = tolerance)
    for i in [1,2, 3,]
        setebc!(dchi, l1, true, i)
    end
    # Roller at the other
    l1 = selectnode(fens; box = Float64[0 0 0 0 ax ax], inflate = tolerance)
    for i in [2,]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    # Edge load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = Float64[0 0 -Inf Inf ax ax], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0.1*phun("MPa")*load_multiplier*thickness, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    l1 = selectelem(fens, bfes, box = Float64[0 0 -Inf Inf 0 0], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, -0.1*phun("MPa")*load_multiplier*thickness, 0, 0, 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    for i in 1:3
        # @show maximum(dchi.values[:, i]) ./phun("mm")
        # @show minimum(dchi.values[:, i]) ./phun("mm")
    end
    # vtkwrite("plate-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    if load_multiplier < 0
        @test minimum(dchi.values[:, 1]) ./phun("mm") ≈ -0.22004349767718365
    else
        @test maximum(dchi.values[:, 1]) ./phun("mm") ≈ 0.22004349767718365
    end
    true
end
end
using .mcompshell9
mcompshell9.test(+1)
mcompshell9.test(-1)

module mcompshell10
# A mixed-enhanced finite-element for the analysis of laminated composite plates
# FERDINANDO AURICCHIO, ELIO SACCO

# Laminated [0/90] cross-ply square simply supported plate. Uniform Distributed
# Loading. Results in Table IV.

# Note: The scaling of the horizontal displacement of the midpoint of the edge
# is wrong in the paper. The power of the aspect ratio should also be 4.
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    a = 250*phun("mm") # note the changed dimension
    nx = ny = 8
    # high modulus orthotropic graphite=epoxy composite material
    E1 = 100*phun("GPa")
    E2 = E1/25
    G12 = G13 = E2*0.5
    nu12 = 0.25;
    nu23 = 0.25
    G23 = E2*0.2
    q = 10.0*phun("kilo*Pa")
    thickness = a/10;
    # With these inputs, the Abaqus-verified solution is -1176.8346659735575

    tolerance = a/nx/100
    CM = CompositeLayupModule

    mater = CM.lamina_material(E1, E2, nu12, G12, G13, G23)
    plies = CM.Ply[
        CM.Ply("ply_0", mater, thickness/2, 0),
        CM.Ply("ply_90", mater, thickness/2, 90)
        ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("sacco", plies, mcsys)

    fens, fes = T3block(a,a,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    # The boundary conditions are a bit peculiar: refer to the paper
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf 0 0], inflate = tolerance)
    for i in [2, 4]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[a a -Inf Inf 0 0], inflate = tolerance)
    for i in [2, 4]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 0 0], inflate = tolerance)
    for i in [1, 5]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf a a 0 0], inflate = tolerance)
    for i in [1, 5]
        setebc!(dchi, l1, true, i)
    end
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [3]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])

    # @show maximum(dchi.values[:, 3]) * (100*E2) / (q*thickness*(a/thickness)^4)
    # @show maximum(dchi.values[:, 1]) * (100*E2) / (q*thickness*(a/thickness)^4)
    @test maximum(dchi.values[:, 3]) * (100*E2) / (q*thickness*(a/thickness)^4) ≈ 1.9529295928169146
    @test maximum(dchi.values[:, 1]) * (100*E2) / (q*thickness*(a/thickness)^4) ≈ 0.12375916912652352

    # vtkwrite("plate--uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

    # @test minimum(dchi.values[:, 3]) ./phun("mm") ≈ -1176.8346659735575
    true
end
end
using .mcompshell10
mcompshell10.test()

module mcompshell11
# From Barbero's Finite Element Analysis using Abaqus ... book Example 3.7
# Here we actually solve original benchmark:
# R0031(3): Three-layer sandwich shell under normal pressure loading
using LinearAlgebra: norm, Transpose, mul!, I
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CompositeLayupModule
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFCompModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using Test

function test()
    formul = FEMMShellT3FFCompModule
    CM = CompositeLayupModule

    a = 2*127*phun("mm")
    nx = ny = 18

    E1 = 68947.57*phun("MPa")
    E2 = 27579.03*phun("MPa")
    G12 = G13 = 12927.67*phun("MPa")
    nu12 = 0.3;
    G23 = 12927.67*phun("MPa")
    thickness_face = 0.7112*phun("mm")
    face = CM.lamina_material(E1, E2, nu12, G12, G13, G23)

    E1 = 0.068948*phun("MPa")
    E2 = 0.068948*phun("MPa")
    nu12 = 0.0;
    G12 = 0.068948*phun("MPa")
    G13 = 209.6*phun("MPa")
    G23 = 82.737*phun("MPa")
    thickness_core = 19.05*phun("mm")
    core = CM.lamina_material(E1, E2, nu12, G12, G13, G23)

    q = 0.689476*phun("MPa")

    tolerance = a/nx/100

    plies = CM.Ply[
        CM.Ply("ply_face_1", face, thickness_face, 0),
        CM.Ply("ply_core", core, thickness_core, 0),
        CM.Ply("ply_face_2", face, thickness_face, 0),
        ]
    mcsys = CM.cartesian_csys((1, 2, 3))
    layup = CM.CompositeLayup("Barbero 3.7", plies, mcsys)

    fens, fes = T3block(a,a,nx,ny);
    fens.xyz = xyz3(fens)
    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)

    femm = formul.make(IntegDomain(fes, TriRule(1), CM.thickness(layup)), layup)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Simple support
    l1 = connectednodes(meshboundary(fes))
    for i in [1,2,3,]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    formul.associategeometry!(femm, geom0)
    K = formul.stiffness(femm, geom0, u0, Rfield0, dchi);

    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])

    # R0031(3): Three-layer sandwich shell under normal pressure loading
    # lists 0.123'' = 3.1242 mm
    @test maximum(dchi.values[:, 3])./phun("mm") ≈ 3.2024601774900674

    true
end
end
using .mcompshell11
mcompshell11.test()


text:
#==============================================================================#
# names.jl
#
# AWS Endpoint URLs and Amazon Resource Names.
#
# http://docs.aws.amazon.com/general/latest/gr/rande.html
# http://docs.aws.amazon.com/general/latest/gr/aws-arns-and-namespaces.html
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

export aws_endpoint, arn, arn_region

#------------------------------------------------------------------------------#
# AWS Endpoint URLs
#
# e.g.
#
#   aws_endpoint("sqs", "eu-west-1")
#   "http://sqs.eu-west-1.amazonaws.com"
#
#------------------------------------------------------------------------------#

function aws_endpoint(service, region="", hostname_prefix="")

    protocol = "http"

    # HTTPS where required...
    if service in ["iam", "sts", "lambda", "apigateway", "email"]
        protocol = "https"
    end

    # Identity and Access Management API has no region suffix...
    if service == "iam"
        region = ""
    end

    # SES not available in all regions...
    if service == "ses" && !(region in ["us-east-1", "us-west-2", "eu-west-1"])
        region = "us-east-1"
    end

    # No region sufix for s3 or sdb in default region...
    if region == "us-east-1" && service in ["s3", "sdb"]
        region = ""
    end

    # Append region to service...
    if region != ""
        if service == "s3"
            service = "$service-$region"
        else
            service = "$service.$region"
        end
    end

    # Add optional hostname prefix (e.g. S3 Bucket Name)...
    if hostname_prefix != ""
        service = "$hostname_prefix.$service"
    end

    return "$protocol://$service.amazonaws.com"
end

#------------------------------------------------------------------------------#
# Amazon Resource Names
#------------------------------------------------------------------------------#

function arn(service, resource, region="", account="")

    if service == "s3"
        region = ""
        account = ""
    elseif service == "iam"
        region = ""
    end

    "arn:aws:$service:$region:$account:$resource"
end

function arn(aws::SymbolDict,
             service,
             resource,
             region=get(aws, :region, ""),
             account=aws_account_number(aws))

    arn(service, resource, region, account)
end

arn_region(arn) = split(arn, ":")[4]

#==============================================================================#
# End of file.
#==============================================================================#


text:
"""
    PrimeDivisors(n::T) where {T<:Integer}

Return the sorted array of prime divisors of `n`.
"""
function PrimeDivisors(n::T) where {T<:Integer}
    Vector{T}(GAP.Globals.PrimeDivisors(GAP.julia_to_gap(n)))
end;

function RootInt(n, k = 2)
    # check the arguments and handle trivial cases
    if k <= 0
      error("<k> must be positive")
    elseif k == 1
      return n
    elseif n < 0 && mod(k, 2) == 0
      error("<n> must be positive")
    elseif n < 0 && mod(k, 2) == 1
      return -RootInt(-n, k)
    elseif n == 0
      return 0
    elseif n <= k
      return 1
    end

    # r is the first approximation, s the second, we need: root <= s < r
    r = n
    s = BigInt(2)^(div(LogInt(n, 2), k)+1)-1

    # do Newton iterations until the approximations stop decreasing
    while s < r
      r = s
      t = r^(k-1)
      s = div(n+(k-1)*r*t, k*t)
    end

    # and that's the integer part of the root
    return r
end

"""
    SmallestRootInt(n::T) where {T<:Integer}

The smallest root of an integer `n` is the integer `r` of smallest absolute
value for which a positive integer `k` exists such that `n == r^k`.
"""
function SmallestRootInt(n::T) where {T<:Integer}
    # check the argument
    if n > 0
      k = 2
      s = 1
    elseif n < 0
      k = 3
      s = -1
      n = -n
    else
      return 0
    end

    # exclude small divisors, and thereby large exponents
    if mod(n, 2) == 0
      p = 2
    else
      p = 3
      while p < 100 && mod(n, p) != 0
        p = p+2
      end
    end
    l = LogInt(n, p)

    # loop over the possible prime divisors of exponents
    # use Euler's criterion to cast out impossible ones
    while k <= l
      q = 2*k+1
      while ! PrimesGAP.isprime(q)
        q = q+2*k
      end
      if powermod(n, div(q-1, k), q) <= 1
        r = RootInt(n, k)
        if r ^ k == n
          n = r;
          l = div(l, k)
        else
          k = PrimesGAP.nextprime(k+1)
        end
      else
        k = PrimesGAP.nextprime(k+1)
      end
    end

    return s * n
end

"""
    IsPrimePowerInt(n::T) where {T<:Integer}

Return `true` if `n` is a prime power, and `false` otherwise.
"""
IsPrimePowerInt(n::T) where {T<:Integer} = PrimesGAP.isprime(SmallestRootInt(n))

"""
    Lambda(m::T) where {T<:Integer}

Return the exponent of the group of prime residues modulo the integer `m`.
"""
function Lambda(m::T) where {T<:Integer}
    # make <m> it nonnegative, handle trivial cases
    if m < 0
      m = -m
    end
#   if m < Length(PrimeResiduesCache) then
#     return Length(PrimeResiduesCache[m+1]);
#   end
    # loop over all prime factors $p$ of $m$
    lambda = 1
    for p in PrimeDivisors(m)
        # compute $p^e$ and $k = (p-1) p^(e-1)$
        q = p
        k = p-1
        while mod(m, (q*p)) == 0
          q = q * p
          k = k * p
        end
        # multiples of 8 are special
        if mod(q, 8) == 0
          k = div(k, 2)
        end
        # combine with the value known so far
        lambda = lcm(lambda, k)
    end
    return lambda
end

"""
    LogInt(n::Int, base::Int)

Return the integer part of the logarithm of the positive integer `n`
with respect to the positive integer `base`,
i. e., the largest positive integer ``e`` such that ``base^e \\leq n``.
"""
LogInt(n::T, base::Int) where {T<:Integer} = Int(floor(log(base, n)))

"""
    OrderMod(n::Int, m::T1, bound::T2 = Lambda(m)) where {T1<:Integer} where {T2<=Integer}

Return the multiplicative order of the integer `n` modulo the
positive integer `m`, i. e.,
the smallest positive integer ``i`` such that ``n^i \\equiv 1 \\pmod{m}``.
The function returns `0` if `n` and `m` are not coprime.
"""
function OrderMod(n::Int, m::T1, bound::T2 = Lambda(m)) where {T1<:Integer} where {T2<:Integer}
    # check the arguments and reduce $n$ into the range $0..m-1$
    if m <= 0
      error("<m> must be positive")
    end
    if n < 0
      n = mod(n, m) + m
    end
    if m <= n
      n = mod(n, m)
    end

    # return 0 if $m$ is not coprime to $n$
    if gcd(m, n) != 1
      o = 0

    # compute the order simply by iterated multiplying, $x= n^o$ mod $m$
    elseif m < 100
      x = n
      o = 1
      while x > 1
        x = mod( x * n, m )
        o = o + 1
      end

    # otherwise try the divisors of $\lambda(m)$ and their divisors, etc.
    else
      o = bound
      for d in PrimeDivisors(o)
        while mod(o, d) == 0 && powermod(n, div(o, d), m) == 1
          o = div(o, d)
        end
      end
    end

    return Int(o)
end;

"""
    DuplicateFree(v)

Return a vector that contains the same entries as the array `v`,
but each entry occurs exactly once.
"""
function DuplicateFree(v)
  v1 = Vector{eltype(v)}()
  if length(v) > 0
    push!(v1, v[1])
    for elm in v
      if ! (elm in v1)
        push!(v1, elm)
      end
    end
  end
  return v1
end

# recursive function used by `DivisorsInt`
function DivisorsIntInner(i::Int, m::Int, factors::Vector{Int})
    if length(factors) < i
      return [m]
    elseif mod(m, factors[i]) == 0
      return DivisorsIntInner(i+1, m*factors[i], factors)
    else
      res = DivisorsIntInner(i+1, m, factors)
      append!(res, DivisorsIntInner(i+1, m*factors[i], factors))
      return res
    end
end

"""
    DivisorsInt(n::Int)

Return an array of all divisors of the integer `n`.
This array is sorted, so that it starts with `1` and ends with `n`.
We define that `DivisorsInt( -n ) == DivisorsInt( n )`.
"""
function DivisorsInt(n::Int)
    # make <n> nonnegative, handle trivial cases, and get prime factors
    if n < 0
      n = -n
    end
    if n == 0
      error("DivisorsInt: <n> must not be 0")
    end
#   if n <= length(DivisorsIntCache)
#     return DivisorsIntCache[n]
#   end
    factors = PrimesGAP.factors(n)

    return sort(DivisorsIntInner(1, 1, factors))
end

"""
    Phi(m::T) where {T<:Integer}

Return the value of Euler's totient function of `m`.
"""
function Phi(m::T) where {T<:Integer}
    # make <m> nonnegative, handle trivial cases
    if m < 0
      m = -m
    end

#   if m < Length(PrimeResiduesCache) then
#     return Length(PrimeResiduesCache[m+1]);
#   fi;

    # compute $phi$
    phi = m
    for p in PrimeDivisors(m)
      phi = div(phi, p) * (p-1)
    end

    # return the result
    return phi
end

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
mutable struct WeightMatrix
    seed::Int64
    log_prices::Array{Float64,1}
    block_size::Int64
    num_bootstap::Int64
end

function block_bootstrap_index(block_ind, n, b)
    rand_blocks = sample(block_ind, floor(Int,n/b))
    sample_ind = transpose(repeat(rand_blocks, 1, b))
    sample_ind = sample_ind[:]
    addition_vec = repeat(0:b-1,floor(Int,n/b))
    sample_ind = sample_ind + addition_vec
    return sample_ind
end

function block_bootstrap_estimator(wm::WeightMatrix)
    # Step 1: Apply a Moving Block Bootstrap to the Measured Series
    Random.seed!(wm.seed)

    n = size(wm.log_prices, 1)
    block_ind = 1:n-wm.block_size+1
    # b_samples = [log_prices[block_bootstrap_index(block_ind, n, b)] for i in 1:num_bootstap]
    b_samples = Array{Float64,2}(undef, n, wm.num_bootstap)
    for i in 1:wm.num_bootstap
        b_samples[:,i] = wm.log_prices[block_bootstrap_index(block_ind, n, wm.block_size)]
    end

    # Step 2: Calculate Distributions for Each Moment and Test Statistic
    dist = get_summary_stats(b_samples, wm.log_prices)
    W = inv(cov(dist))
    return W
end

function(wm::WeightMatrix)(algo)
    return algo(wm)
end

function select_moments(log_prices)
    mean_log = mean(log_prices)
    std_log = std(log_prices)
    kurt_log = normal_kurtosis(log_prices)
    ks_stat = 0.0
    hurst_log = generalized_hurts_exp(log_prices)

    return [mean_log, std_log, kurt_log, ks_stat, hurst_log]
end

# function identity_estimator(moments, simulated_moments)
#     return Matrix{Float64}(I, size(moments,1), size(moments,1))
# end
#
# function two_step_cov_estimator(moments, simulated_moments)
#     weight_matrix_1 = identity_estimator(moments, simulated_moments)
#     smm_estimated_parameters =
#     moment_error =
#     Ω = 1/N * (moment_error * moment_error')
#     weight_matrix_2 = inv(Ω)
#     return weight_matrix_2
# end
#
# function iterated_var_cov_estimator(moments, simulated_moments)
#
# end
#
# function newey_west_consistent_estimator(moments, simulated_moments)
#
# end